U(1) -symmetric infinite projected entangled-pair state study of the spin- 1/2 square J 1 − J 2 Heisenberg model
UU (1) -symmetric infinite projected entangled-pair state study of the spin- / square J − J Heisenberg model
R. Haghshenas and D. N. Sheng Department of Physics and Astronomy, California State University, Northridge, California 91330, USA (Dated: June 1, 2018)We develop an improved variant of U (1) -symmetric infinite projected entangled-pair state (iPEPS) ansatzto investigate the ground state phase diagram of the spin- / square J − J Heisenberg model. In order toimprove the accuracy of the ansatz, we discuss a simple strategy to select automatically relevant symmetricsectors and also introduce an optimization method to treat second-neighbor interactions more efficiently. Weshow that variational ground-state energies of the model obtained by the U (1) -symmetric iPEPS ansatz (fora fixed bond dimension D ) set a better upper bound, improving previous tensor-network-based results. Bystudying the finite- D scaling of the magnetically order parameter, we find a N´eel phase for J /J < . . For . < J /J < . , a non-magnetic columnar valence bond solid (VBS) state is established as observed bythe pattern of local bond energy. The divergent behavior of correlation length ξ ∼ D . and vanishing orderparameters are consistent with a deconfined N´eel-to-VBS transition at J c /J = 0 . , where estimatedcritical anomalous exponents are η s ∼ . and η d ∼ . for spin and dimer correlations respectively. We showthat the associated VBS order parameter monotonically increases with J /J and finally a first-order quantumphase transition takes place at J c /J = 0 . to the conventional Stripe phase. We compare our resultswith earlier DMRG and PEPS studies and suggest future directions for resolving remaining issues. PACS numbers: 75.40.Mg, 75.10.Jm, 75.10.Kt, 02.70.-c
I. INTRODUCTION
Understanding of quantum many-body systems is of funda-mental importance. These systems, even in the simplest form,reveal fascinating quantum collective behavior distinctly dif-ferent from noninteracting particles. For instance, frustratedquantum spin systems, defined by a simple spin model, areconsidered one of the most important playgrounds to observeexotic phenomena. Quantum spin liquid with a topolog-ically order, valence bond solid (VBS) order and de-confined quantum criticality are some of well-known ex-amples manifested in such systems. Specifically, searchingfor the quantum spin-liquid states has received much atten-tion due to their distinct characteristics, such as long-rangeentanglement and nontrivial anyon statistics. A compre-hensive characterization of them might lead to new under-standing in physics of frustrated quantum magnetism and pro-viding ‘a new theoretical framework’ for characterizing ex-otic phases of matter. J ∼ J The frustrated spin- / J − J square Heisen-berg model (SHM) is one of the candidate models fea-turing aforementioned exotic phases. The J − J SHMhas stimulated extensive theoretical studies over the lasttwo decades, due to its simplicity and its experimental re-alization in several materials, such as vanadium Lay-ered oxides Li VO ( Si,Ge ) O and polycrystalline samplesBaCdVO ( PO ) . In particular, these studies have establishedthat the second-neighbor J coupling controlling frustrationinduces non-magnetic phases around the highly frustratedpoint J /J = 0 . . Despite that, depending on numer-ical approaches, several scenarios have been proposed aroundthis point: the earlier studies based on small-size exact di-agonalization, spin-wave theory, series expansion and cou-pled cluster methods find different candidate states, such ascolumnar, plaquette VBS states and resonating va- lence bond spin liquid states.The recent SU (2) -symmetric density matrix renormaliza-tion group (DMRG) study has demonstrated an intermediateplaquette VBS phase between a N´eel and Stripe magneticallyordered phases, which does not support the previous DMRGresults of gapped Z spin liquid as the intermediate phase. However, in a small window of . < J /J ≤ . , the SU (2) DMRG results cannot distinguish between two pos-sible scenarios, between a true deconfined quantum-criticalpoint or a gapless spin-liquid phase. A very recent DMRGstudy further supports the possibility of a gapless spin liq-uid between the N´eel and the VBS phases by following theenergy level crossings between different low energy excitedstates. On the other hand, variational Monte Carlo (VMC)results predict a gapless Z spin liquid in the whole region . ≤ J /J ≤ . , while a very recent VMC study chal-lenged this result by predicting a columnar VBS order for . ≤ J /J ≤ . . The critical exponents reported in thisstudy show small deviation from those of the J - Q models.However, understanding the true nature of quantum criticalpoints and the corresponding universality classes turn out tobe even more challenging using unbiased methods. Recently, tensor-network-based methods have also been ap-plied to study the J − J SHM. An early plaquette renor-malized tensor-network study has predicted a possible pla-quette VBS order for the intermediate phase. They estimatedthe second-order phase transition between N´eel and plaque-tte VBS phase to occur around J c /J ≈ . . On the otherhand, finite-size projected entangled pair states (PEPS) ansatzwith the cluster-update scheme finds a direct N´eel-to-VBStransition occurring at J c /J ≈ . . The finite-size PEPSresults did not identify the true nature of VBS order, specifi-cally between plaquette and columnar. They also find corre-sponding critical exponents are consistent with the J - Q mod-els. A very recent SU (2) -symmetric infinite PEPS (iPEPS) a r X i v : . [ c ond - m a t . s t r- e l ] M a y ansatz suggests a quantum critical point at J c /J (cid:39) . ,where in contrast to the finite-size PEPS results, the ex-tracted critical exponents seem to deviate from those of the J - Q models. In this paper, we aim to develop a fully U (1) -symmetriciPEPS ansatz with an ‘improved’ update scheme to reexaminethe phase diagram of the J − J SHM. So far, the iPEPS up-date algorithms have been able to treat the first-neighborinteractions with high efficiency. They have been shown inpractice to be quite accurate and stable providing reliable re-sults. However, in the case of longer-range interactions (e.g.second-neighbor interactions) a similarly efficient scheme isstill highly desired. To this end, we present a new up-date method based on the so-called positive approximant andreduced-tensor application to treat second-nearest neigh-bor interactions more accurately and efficiently. We find thatthe new update scheme significantly improves efficiency andprovides more accurate results in comparison with previousschemes.
In addition, we also investigate the implementa-tion of U (1) symmetry into the iPEPS ansatz by introducinga general scheme to pick up relevant symmetry sectors. Weshow that it solves the loss of accuracy observed when ap-plying continuous symmetry groups and provides the sameaccuracy as non-symmetric iPEPS.By using the U (1) -symmetric iPEPS ansatz, we clarify thequantum phase diagram and the nature of phase transitionsfor the J - J SHM with substantially improved variationalwave function (of the ground state), and bridge the gap be-tween the previous tensor-network and DMRG studies. Weshow that the non-magnetic phase appears in the range of . < J /J ≤ . . The critical point J c /J (cid:39) . is of the deconfined type confirmed by continuously vanish-ing the N´eel order parameter and the divergence of the cor-relation length ξ ∼ D . . By extrapolating dimer-dimerand spin-spin correlation functions in the D → ∞ limit,we estimate the critical anomalous exponents η s ∼ . and η d ∼ . . The pattern of the local nearest neighboring bondenergies shows that a columnar VBS phase is established upto J c /J (cid:39) . . However, the observed (variational) ener-gies from different approaches indicate both columnar andplaquette VBS phases are competitive candidates for the in-termediate phase. With further increasing J /J , a first-orderphase transition takes place from VBS phase to the conven-tional Stripe phase.The paper is organized as follows. We first introduce themodel and briefly summarize different types of the phases andthe resulting phase diagram obtained by our iPEPS studies inSec. II. In Sec. III, we briefly introduce the U (1) -symmetriciPEPS ansatz and discuss a general scheme to select automat-ically relevant symmetric sectors (Sec. III B). We then presenta new iterative scheme in detail and compare it with previousschemes (Sec. III C). Sec. IV provides the main simulation re-sults. The variational ground-state energy and N´eel order pa-rameter are presented in Secs. IV C. We show that the interme-diate phase is a columnar VBS represented in Sec. IV D. Thecritical properties of the deconfined quantum-critical pointare discussed by studying correlation function and correla-tion length in Sec. IV E—further plots of the correlation func- FIG. 1. (Color online) Phase diagram of the J − J SHM as afunction of coupling J . The arrows show pattern of magnetic orderappeared in AFM N´eel and Stripe phases. The eclipses in intermedi-ate phase (columnar VBS) stand for entangled spins (singlet states). tions are presented in Appendix. A. Using different initial ten-sors representing different symmetry breaking states, we de-termine the boundary of columnar VBS and the conventionalStripe phase in Sec. IV F. Finally, we summarize our workwith some discussions in Sec. V. II. MODEL
The J − J SHM is defined by the Hamiltonian H = J (cid:88) (cid:104) i,j (cid:105) S i · S j + J (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) S i · S j , where S i ≡ ( S xi , S yi , S zi ) are spin- / operators. The cou-plings J and J stand for the first- and second-neighbor anti-ferromagnetic (AFM) interactions. We set J = 1 throughoutthe paper and consider the frustrated interaction J > .In the extreme cases J ≈ or J (cid:29) , the ground statesare respectively defined by two magnetically ordered phases,i.e., AFM N´eel and Stripe. The patterns of magnetic orders forthese phases have been shown in Fig. 1. All the earlier studiessuggest that these two phases are separated by an (or several)intermediate phase(s). Our goal is to locate and characterizethe intermediate phase.The obtained phase diagram has been illustrated in Fig. 1.We find that the intermediate phase is a paramagnetic phasethat breaks lattice symmetry, i.e., a columnar VBS. As seenin Fig. 1, columnar VBS order (in which vertical spins arestrongly entangled) only breaks lattice symmetry in the y -direction. The columnar VBS phase is separated from theN´eel one by a continuous phase transition occurred at J c =0 . . In addition, the quantum phase transition betweenVBS and AFM Stripe phases takes place at J c = 0 . ,which is of the first-order type. left do w nup right FIG. 2. (Color online) Tensor-network representation of the iPEPSansatz. ( a ) U (1) -invariant five-rank tensor (particle numbers asso-ciated to incoming and outgoing arrows are equal). Virtual bondsare labeled by { left , down , right , up } . ( b ) A virtual bond is labeledby vectors ( −→ n , −→ t n ) = ( · · · , ( n ( i ) , t ( i ) n ) , · · · ) , where i th components n ( i ) and t ( i ) n represent a particle number and its associated dimen-sion. ( c ) U (1) -symmetric × unit cell iPEPS. III. METHODA. U (1) -symmetric iPEPS ansatz An iPEPS is constructed by building-block tensors that aresitting on sites of the physical lattice. The tensors are con-nected to each other by the so-called virtual bonds (graph-ically drawn by arrows) constructing a specific geometricalpattern (usually similar to physical lattice). For instance, asdepicted in Fig.2-(a-c), we have constructed a × unit celliPEPS on the infinite two-dimensional square lattice by re-peating periodically five-rank tensors { a, b, c, d } . The geo-metrical structure produced by the connections of tensors hasan important feature: the iPEPS could reproduce entangle-ment area law. The amount of this entanglement is con-trolled by the so-called bond dimension (number of elements)of virtual bonds, denoted by D . By increasing D , the iPEPSis able to represent highly entangled stats.We aim to use the iPEPS as a variational ansatz to obtainthe approximate ground state of the model. The accuracyof this variational method is controlled by the bond dimen-sion D (variational parameters are of order O ( D ) ). To cap-ture the physics of highly entangled states, one needs to con-sider larger D and does the finite- D analysis (extrapolating D → ∞ ). By exploiting U (1) symmetry, one can studythe iPEPS with a larger D . In the presence of this symme-try, each tensor takes a block diagonalized form (each blockis corresponding to a specific symmetric sector) which cor-respondingly reduces computational costs. The symmetricsectors in the case of U (1) symmetry are labeled by the con-served particle numbers −→ n . To implement this symmetry intothe iPEPS, we label each arrow by some particle numbers ( −→ n , −→ t n ) = ( · · · , ( n ( i ) , t ( i ) n ) , · · · ) , as depicted in Fig.2-(b), sothat their sign ( ± ) being specified by outgoing and incomingarrows. For example, a virtual bond with the associated label ( −→ n , −→ t n ) = (( − , , (1 , would have the bond dimension D = 5 with particle numbers − , and associated dimen- -3 -2 ∆ E DNo symmetryU(1)Bauer et al. -4 -3 -2 D U(1)No-symmetry
FIG. 3. (Color online) Relative error ∆ E in the ground-state energyof Heisenberg model on square lattice ( J = 0 ) for ( a ) simple- and ( b ) full-update schemes. The data shown by blue circles are obtainedby homogeneous U (1) -symmetric iPEPS (all virtual bond have thesame symmetry sector) reported in Ref. 46. The green square sym-bols show our results (a non-homogeneous structure, see for exampleTab. I) obtained by the scheme explained in the main text. sions , respectively. Each tensor is U (1) invariant as thesum of incoming particle numbers equals to the sum of out-going ones. In this case, when each individual tensor is U (1) invariant, the iPEPS automatically respects U (1) symmetry. B. Selection of relevant symmetric sectors
For infinite symmetric groups, like U (1) , it is not possibleto uniquely specify symmetric sectors −→ n . In addition, numberof states in each sector −→ t n should be manually chosen. Fur-thermore, the virtual bonds could possess non-homogeneousstructures: each virtual bond takes different symmetric sec-tors from another one. These possibilities in selecting thesymmetric sectors impede the U (1) -symmetric iPEPS ansatzfrom providing accurate results as reported in Ref. 46—note,this loss of accuracy is not observed in the case of finite sym-metric groups (even for homogeneous-bond structure) due tothe finite number of the symmetric sectors.We introduce a simple strategy to select automati-cally relevant symmetric sectors ( −→ n , −→ t n ) by using simple-update simulation. In this scheme, we assume a non-homogeneous structure for the virtual bonds (each of themcould take different symmetric sectors). We then performsimple-update simulation: ( i ) randomly initialize iPEPS(picking up a random set of symmetric sectors), ( ii ) applythe local imaginary time-evolution operator to the virtualbonds and ( iii ) use high-order singular value decomposition to keep the largest singular values, which also determines thesymmetric sectors. In addition, to obtain the expectation val-ues, we similarly assume a non-homogeneous structure forthe so-called environment tensors and pick up the sym-metric sectors by using the singular values decomposition ap-peared in the corner transfer matrix (CTM) approach. Fur-thermore, to do the full-update simulation, we first fix thesymmetric sectors for all virtual bonds (obtained by the afore- tensor left down right up a ( − , ,
2) ( − , − ,
1) ( − , ,
3) ( − , − , b ( − , ,
3) ( − , ,
3) ( − , ,
2) ( − , , c ( − , ,
3) ( − , − ,
1) ( − , ,
2) ( − , − , d ( − , ,
2) ( − , ,
3) ( − , ,
3) ( − , , TABLE I. The resulting particle numbers −→ n for the virtual bonds oftensors { a, b, c, d } for the Heisenberg model J = 0 . The associateddegeneracy is always −→ t n = (2 , , for all virtual bonds. The labels { left , down , right , up } show four virtual bonds of each tensor. Theparticle numbers have obtained by the scheme explained in the maintext. mentioned scheme), then randomly initialize tensors, and fi-nally apply the optimization schemes (e.g., as explained inSec. III C).We apply this scheme to the Heisenberg model on thesquare lattice ( J = 0 ) to compare its accuracy with thatof previous ones presented in Ref. 46. We use the rela-tive error ∆ E in the ground-state energy to provide bench-marks: ∆ E = E D − E MC E MC , where E D and E MC are respec-tively the iPEPS energy with finite bond-dimension D andthe precise Monte-Carlo energy from Ref. 57. As shown inFig. 3-(a, b), a proper choice of symmetric sectors makesthe U (1) iPEPS highly accurate. We observe that a U (1) -symmetric iPEPS ansatz produces the same accuracy as thenon-symmetric ones in both full- and simple-update methodsfor the same bond dimension D . The resulting symmetry sec-tors for virtual bonds (with D = 7 ) are shown in Tab. I, wherewe have started from a homogeneous structure. The particlenumbers −→ n dynamically vary during the simulation and (fi-nally) take a non-homogeneous structure—as −→ n is differentfor each virtual bond. Specifically, our results imply that thenon-homogeneous structures are crucial in the case of U (1) symmetry. C. Optimization method
In order to do a full-update simulation for the models in-cluding second-neighbor interactions, we introduce a new it-erative scheme to optimize the tensors. We use both applica-tions of positive approximant and reduced tensors in aniterative way to improve accuracy and convergence rate of theoptimization algorithm. We discuss the general ideas here,but for computational implementation and more details seeRefs. 47, 48, and 58.To perform a full-update simulation, we need tostudy imaginary-time evolution of an initial (random)iPEPS | ψ ( a, b, c, d ) (cid:105) . We use first-order Suzuki-Trotterdecomposition to split the imaginary time-evolution oper-ator into a sequence of local terms. Such local operators areacting on specific bonds, increasing the corresponding bonddimension D → D (cid:48) . In order to have a tractable algorithm,we need to reduce the bond dimension (approximating the re-sulting iPEPS). We explain this procedure by considering lo- FIG. 4. (Color online) Tensor-network diagram of reduced-tensorapplication. Tensors { b (cid:48) , c (cid:48) } are decomposed to low-rank tensors { l, r, Q, ¯ Q } by the means of LQ and QR decomposition. For exam-ple in ( a ) , tensors L and Q are obtained by fusing { left , up , down } indices and { right , physical } indices of tensor b (cid:48) (to make a matrix)and performing LQ decomposition of that matrix. cal imaginary-time operators acting on, e.g., tensors { a, b, c } U = e − δ ( J S c · S a + J S c · S b + J S a · S b ) , where δ stand for small time steps and S a is acting on tensor a (analogous for other operators). After applying U , the result-ing wave function should be approximated by a new iPEPSwith the bond dimension D | ψ (cid:48) ( a (cid:48) , b (cid:48) , c (cid:48) , d (cid:48) ) (cid:105) ≈ U | ψ ( a, b, c, d ) (cid:105) , where tensors { a (cid:48) , b (cid:48) , c (cid:48) , d (cid:48) } are determined by minimizing thesquare distance. We consider tensors { a (cid:48) , b (cid:48) , c (cid:48) , d (cid:48) } as varia-tional parameters and accordingly find them by minimizingthe square distance min { a (cid:48) ,b (cid:48) ,c (cid:48) ,d (cid:48) } f ( | ψ (cid:48) (( a (cid:48) , b (cid:48) , c (cid:48) , d (cid:48) )) (cid:105) , U | ψ ( a, b, c, d ) (cid:105) ) , where f = (cid:104) ψ | U † U | ψ (cid:105) + (cid:104) ψ (cid:48) | ψ (cid:48) (cid:105) − (cid:104) ψ (cid:48) | U | ψ (cid:105) − (cid:104) ψ | U † | ψ (cid:48) (cid:105) . To minimize the cost function f , we use positive approximantand reduced-tensor schemes:(a) reduced-tensor application: We use QR and LQ de-composition to split tensors { b (cid:48) , c (cid:48) } to sub-tensors { l, r, Q, ¯ Q } as depicted in Fig. 4-(a). We aim to mini-mize the cost function with respect to tensors { l, a (cid:48) , r } ,thus we rewrite the cost function as following min { a (cid:48) ,r,l } f = const + r † a (cid:48)† l † N la (cid:48) r − r † a (cid:48)† l † ¯ N − ¯ N † ra (cid:48) l, where N is called norm tensor. Tensor-network repre-sentations of the norm tensor and the term r † a (cid:48)† l † N la (cid:48) r are shown in Fig. 5-(a, b). Note that the first term doesnot play any role in the optimization procedure. = FIG. 5. (Color online) . ( a ) Tensor-network representation of normtensor N . Tensors { A , · · · , A , B , · · · , B , T , · · · , T } are theenvironment tensors (with bond dimenstion χ ) obtained by CTMrenormalization group approch. ( b ) tensor-network representa-tion of the term r † a (cid:48)† l † N la (cid:48) r appeared in cost function f . Note thatby taking conjugate transpose the direction of arrows changes. ( c ) Tensor-network representation of transfer matrix. For instance, ten-sor A is obtained by contracting physical index of tensors a † and a ,and then fusing (combining) corresponding virtual bonds—so, thebond dimension of each virtual bond is D . (b) positive approximant: In principle, the norm tensor N should be positive and Hermitian. But mainly due to theCTM approximation, it has some small negative parts.We explicitly eliminate that part by enforcing N to bepositive. We also replace N by its Hermitian positivecounterpart ( N + ) in the cost function: N + = (cid:112) N ,where N = ( N + N † ) / .(c) alternating-least-squares (ALS) sweep: We then itera-tively optimize the cost function by finding the optimumtensors { l, a (cid:48) , r } : e.g., we minimize the cost functionwith respect to a (cid:48) by solving equation ∂ a (cid:48)† f = 0 byholding fixed tensors l, r . Then we repeat this proce-dure for another tensor with holding rest fixed until costfunction converges.(d) recovering: After finding optimum tensors { l, a (cid:48) , r } , weabsorb tensor { l, r } to { Q, ¯ Q } to recover the final opti-mum tensors { b (cid:48) , a (cid:48) , c (cid:48) } .The positive approximant in our scheme is crucial in mak-ing the algorithm highly stable and accelerating its conver-gence. The computational cost of the norm tensor and ALSsweep are respectively O ( D χ ) and O ( D ) , where χ is thebond dimension of the environment tensors. Since we onlyneed to calculate the norm tensor once, the dominant compu-tational cost belongs to ALS sweep, i.e., O ( D ) . We should also notice in the case that χ > D , that is suppressed by O ( D χ ) (as occurs in our calculations).Although, in practice, the steps-( a - d ) provide proper ac-curacy and approximates the iPEPS wave function well, butthere is still room to improve it. Specifically, we did not takeinto account tensor d (cid:48) in the optimization procedure (left un-touched). In addition, sub-tensor application might reduce theaccuracy. The main idea is to use sub-tensor application ina different way (see Fig. 4-(b)) to design an efficient strat-egy to include tensor d (cid:48) in the optimization procedure. Thesteps are as follows: ( i ) we decompose tensors { b (cid:48) , c (cid:48) } to sub-tensors { l, r, Q, ¯ Q } as shown in Fig. 4-(b) and rewrite the costfunction f = const + r † d (cid:48)† l † N ld (cid:48) r − r † d (cid:48)† l † ¯ N − ¯ N † rd (cid:48) l , ( ii ) use positive approximation N → N + , ( iii ) optimizethe cost function by finding the optimum tensors { l, d (cid:48) , r } and ( iv ) after finding optimum tensors { l, d (cid:48) , r } , we absorbtensor { l, r } to { Q, ¯ Q } to recover the final optimum tensors { b (cid:48) , d (cid:48) , c (cid:48) } . The computational costs for the steps-( i - iv ) aresimilarly O ( D χ ) and O ( D ) .The optimization procedure is completed by iteratively re-peating steps-( a - d ) and -( i - iv ) until the cost function does notchange up to the desired threshold. In Fig. 6-(a, b) we haveplotted the typical behavior of the cost function f and its meanvalue of the relative change ¯ f = | f u +1 − f u f init | versus consec-utive iteration number u for different optimization schemes.It is seen that our scheme significantly improves convergencerate and provides better accuracy than previous schemes. In the full CG method, all tensors { a (cid:48) , b (cid:48) , c (cid:48) , d (cid:48) } are entirelyoptimized which makes its final result highly accurate. Weempirically observe that CG method eventually provides bet-ter accuracy than our scheme after ∼ iterations. In Fig. 6-(c), we have plotted the effect of these optimization schemeson the ground-state energy. The ground-state energy is calcu-lated at J = 0 . by using U (1) -symmetric iPEPS with bonddimension D = 6 . Similarly, it shows that our scheme im-proves the ground-state energy as expected. -0.4945-0.494-0.4935 0 5 10 15 20 25 30 35 E iterationfull CGsteps-(a-iv)Corboz et al -5 -4 -3 u full CGsteps-(a-iv)steps-(a-d)Corboz et al -8 -7 -6 -5 -4 -3 -2 u steps-(a-d)Corboz et alsteps-(a-iv) FIG. 6. (Color online) ( a, b ) Log-linear plots of the cost function f and its mean value of the relative change ¯ f as a function of con-secutive iterations u for different optimization schemes with bonddimension D = 6 at J = 0 . . The data were generated for a fixedtime step δ = 0 . . ( c ) The ground-state energy versus loop itera-tions for the optimization schemes. The data shown by red triangularare obtained by the optimization method introduced in Ref. 45.
IV. RESULTSA. Simulation remarks
In our simulation, we run several full-update simulationsinitialized by random and/or ordered states (such as N´eel,VBS and Stripe) to find the lowest variational ground-state en-ergy. We first pick up the symmetric sectors with the schemeexplained in Sec. III and then start the optimization procedureby performing the iterative scheme. At the end, a few stepsof the full CG method is used to improve the results evenmore. All data points reported here correspond to the low-est variational ground-state energy that we have been able toobtain. The largest bond dimensions that we could afford are ( D, χ ) = { (8 , , (9 , } .We always check the behavior of the ground-state energywith respect to χ to make sure that the error due to the en-vironment approximation is negligible. The expectation val-ues are calculated by a modified CTM renormalization groupapproach. We find that this approach produces much bet-ter convergence rate and more accurate results in comparisonwith other variants of CTM ones.
B. Order parameters
We need to define some order parameters to establish dif-ferent ordered phases appeared in the J − J SHM. Mag-netically ordered phases could be addressed by magnetizationparameter m = (cid:80) i |(cid:104) S zi (cid:105)| , where the index i runs over the sitescorresponding to the building-block tensors { a, b, c, d } . In ad-dition, we use the local the nearest neighboring bond energyto detect the translational lattice symmetry breaking. To thisend, we define the following order parameters ∆ T x = max( E x ) − min( E x ) , ∆ T y = max( E y ) − min( E y ) , where E x and E y stand for local nearest neighboring bondenergy in the x - and y -directions. The order parameters ∆ T x and ∆ T y stand for different type of lattice symmetrybreaking. The lattice order parameters plus the magnetiza-tion are capable of distinguishing between the phases ap-peared in the J − J SHM. In the N´eel phase, we expect { m (cid:54) = 0 , ∆ T x = 0 , ∆ T y = 0 } as the local bond energy re-mains the same in different directions. In the columnar VBSphase orientated in y -direction (analogous to one in Fig. 1),we expect { m = 0 , ∆ T x = 0 , ∆ T y (cid:54) = 0 } , while in AFMStripe phase it becomes { m (cid:54) = 0 , ∆ T x = 0 , ∆ T y = 0 } (seeFig. 1). Thus, by studying m and ∆ T y , we are able to detectN´eel-to-VBS and VBS-to-Stripe quantum phase transitions.In order to study quantum critical points, we use (con-nected) transverse correlation function defined by C t ( r ) = (cid:104) (cid:98) O ( x,y ) (cid:98) O ( x + r,y ) (cid:105) − (cid:104) (cid:98) O ( x,y ) (cid:105)(cid:104) (cid:98) O ( x + r,y ) (cid:105) , where indices ( x, y ) show spatial coordinate and subindex t stand for word ‘transverse’. The operators (cid:98) O ( x,y ) are chosento be S ( x,y ) and S ( x,y ) · S ( x,y +1) , respectively, for the spin-spin ( C st ) and dimer-dimer ( C dt ) correlation functions. The corre-lation function could determine universality class of a criticalphase, revealed in the power-law behavior. It algebraicallyfalls off at critical point as C st ∼ r − (1+ η s ) ,C dt ∼ r − (1+ η d ) , where { η s , η d } are anomalous spin and dimer exponents, re-spectively. A finite bond dimension D (usually) induces ex-ponential decay ( C t ∼ e − rξ ) for large distances r (cid:29) . Thus,we find the correlation length ξ by obtaining the slopes of thefollowing function log( C t ( r )) = ( − ξ ) r + const r (cid:29) . and obtain the critical behavior through scaling to large bonddimension limit. Note that in this method, there is one as-sociated correlation length for each correlation function. Wecould also define characteristic correlation length ξ by usingthe eigenvalues of the transfer matrix as shown in Fig. 5-(c).It is given by ξ = − | λ λ | ) where λ and λ ( | λ | > | λ | ) arerespectively the largest eigenvalues of the transfer matrix.In some cases, when the system reveals different correlationlengths in x - and y -directions, we also need to define longitu-dinal correlation functions given by C l ( r ) = (cid:104) (cid:98) O ( x,y ) (cid:98) O ( x,y + r ) (cid:105) − (cid:104) (cid:98) O ( x,y ) (cid:105)(cid:104) (cid:98) O ( x,y + r ) (cid:105) . where subindex l stand for word ‘longitudinal’. C. N´eel phase
In Fig. 7-(a), we have compared the ground-state energy ob-tained by U (1) -symmetric iPEPS with the extrapolated value( N → ∞ ) of the finite-size PEPS at J = 0 . . Ourground-state energies for bond dimensions D ≥ are lowerthan that of the finite-size PEPS with a larger bond dimension D = 9 . The PEPS results are obtained by the cluster-updatescheme which is considered as an intermediate optimiza-tion approach between simple and full update—it is computa-tionally cheaper than full update. Our best variational energy E D =8 iPEPS = − . sets an upper bound to the true ground-state energy at J = 0 . .We have also compared the result of the SU (2) -symmetriciPEPS ansatz applied recently to the J − J SHM at J =0 . . As seen in Fig. 7-(b), U (1) iPEPS with bond dimensions ( D, χ ) = (5 , provides the same ground-state energy as SU (2) iPEPS with bond dimensions ( D, χ ) = (7 , χ → ∞ ) .The reason might be due to the effect of finite bond dimension;as the SU (2) iPEPS provides an efficient representation foronly symmetric phases. We argue that the system at J = 0 . is still magnetically ordered, thus, the SU (2) iPEPS probablypicks up a superposition of the states requiring larger bond di-mensions. Our best variational energy at the highly frustratedpoint J = 0 . is E D =9 iPEPS = − . , which is quite close tothe DMRG extrapolated value, E DMRG = − . . -0.5097-0.5091-0.5084-0.505 1/8 1/5 1/3 E -0.491-0.49-0.488-0.486-0.484 1/8 1/5 1/3 E -0.498-0.4968-0.495-0.4886 1/9 1/5 1/3 m J =0.0J =0.45J =0.50J =0.53 FIG. 7. (Color online) ( a, b ) A comparison between the U (1) -symmetric iPEPS variational energy with that of the finite-size PEPS, DMRG and SU (2) -symmetric iPEPS at J = { . , . } , respectively. The blue (purple) dashed line repre-sent polynomial (linear) fit up to the fourth order. ( c ) The U (1) -symmetric iPEPS variational energy at the deconfined quantum-critical point J c = 0 . . ( d ) The N´eel order parameter m as afunction of /D . A linear fit in large-D limit reveals m vanishes atpoint J = 0 . . We study the N´eel order parameter m as a function of thebond dimension D to find the critical point, where the N´eelphase disappears. At point J = 0 , we find that a linearextrapolation with the large bond dimensions ( D ≥ ) pro-vides a proper estimation of m . The relative error of ourestmation with that of the Monte-Carlo result is of order ∆ m = m D →∞ − m MC m MC < × − . In Fig. 7-(d), we have plot-ted the N´eel order parameter m versus /D for different val-ues of J . A linear extrapolation (dashed lines) for the largerbond dimensions ( D ≥ ) reveals that m remains finite in therange of ≤ J ≤ . . In this interval, the order parame-ters ∆ T x < . and ∆ T y < . are both small consistentwith that the N´eel phase persists up to point J = 0 . . Atthis point, m is almost zero ( < − ), thus, we conclude thequantum critical point occurs at J c = 0 . . The ground-state energy at this point has been shown in Fig. 7-(c): thebest upper bond on the ground-state energy and the extrapo-lated value are E D =8 iPEPS = − . and E D →∞ iPEPS = − . (from polynomial fit), respectively. D. Columnar VBS phase
We study the order parameters ∆ T x , ∆ T y and correlationfunctions for J > . to find the true nature of the non- magnetic phase. We plot ∆ T x and ∆ T y for points J = { . , . , . } as depicted in Fig. 8-(a). It suggests a colum-nar VBS order for no-magnetic phase: in the large-D limit, ∆ T y remains finite, while ∆ T x is one order of magnitudesmaller than ∆ T y . By increasing J , the order parameter ∆ T y monotonically increases and reaches its maximum valuearound J ≈ . —we will show later that a first-order phasetransition takes place at this point. To check the validity ofthe result, we compare the ground-state energy with previousstudies at J = 0 . , as depicted in Fig. 8-(b). At this point,DMRG and finite-size PEPS study respectively predicteda plaquette VBS order and a critical behavior (algebraic fall-off of the correlation function up to N = 24 × ). We expectthat this critical behavior, in the PEPS calculation, eventuallydisappears in the thermodynamic limit. We notice that the es-sential difference between the iPEPS and PEPS anstaz lies inthe finite-size boundary effects, as both are using the same un-derlying tensor-network wave function. Since our variationalenergy is quite compatible with that of finite-size PEPS andour result predicts a VBS order, we might conclude that asthe system size increases, algebraic fall-off of the correlationfunction get eventually dominated by an exponential behavior.In order to gain more insight, we investigate the correla-tion functions at the point J = 0 . . In Fig. 8-(c), we haveplotted the the transverse and longitudinal correlation func-tions. The longitudinal and transverse correlation functions Δ T y J2=0.60J2=0.55J2=0.54 Δ T x
123 4 5 6 7 8 ξ D transfer-matrix ξ y transfer-matrix, ξ x -0.481-0.479-0.477-0.475 1/16 1/9 1/6 E y DMRGCluster PEPS -12 -10 -8 -6 -4 -2 C s , C d rC s C s C d C d FIG. 8. (Color online) ( a ) Order parameters ∆ T y and ∆ T x (inset)as a function of bond dimension /D . ( b ) The variational ground-state energies at J = 0 . . The DMRG and the finite-size PEPS dataare respectively obtained on the tilted cylinder with width L y and L y × L y torus. ( c ) Log-log plot of the spin-spin and dimer-dimer cor-relation functions versus distance r in the x - and y -directions. Thedate has been plotted for bond dimension D = 8 , ( d ) The correlationlength versus bond dimension D . The dashed lines are power-law fits ξ ∼ D α . ξ D spin-spindimer-dimertransfer-matrix -8 -7 -6 -5 -4 -3 -2 -1 C s rD=4D=5D=6D=7D=8D → ∞ r -1.6 -9 -8 -7 -6 -5 -4 -3 -2 C d rD=5D=6D=7D=8D → ∞ r -2.9 -12 -10 -8 -6 -4 -2
105 15 20 25 30 35 40 C s rD=4D=5D=6D=7D=8 FIG. 9. (Color online) Behavior of the correlation functions at de-confined quantum-critical point . ( a, b ) Log-log plot of the spin-spinand dimer-dimer correlation functions versus distance r . ( c ) Log-linear plots of the spin-spin correlation function versus distance r ,where slopes show inverse of the correlation length ξ − s . ( d ) Thecorrelation lengths as a function of bond dimension D ; dashed linesare power-law fits ξ ∼ D α . show different correlation lengths as expected from the na-ture of the columnar VBS ordered state. We observe that theyexponentially fall off, as confirmed by the behavior of cor-relation lengths, shown in Fig. 8-(d). The characteristic cor-relation lengths ξ x,y increase slowly with bond dimension D and seem to saturate in the large- D limit. A power-law fit ξ ∼ D α to the largest bond dimensions D = { , , } reveals α < . . E. Deconfined quantum criticality
In this section, we investigate critical properties of the de-confined quantum-critical point by studying the correlationfunctions and the associated correlation lengths. We studycorrelation functions at the critical point J = 0 . and com-pare the results with the previous studies. In Fig. 9-(a, b),we have plotted correlation functions C st and C dt as a func-tion of distance r . The data for each bond dimension D areobtained by the largest environment bond dimension χ , al-though in contrast to Ref. 41 we do not observe any strongdependency on χ . In order to understand the true behav-ior of the correlation functions, we need to study the asso-ciated correlation lengths as a function of D . In Fig. 9-(c), weshow the log-linear plot of the spin-spin correlation functionversus large distance r (cid:29) . The slopes reveal the inverseof the spin correlation length ξ − s . ξ s increases significantlyby increasing the bond-dimension D as expected in a critical regime. They follow an empirical power-law relation ξ ∼ D α as shown in Fig. 9-(d). The spin correlation length ξ s andcharacteristic correlation length ξ (extracted from the trans-fer matrix) diverges similarly as ξ s , ξ ∼ D . . Instead, dimercorrelation length is governed by different scaling exponent as ξ d ∼ D . . This divergent behaviors suggest that J ≈ . is a critical point which is consistant with vanishing the orderparameters. The divergent behavior of correlation length ξ ( D ) impliesan algebraic fall-off of the correlation function in the rangeof < r < ξ ( D ) . An accurate estimation of the criticalanomalous exponents requires ξ ( D ) to be large enough. Par-ticularly, in our case, ξ ( D ) is still small even for the largestbond dimension. Thus, we need to rely on the extrapolateddata in the D → ∞ limit, which correspondingly represent alarge correlation length ξ ( D → ∞ ) >> . We use a linearextrapolation D > to obtain the extrapolated data of corre-lation function up to r ∼ , where error-bars are still small(see Appendix. A). As shown in in Fig. 9-(a, b), we have fit-ted the data (in D → ∞ ) to a power-law function r − (1+ η ) toestimate the exponents. The critical exponents for spin-spinand dimer-dime correlations are, respectively, η s ∼ . and η d ∼ . , which are in agreement with Ref. 41. Our resultsshow that dimer-dimer correlation falls off more rapidly thanpredicted by J - Q models ( . < η J - Q < . ). Thismay indicate different universality classes of the deconfinedcriticality for different models.Therefore, our results predict a continuous N´eel-to-VBStransition, which is forbidden in Landau-Ginzburg theory dueto the different types of broken symmetry—unless it wouldbe of the first-order type. So, we conclude that this quantumphase transition fits well in the paradigm of ‘deconfined quan-tum criticality’. However, the field-theory description of thisdeconfined quantum critical point might be different from thatof the J - Q models, as seen by different scaling behavior ofcorrelation functions. F. First-order quantum phase transition
We expect a quantum phase transition to occur betweencolumnar VBS and AFM Stripe phases as J increases. Tolocate the quantum phase transition point, we sketch dia-grammatically local nearest neighboring ( J ) bond energy at J = { . , . } , see Fig. 10-(a, b). The pattern of bondenergy shows the lattice symmetry breaking in the y -direction( ∆ T y ∼ . and ∆ T x ∼ . ) disappears at J = 0 . , wherethe AFM Stripe phase emerges. Since order parameter ∆ T y has been monotonically increased from the point J = 0 . ,we expect the quantum phase transition to be the first-ordertype rather than continuous one.We use hysteresis analysis as explained in Ref. 45 to findwhether the quantum phase transition is the first-order type: ( i ) we initialize the iPEPS ansatz by competitive orderedstates in the vicinity of the critical point, ( ii ) find where theenergies become equal for different values of bond dimension D and ( iii ) check whether the order parameters remain non-zero. As shown in Fig. 10-(c, d, e), we have compared the en-ergies internalized by columnar VBS and AFM Stripe statesat J = { . , . , . } . We observe that energies of stateswith different initializations at the point J = 0 . becomealmost equal, but for J = { . , . } , the columnar VBSand the AFM Stripe respectively provide lower energy. Thus,we conclude that they cross around J = 0 . . At this point,the order parameters for both states remain finite, as shown inin Fig. 10-(f). For columnar VBS and AFM Stripe, we respec-tively obtain in the large- D limit (∆ T y , ∆ T x ) ≈ (0 . , . and m ≈ . . Therefore, the transition occurring J c =0 . is of the first order. V. DISCUSSION AND CONCLUSION
In this paper, we have addressed two main obstacles re-garding the iPEPS ansatz: ( i ) how to improve iPEPS updateschemes in the presence of second-neighbor interactions and ( ii ) how to automatically select relevant symmetry sec-tors (in the case of continuous symmetry) without losingaccuracy. We considered the first issue by introducing an‘improved’ update scheme based on positive approximant andreduced-tensor application. The update scheme significantlyaccelerates the convergence rate and also improves the accu-racy/stability in comparison with previous schemes.
Forthe second issue, a simple strategy is introduced to pick uprelevant symmetry sectors so that the accuracy remains thesame as non-symmetric cases. We also showed that taking anon-homogeneous structure for all virtual bonds is crucial inthe case of the U (1) -symmetric iPEPS ansatz—which doesnot seem to be the case for finite symmetry groups.We utilize our U (1) -symmetric iPEPS ansatz to investigatethe ground-state phase diagram of the J − J SHM on thesquare lattice. A N´eel phase is found for J < . by ob-serving a non-zero value of the magnetically order parame-ter in the large-D limit. In the range . < J < . ,by studying the lattice symmetry breaking order parameters,we find that a columnar VBS phase is established. The point J c = 0 . represents a deconfined N´eel-VBS quantum criti-cal point, as confirmed by vanishing the order parameters anddivergent behavior of the characteristic correlation length andspin correlation length, i.e., ξ ∼ D . . This result is consis-tent with that of DMRG studies: accurate SU (2) -symmetricDMRG estimates the transition point ≈ > . , while a veryrecent U (1) -symmetric DMRG study based on level spec-troscopy has predicted the transition to be ≈ . , although asmall window of possible gapless spin-liquid is suggested inthis work. Our findings improve the result of finite-size PEPSstudy which obtained a critical point around ≈ . . Themain reason for such difference may come from the lack ofthe finite- D extrapolation in Ref. 40.We have studied dimer-dimer and spin-spin correlationfunctions to compare the associated critical exponents withthat of the J - Q model, i.e. η J - Q ∼ . . Our estimateddimer and spin anomalous exponents, η s ∼ . and η d ∼ . ,show deviation from that value. That observation is also man-ifested in the divergent behavior of the correlation lengths:the spin and dimer correlation lengths diverge as ξ s ∼ D . - . -0.16 -0.16 -016-0.16 -0.16 -0.16-0.16 -0.16 -0.16-0.16 -0.16 -0.16 - . - . - . - . - . - . - . - . - . - . - . - . -0.26 - . - . -0.26-0.26-0.26 - . - . - . - . - . - . - . - . - . - . - . - . - . - . - . - . - . -0.478-0.477-0.475-0.471 1/8 1/5 1/3 E -0.477-0.475-0.473-0.471 1/8 1/5 1/3 -0.48-0.479-0.477-0.475-0.473 1/8 1/5 1/3 E Δ T , m Δ T y Stripe:m
FIG. 10. (Color online) ( a, b ) Schematic pattern of J bond-energyat J = 0 . and J = 0 . for bond dimension D = 8 . ( c, d, e ) The U (1) -symmetric iPEPS ground-state energy initialized by N´eel andStripe states at J = { . , . , . } , respectively. ( d ) The orderparameters ∆ T y and m at the point J = 0 . . It shows both N´eeland AFM Stripe states exists at this point. and ξ d ∼ D . , respectively. A very recent SU (2) -iPEPSstudy has suggested that spin correlation length diverges lin-early with environment bond dimension χ , ξ s ∼ χ (although,in contrast, we do not observe any strong dependency on χ inour calculations).The pattern of local nearest neighboring bond energy re-veals that the nature of the VBS order is of the columnar type.The associated VBS order parameter increases monotonicallyup to the point J ≈ . , where a first-order phase transitionoccurs. In comparison with the plaquette VBS order predictedby DMRG simulations, both phases seem to be quite compet-itive. We have estimated transition point at J = 0 . based on hysteresis analysis. At this point both associated or-der parameters of the columnar VBS and the AFM Stripe arenon-zero.Our study clearly shows that the iPEPS ansatz finds a non-0zero N´eel order parameter in the range of . < J < . ,where DMRG studies predict a possible gapless phase. It isan interesting direction to improve both methods further to ob-tain more accurate estimation of the relevant order parametersand reach a rigorous conclusion for that phase. A natural nextstep is to apply the method, determining relevant symmetricsectors, to the SU (2) -symmetric iPEPS ansatz, which mightimprove the accuracy similar to the U (1) case. In addition,for models with long-range interactions, such as J − J − J Heisenberg models, an efficient generalization of the updatescheme is needed. Furthermore, using the U (1) -symmetriciPEPS ansatz for larger-spin systems (defined on different ge-ometries) to characterize different quantum phases is another direction of further studies. ACKNOWLEDGMENTS
We thank Shou-Shu Gong for stimulating discussions. Wealso acknowledge Mac Lee for reading the manuscript. Thisresearch is supported by National Science Foundation GrantsPREM DMR-1205734 and DMR-1408560. We use
Uni10library , an open-source library, to build and perform theiPEPS algorithms introduced in this paper. L. Balents, Nature , 199 (2010). L. Savary and L. Balents, Reports on Progress in Physics ,016502 (2016). X. G. Wen and Q. Niu, Phys. Rev. B , 9377 (1990). N. Read and S. Sachdev, Phys. Rev. Lett. , 1773 (1991). I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. , 799 (1987). N. Read and S. Sachdev, Phys. Rev. Lett. , 1694 (1989). N. Read and S. Sachdev, Phys. Rev. B , 4568 (1990). T. Senthil, A. Vishwanath, L. Balents, S. Sachdev,and M. P. A. Fisher, Science , 1490 (2004),http://science.sciencemag.org/content/303/5663/1490.full.pdf. X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B , 155138(2010). X.-G. Wen,
Quantum field theory of many-body systems: fromthe origin of sound to an origin of light and electrons (OxfordUniversity Press, Oxford, 2007). X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B , 035107(2011). R. Melzi, P. Carretta, A. Lascialfari, M. Mambrini, M. Troyer,P. Millet, and F. Mila, Phys. Rev. Lett. , 1318 (2000). R. Nath, A. A. Tsirlin, H. Rosner, and C. Geibel, Phys. Rev. B , 064422 (2008). T. Koga, N. Kurita, M. Avdeev, S. Danilkin, T. J. Sato, andH. Tanaka, Phys. Rev. B , 054426 (2016). L. B. IOFFE and A. I. LARKIN, InternationalJournal of Modern Physics B M. E. Zhitomirsky and K. Ueda, Phys. Rev. B , 9007 (1996). L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173 (2000). K. Takano, Y. Kito, Y. ¯Ono, and K. Sano, Phys. Rev. Lett. ,197202 (2003). M. Mambrini, A. L¨auchli, D. Poilblanc, and F. Mila, Phys. Rev.B , 144422 (2006). L. Capriotti, F. Becca, A. Parola, and S. Sorella, Phys. Rev. Lett. , 097201 (2001). G.-M. Zhang, H. Hu, and L. Yu, Phys. Rev. Lett. , 067201(2003). J. Sirker, Z. Weihong, O. P. Sushkov, and J. Oitmaa, Phys. Rev. B , 184420 (2006). R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Kr¨uger,and J. Richter, Phys. Rev. B , 214415 (2008). K. S. D. Beach, Phys. Rev. B , 224431 (2009). J. Richter and J. Schulenburg, The European Physical Journal B , 117 (2010). L. Wang, D. Poilblanc, Z.-C. Gu, X.-G. Wen, and F. Verstraete,Phys. Rev. Lett. , 037202 (2013). S. Sachdev and R. N. Bhatt, Phys. Rev. B , 9323 (1990). R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys.Rev. B , 7278 (1999). A. V. Chubukov and T. Jolicoeur, Phys. Rev. B , 12050 (1991). V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. B , 195119(2009). M. Sadrzadeh, R. Haghshenas, S. S. Jahromi, and A. Langari,Phys. Rev. B , 214419 (2016). S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P. A.Fisher, Phys. Rev. Lett. , 027201 (2014). H.-C. Jiang, H. Yao, and L. Balents, Phys. Rev. B , 024424(2012). L. Wang and A. W. Sandvik, arXiv preprint arXiv:1702.08197(2017). W.-J. Hu, F. Becca, A. Parola, and S. Sorella, Phys. Rev. B ,060402 (2013). S. Morita, R. Kaneko, and M. Imada, Journal ofthe Physical Society of Japan , 024720 (2015),http://dx.doi.org/10.7566/JPSJ.84.024720. A. W. Sandvik, Phys. Rev. B , 134407 (2012). L. Wang, Y.-J. Kao, and A. W. Sandvik, Phys. Rev. E , 056703(2011). M. Lubasch, J. I. Cirac, and M.-C. Bauls, New Journal of Physics , 033014 (2014). L. Wang, Z.-C. Gu, F. Verstraete, and X.-G. Wen, Phys. Rev. B , 075143 (2016). D. Poilblanc and M. Mambrini, Phys. Rev. B , 014414 (2017). P. Corboz, R. Or´us, B. Bauer, and G. Vidal, Phys. Rev. B ,165104 (2010). P. Corboz, J. Jordan, and G. Vidal, Phys. Rev. B , 245119(2010). H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Or´us,Phys. Rev. B , 035142 (2015). P. Corboz and F. Mila, Phys. Rev. B , 115144 (2013). B. Bauer, P. Corboz, R. Or´us, and M. Troyer, Phys. Rev. B ,125106 (2011). F. Verstraete, V. Murg, and J. Cirac, Advances in Physics , 143(2008), http://dx.doi.org/10.1080/14789940801912366. R. Ors, Annals of Physics , 117 (2014). S. Singh, R. N. C. Pfeifer, and G. Vidal, Phys. Rev. A , 050301(2010). S. Singh, R. N. C. Pfeifer, and G. Vidal, Phys. Rev. B , 115125(2011). -0.00200.0020.0040.0060.0080.010.0120.0140.016 1/9 1/5 1/3 C s -15-13-11-9-7-5-3 5 15 20 25 30 35 40 l og ( C s )(r) rD=4D=5D=6D=7D=8 FIG. 11. (Color online) The extrapolated data points for the corre-lation function. ( a ) The spin-spin correlation function versus bonddimension /D for different value of r . A linear fit is used to es-timate the extrapolate data in the D → ∞ . ( b ) Log-linear plot ofthe dimer-dimer correlation function versus distance r , where slopesshow inverse of the correlation length ξ − d . H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. ,090603 (2008). Z.-C. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B , 205116(2008). Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang,Phys. Rev. B , 045139 (2012). R. Or´us, Phys. Rev. B , 205117 (2012). R. Or´us and G. Vidal, Phys. Rev. B , 094403 (2009). T. Nishino and K. Okunishi, Journal of the Physical Society ofJapan , 3040 (1997), https://doi.org/10.1143/JPSJ.66.3040. A. W. Sandvik, Phys. Rev. B , 11678 (1997). Y.-J. Kao, Y.-D. Hsieh, and P. Chen, Journal of Physics: Confer-ence Series , 012040 (2015). M. Suzuki, Physics Letters A , 319 (1990). P. Corboz, T. M. Rice, and M. Troyer, Phys. Rev. Lett. ,046402 (2014). Y.-K. Huang, P. Chen, and Y.-J. Kao, Phys. Rev. B , 235102(2012). Our results support a divergent behavior for the dimer correlationlength at the critical point J ∼ . . However, in order to obtaina more rigorous conclusion, to show that it does not tend to a finitevalue in the D → ∞ limit, the larger bond dimension is required,which is beyond our current capability. The typical behavior of correlation function, critical exponentsand correlation length are almost the same in the range J ∼ . − . . So, we mainly rely on vanishing N´eel order parameterto estimate the critical point. J. Lou, A. W. Sandvik, and N. Kawashima, Phys. Rev. B ,180414 (2009). Appendix A: Extrapolated data for the correlation functions
In this section, we provide further data points of the corre-lation functions and discuss the extrapolation procedure usedin the estimation of the critical exponents. In order to esti-mate, e.g., the spin critical exponent, we first obtain the spin-spin correlation function C st ( r ) for the large bond dimensions D ∼ − . Then, we use a linear fit (in /D ) to extrapolate C st ( r ) in the D → ∞ limit; As depicted in Fig. 11-(a), wehave plotted C st ( r ) as a function of /D and have shown thelinear fits for different values of distance r ∼ − . A lin-ear fit seems to provide reliable estimation of the extrapolateddata points. We finally use the these data points to estimatethe exponents as shown in Fig. 9-(a).We also report in Fig. 11-(b), the log-linear plot of thedimer-dimer correlation function versus large distance r (cid:29) .The slopes reveal the inverse of the spin correlation length ξ − d , which seems to weakly depend on D . We use the slopesin Fig. 9-(d) to study D-dependence behavior of ξ dd