# A scaling hypothesis for projected entangled-pair states

Bram Vanhecke, Juraj Hasik, Frank Verstraete, Laurens Vanderstraeten

AA scaling hypothesis for projected entangled-pair states

Bram Vanhecke, ∗ Juraj Hasik, † Frank Verstraete, and Laurens Vanderstraeten Department of Physics and Astronomy, University of Ghent, Krijgslaan 281, 9000 Gent, Belgium Laboratoire de Physique Th´eorique UMR5152, C.N.R.S. and Universit´e de Toulouse,118 rte de Narbonne, 31062 Toulouse, France

We introduce a new paradigm for scaling simulations with projected entangled-pair states (PEPS)for critical strongly-correlated systems, allowing for reliable extrapolations of PEPS data with rela-tively small bond dimensions D . The key ingredient consists of using the eﬀective correlation length ξ for inducing a collapse of data points, f ( D, χ ) = f ( ξ ( D, χ )), for arbitrary values of D and theenvironment bond dimension χ . As such we circumvent the need for extrapolations in χ and can usemany distinct data points for a ﬁxed value of D . Here, we need that the PEPS has been optimizedusing a ﬁxed- χ gradient method, which can be achieved using a novel tensor-network algorithm forﬁnding ﬁxed points of 2-D transfer matrices, or by using the formalism of backwards diﬀerentiation.We test our hypothesis on the critical 3-D dimer model, the 3-D classical Ising model, and the 2-Dquantum Heisenberg model. Although variational algorithms over the class of pro-jected entangled-pair states (PEPS) are by now wellestablished for simulating ground states of 2-D quan-tum spin systems and free energies of 3-D statisticalmechanical systems , the accurate simulation of criti-cal points and gapless phases in (2+1)-D remains elusivedue to the prohibitive computational cost of simulatingPEPS with a large bond dimension. Indeed, PEPS arethe natural higher-dimensional generalization of matrixproduct states (MPS), but, unlike for MPS, expectationvalues cannot be computed exactly as that would requirethe contraction of an inﬁnite two-dimensional tensor net-work. Instead, such networks are typically contracted byconstructing eﬀective environments – either as a bound-ary MPS or with corner-transfer matrices – bothof which carry a bond dimension χ controlling the errorin the environment.The investigation of critical systems by PEPS has beenseverely hampered by the presence of these two controlparameters. In the case of MPS, it has been realizedthat working at ﬁnite bond dimension induces a relevantperturbation in the system, much in the same way as ﬁ-nite sizes enter in exact diagonalization or Monte-Carlotechniques. As an analogue to ﬁnite-size scaling , thetheory of ﬁnite-entanglement scaling identiﬁes an in-duced length scale due to the entanglement compressioncaused by ﬁnite χ . This has allowed MPS simulationsto reach a high precision in the determination of criti-cal points and exponents, both for 1-D quantum systemsand problems in 2-D statistical mechanics.Recently, there have been two crucial developmentsin realizing the same program for PEPS simulations.First, variational algorithms were devised that allow toﬁnd an optimal PEPS approximation for a given modelHamiltonian or transfer matrix at a certain valueof the bond dimension D – in contrast to the simple-update or full-update schemes that were shown to yieldsuboptimal PEPS results. Second, the correlation lengthof the optimal PEPS at a certain bond dimension, D ,was identiﬁed as a suitable length scale for extrapolatingvariational energies and other observables . In both these advances, the role of the environment bond dimen-sion χ was not given any meaningful content: one simplyassumes that all optimizations are done for large enoughvalues of χ , and correlation lengths are then extrapolatedto get rid of any χ dependence. This makes PEPS simula-tions very costly, as the χ -limit is often very tough if notimpossible to reach, especially for large D . Moreover,as an implication, only a very limited number of datapoints can contribute to the extrapolation – i.e., a singlepoint for every value of D accessible by computationalresources. Contrary to MPS, which can be optimized upto bond dimensions of order O (10 ), the gradient-basedPEPS optimizations have so far been limited to D < χ and ﬁnite- D data on anequal footing in PEPS extrapolation schemes. In orderfor this to be possible, the optimized PEPS must be anextreme point of the cost function for a given value of D and χ . That is, we need to solve the optimizationproblem A = arg min A f { D,χ } ( A, ¯ A ) (1)for tensor(s) A forming the PEPS ansatz where the costfunction is the value of the energy (quantum) or free en-ergy (classical) as a result from the approximate PEPScontraction at a ﬁxed value of χ . In the case of two-dimensional PEPO transfer matrices, we show that thiscan be achieved by gradient optimization, since the χ -dependent gradient of the cost function is obtained by asimple contraction (see Appendix). For quantum Hamil-tonians, one can construct a complex time evolutionPEPO and do the same, or the optimization can bedone by supplementing existing contraction algorithmswith a backwards diﬀerentiation step, as e.g. done usingautomatic diﬀerentiation .The crucial contribution of this paper is the followinghypothesis: close to a critical point the dependence ofexpectation values (cid:104) O (cid:105) of an optimized PEPS wavefunc-tion – optimized in the sense of Eq. 1 – as a function of D and χ exhibits a data collapse when expressed throughthe induced correlation length ξ of the environment –the a r X i v : . [ qu a n t - ph ] F e b FIG. 1. Residual entropy per site for the dimer cov-ering problem on the cubic lattice. We have opti-mized PEPS tensors for D = (2 , , , ,

6) and χ =(20 , , , , , , , , , , D > x power law (see inset) forwhich the origin is to us an enigma. one used during optimization– with bond dimension χ : (cid:104) O (cid:105) ( D, χ ) = (cid:104) O (cid:105) ( ξ ( D, χ )) . (2)The consequences of the above hypothesis are twofold.On the one hand, it implies that even optimized PEPSwith relatively small χ can be used for accurate extrapo-lations, provided that they are in the scaling regime. Onthe other hand, it allows for a substantial increase in theamount of data relevant for extrapolation as much moredata points can be constructed by varying χ for each D considered. We want to stress that this hypothesis onlyholds for PEPS algorithms that yield extremal points forthe cost function in Eq. 1, a condition that is not satisﬁedfor simple and full-update algorithms and for earlier vari-ational algorithms , but which is satisﬁed for the exact( D, χ )-gradient algorithms such as the one presented inthe Appendix and the ones based on backwards diﬀeren-tiation.We ﬁrst illustrate this hypothesis by scaling the free en-ergy and order parameters of statistical mechanical prob-lems in three dimensions, where we optimize PEPS ﬁxedpoints of the 2-D transfer matrix with a gradient op-timization at ﬁxed values of D and χ (see Appendix).First, we consider the partition function of dimer cov-erings on a cubic lattice. This problem was ﬁrst tackledwith PEPS in Ref. 7, but here we optimize at ﬁxed valuesof D and χ allowing for better convergence of the PEPSoptimization in addition to enabling us to test our PEPSscaling scheme. The data is condensed on a single plotin Fig. 1. We observe that, indeed, the PEPS data fordiﬀerent ( D, χ ) all lie on the same curve, even for rela-

FIG. 2. Magnetization of the classical 3-D Ising model.We have optimized PEPS tensors for D = (2 , ,

4) and χ = (10 , , , , , ,

50) for 100 equally spaced tempera-tures around the critical temperature. The same MPS bonddimension was used for the contraction of the double- andtriple-layer (see appendix). To plot the data, we shift thetemperature by T c = 4 . ) and rescaled both axes with appropriate pow-ers β = 0 . ν = 0 . ξ (obtained from theconformal bootstrap method ), with ξ the correlation lengthof the double-layer MPS environment. In the inset we providea partial data collapse (also showing the abundance of datapoints), showing the crossing of the rescaled magnetization atthe critical temperature. tively small values of χ . An extrapolation in 1 /ξ (mo-tivated by the inset) gives rise to an extremely precisevalue for the residual entropy, S = − . of S = − . ± . D, χ data for the magnetization, shownin Fig. 2. To perform the data collapse, we use the criticalexponents calculated by bootstrap and ﬁnd the temper-ature T c that optimizes the collapse: T c = 4 . T c = 4 . .Finally we consider the square-lattice quantum Heisen-berg model. Here, the PEPS wavefunction for the groundstate is completely described by a single real tensor A with C v point group symmetry , so we can use theoriginal straightforward and highly-eﬃcient formulationof the corner-transfer matrix algorithm . In order to eval-uate the gradient of the variational energy, the wholeprocess of evaluating the energy expectation value is dif-ferentiated by backwards diﬀerentiation . Again,we optimize PEPS tensors for diﬀerent ( D, χ ), and plotthe energy and magnetization as a function of inversecorrelation length in Fig. 3. The staggered magnetiza-tion can be extrapolated as a linear function of 1 /ξ , as FIG. 3. Data for the 2-D Heisenberg model for D = (6 , , χ from 17 up to 200 for D = 6 , D = 8. Top: energy vs. inverse correlation length of theenvironment MPS. Bottom-left: energy vs. the cube of theinverse correlation length. Bottom-right: magnetization vsthat same inverse correlation length. Fits are performed onall the data and are shown in black. was noted in Ref. 21, and is here motivated by a plot inFig. 3; we ﬁnd the extrapolated value m s = 0 . /ξ , as shown in theleft-bottom plot of Fig. 3, to e = − . values e = − . m s = 0 . U (1) symmetry onthe PEPS (not on the environment) and select the classthat previously proved to be most eﬀective . Notice that the D = 2 data (blue indicators) for thedimer model are clearly not in the scaling regime. Wehave included them as an illustration of the fact that highenough bond dimensions must be chosen to reach thescaling regime. For all the above results we implementsome ﬁltering of the data, imposing that the magneti-zation decreases monotonically in D and χ separately .This allows us to identify data points that seemed com-pletely converged, based on the norm of the gradient,but are in fact stuck at a local minimum and also do notsatisfy the scaling hypothesis. This can be done morethoroughly for the 3-D Ising model, as we can also im-pose monotonicity in temperature in addition to D and χ . Finally, we always use the largest length scale in theboundary MPS. In the case of Heisenberg model it coin-cides with the spin-spin correlation length, but in othermodels it might prove more fruitful to use some smallercorrelation length associated to a more relevant correla-tor; this can be achieved by using symmetric tensor net-works or direct ﬁtting of the relevant correlation function.In conclusion, we have formulated a scaling hypothe-sis for PEPS and identiﬁed an induced correlation length ξ ( D, χ ) as the only relevant parameter when the PEPSis optimized using variational gradient methods at ﬁxedbond dimensions D and χ . This allows for simulatingquantum critical points in 2+1D with much greater accu-racy and smaller cost. We conjecture that a similar scal-ing hypothesis holds in any dimension, therefore openingup the possibility of simulating quantum critical pointsin 3-D and higher. Acknowledgments. — The authors would like to thankFabien Alet and Andreas L¨auchli for inspiring discus-sions. This work was supported by the Research Foun-dation Flanders (G0E1820N) and the ERC grant QUTE(647905). Part of the computations were carried out onthe HPC resources of CALMIP supercomputing centerunder the allocation 2020-P1231. JH is supported bythe TNSTRONG ANR-16-CE30-0025 and the TNTOPANR-18-CE30-0026-01 grants awarded by the French Re-search Council. ∗ [email protected] † [email protected] F. Verstraete and J. I. Cirac, “Renormalization algorithmsfor Quantum-Many Body Systems in two and higher di-mensions,” cond-mat/0407066 (2004). J. I. Cirac, D. Perez-Garcia, N. Schuch, and F. Ver-straete, “Matrix product states and projected entan-gled pair states: Concepts, symmetries, and theorems,”arXiv:2011.12127 (2020). J. Jordan, R. Or´us, G. Vidal, F. Verstraete, and J. I.Cirac, “Classical Simulation of Inﬁnite-Size Quantum Lat-tice Systems in Two Spatial Dimensions,” Phys. Rev. Lett. , 250602 (2008). L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Ver-straete, “Gradient methods for variational optimization of projected entangled-pair states,” Phys. Rev. B , 155123(2016). P. Corboz, “Variational optimization with inﬁnite pro-jected entangled-pair states,” Phys. Rev. B , 035133(2016). T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima,Y. Akutsu, and A. Gendiar, “Two-dimensional tensorproduct variational formulation,” Progress of theoreticalphysics , 409–417 (2001). L. Vanderstraeten, B. Vanhecke, and F. Verstraete,“Residual entropies for three-dimensional frustrated spinsystems with tensor networks,” Phys. Rev. E , 042145(2018). M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber,J. Haegeman, and F. Verstraete, “Faster methods for con- tracting inﬁnite two-dimensional tensor networks,” Phys.Rev. B , 235148 (2018). T. Nishino and K. Okunishi, “Corner Transfer MatrixRenormalization Group Method,” Journal of the PhysicalSociety of Japan , 891 (1996). R. Or´us and G. Vidal, “Simulation of two-dimensionalquantum systems on an inﬁnite lattice revisited: Cornertransfer matrix for tensor contraction,” Phys. Rev. B ,094403 (2009). P. Corboz, S. R. White, G. Vidal, and M. Troyer, “Stripesin the two-dimensional t-J model with inﬁnite projectedentangled-pair states,” Phys. Rev. B , 041108 (2011). E. Br´ezin, “An investigation of ﬁnite size scaling,” in

Finite-Size Scaling , Current Physics–Sources and Com-ments, Vol. 2, edited by J. L. Cardy (Elsevier, 1988) pp.12 – 19. J. L. Cardy, “Conformal invariance and universality inﬁnite-size scaling,” Journal of Physics A: Mathematicaland General , L385–L387 (1984). M. E. Fisher and M. N. Barber, “Scaling theory for ﬁnite-size eﬀects in the critical region,” Phys. Rev. Lett. ,1516–1519 (1972). T. Nishino, K. Okunishi, and M. Kikuchi, “Numericalrenormalization group at criticality,” Physics Letters A , 69–72 (1996). F. Pollmann, S. Mukerjee, A. M. Turner, and J. E.Moore, “Theory of Finite-Entanglement Scaling at One-Dimensional Quantum Critical Points,” Phys. Rev. Lett. , 255701 (2009). L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I. La-torre, “Scaling of entanglement support for matrix productstates,” Phys. Rev. B , 024410 (2008). B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo,“Matrix product states for critical spin chains: Finite-size versus ﬁnite-entanglement scaling,” Phys. Rev. B ,075117 (2012). B. Vanhecke, J. Haegeman, K. Van Acoleyen, L. Vander-straeten, and F. Verstraete, “Scaling hypothesis for matrixproduct states,” Phys. Rev. Lett. , 250604 (2019). B. Vanhecke, L. Vanderstraeten, and F. Verstraete,“Symmetric cluster expansions with tensor networks,”arXiv:1912.10512 (2019). M. Rader and A. M. L¨auchli, “Finite Correlation LengthScaling in Lorentz-Invariant Gapless iPEPS Wave Func-tions,” Phys. Rev. X , 031030 (2018). P. Corboz, P. Czarnik, G. Kapteijns, and L. Tagliacozzo,“Finite Correlation Length Scaling with Inﬁnite ProjectedEntangled-Pair States,” Phys. Rev. X , 031031 (2018). P. Czarnik and P. Corboz, “Finite correlation length scal-ing with inﬁnite projected entangled pair states at ﬁnitetemperature,” Phys. Rev. B , 245107 (2019). H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, “Diﬀer-entiable programming tensor networks,” Phys. Rev. X ,031041 (2019). J. Hasik, D. Poilblanc, and F. Becca, “Investigationof the N´eel phase of the frustrated Heisenberg antifer-romagnet by diﬀerentiable symmetric tensor networks,”arXiv:2009.02313 (2020). A. L. Talapov and H. W. J. Bl¨ote, “The magnetization ofthe 3D Ising model,” Journal of Physics A MathematicalGeneral , 5727–5733 (1996). D. Simmons-Duﬃn, “A semideﬁnite program solver forthe conformal bootstrap,” Journal of High Energy Physics , 174 (2015). I. Beichl and F. Sullivan, “Approximating the permanentvia importance sampling with application to the dimer cov-ering problem,” Journal of Computational Physics ,128 – 147 (1999). J. Hasik and G. B. Mbeng, “peps-torch: A diﬀerentiabletensor network library for two-dimensional lattice models,”(2020). A. W. Sandvik, “Finite-size scaling of the ground-state pa-rameters of the two-dimensional heisenberg model,” Phys.Rev. B , 11678–11690 (1997). A. W. Sandvik and H. G. Evertz, “Loop updates for vari-ational and projector quantum monte carlo simulations inthe valence-bond basis,” Phys. Rev. B , 024407 (2010). Note that such a ﬁltering on the cost function (residualentropy for the dimers, free energy for the 3D Ising model,and energy for the Heisenberg model) is not allowed, as itis not guaranteed to be monotonic in χ or D , and there areindeed data points with a lower value of the cost functionthan the extrapolated value. L. Vanderstraeten, J. Haegeman, and F. Verstraete,“Tangent-space methods for uniform matrix productstates,” SciPost Phys. Lect. Notes , 7 (2019).

APPENDIX: DIFFERENTIATINGVARIATIONAL COST FUNCTIONS WITH PEPS

In the main body of the text, we have explained thatit is crucial to diﬀerentiate the PEPS cost function eval-uated at a ﬁxed value of environment bond dimension χ .In the context of two-dimensional quantum systems, dif-ferentiating the χ -dependent cost function can be donethrough the technique of automatic diﬀerentiation, asshown in Refs. 24 and 25. Here, we show that thesame procedure is possible for PEPS approximations oftransfer-matrix ﬁxed points, which arise in the tensor-network representation of 3-D partition functions orthe imaginary-time evolution operator for 2-D quantumHamiltonians , under certain symmetry conditions onthe involved tensors.We start from a uniform projected entangled-pair op-erator (PEPO) for an inﬁnite system, described by a 6-legtensor T as O ( T ) = TTT TTT TTT , (3)which represents a layer of a three-dimensional tensornetwork. Evaluating this partition function amounts toﬁnding the leading eigenvector | Ψ (cid:105) of this transfer matrix,such that the free energy density is given by f = − lim N →∞ N log (cid:18) (cid:104) Ψ | O ( T ) | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) (cid:19) , (4)where N is the diverging system size – we will show thatwe can give meaning to the density in the thermodynamiclimit directly, i.e. without taking the N → ∞ limit.Assuming this PEPO is Hermitian, we will approximatethe ﬁxed point as a uniform PEPS, parametrized by asingle complex-valued tensor A , | Ψ( A ) (cid:105) = A A AAAAA A A , (5)and solve the optimization problem A = arg min A f ( A, ¯ A ) , (6)with f ( A, ¯ A ) = − N log (cid:18) (cid:104) Ψ( ¯ A ) | T | Ψ( A ) (cid:105)(cid:104) Ψ( ¯ A ) | Ψ( A ) (cid:105) (cid:19) (7)to ﬁnd an optimal tensor A . Here, ¯ A is the complex con-jugate of the tensor A . In order to solve this optimization problem, we have to evaluate numerator and denomina-tor in the above expression for a given PEPS, which boilsdown to contracting an inﬁnite two-layer and three-layertensor network. Here we use boundary MPS methods forthe double-layer transfer matrix T = E E E , (8)with A ¯ A = E , (9)and triple-layer transfer matrix T = E E E , (10)with E = A ¯ AT , (11)Here we assume that the tensor T has the symmetry ¯ T = T XX † . (12)with a given unitary matrix X , and that we have chosenthe PEPS tensor with the symmetry constraint ¯ A = A X (13)with X the same unitary matrix. These symmetry con-straints can be generalised a bit, but the important thingis that the double-layer and triple-layer transfer matrixare both Hermitian operators. In the examples of thedimer and Ising model the X matrix is a simple unitmatrix. Vumps equations

We characterize the boundary MPS as the ﬁxed pointof the vumps algorithm; for more details on uniformMPS and the derivation of the ﬁxed-point equations,see Ref. 33. We ﬁrst list the equations for a normal-ized MPS in the center gauge, described by four tensors { A l , A r , C, A c } , A c = A l C = C A r , (14)with the conditions A l ¯ A l = , A r A ∗ r = (15)and the normalization A c ¯ A c = 1 , C ¯ C = 1 . (16)We characterize the ﬁxed point of a given hermitian MPOby the following ﬁxed-point equations A c O G r G l = λ A c (17)and G r G l C = C . (18)Note that the second equation for C follows from the onefor A c . The left and right ﬁxed points are determined bythe following eigenvalue equations OG l ¯ A l A l = λ G l , (19) O ¯ A r A r G r = λ G r , (20)with the normalization G r G l C ¯ C = 1 . (21)Importantly, through the hermiticity of the transfer ma-trix, we also have the ﬁxed-point equations for ¯ A c and ¯ C , O G r G l ¯ A c = λ ¯ A c , (22)and G r G l ¯ C = ¯ C . (23)Given these ﬁxed-point equations and normalizations,the transfer-matrix eigenvalue per site of the MPO isgiven by λ = ¯ A c A c O G r G l (24) Diﬀerentiating the ﬁxed-point equations

All tensors { A l , A r , A c , C, A ∗ l , A ∗ r , A ∗ c , C ∗ , G l , G r } canbe thought of being functionally dependent on O throughthe above ﬁxed-point equations. If we change thetransfer-matrix tensor as O → O + dO , all these ten-sors have a ﬁrst-order change as well. We can ﬁnd theseﬁrst-order tensors by diﬀerentiating the ﬁxed-point equa-tions. In the following we denote the ﬁrst-order tensorswith small letters.We ﬁrst diﬀerentiate the normalization conditions, ar-riving at the equalities a c ¯ A c + A c ¯ a c = 0 , (25) c ¯ C + C ¯ c = 0 . (26)The normalization of G l and G r can be diﬀerentiated togive G r g l C ¯ C + g r G l C ¯ C + G r G l c ¯ C + G r G l C ¯ c = 0 , (27)which, after using the eigenvalue equation for C and C ∗ and the equation 26, gives G r g l C ¯ C + g r G l C ¯ C = 0 . (28)Next, we ﬁnd the ﬁrst-order change in the transfer-matrixeigenvalue, dλ = dO G r G l ¯ A c A c + O G r g l ¯ A c A c + O g r G l ¯ A c A c + O G r G l ¯ A c a c + O G r G l ¯ a c A c . (29)The four last terms can be simpliﬁed using the eigenvalue equations above and, therefore, cancel, so that we ﬁnd dλ = dO G r G l ¯ A c A c . (30)This implies that the ﬁrst-order change in the eigenvalueof the transfer matrix is given by the simple diagramwhere we diﬀerentiate the O tensor. Optimizing the free energy

We now return to the objective function that we wishto optimize. Using the boundary MPS for double- andtriple-layer transfer matrix, characterized by the vumpsﬁxed-point equations, we approximate the free-energydensity as f ( A, ¯ A ) = − log λ ( A, ¯ A ) , (31)and λ = E A c ¯ A c G l G r / G l E A c G r ¯ A c . (32)The function λ ( A, ¯ A ) is a real function of complex pa-rameters. The gradient of the objective function is givenby g = − λ ( A, ¯ A ) G l E A c G r ¯ A c − e A c ¯ A c G l G r − λ ( A, ¯ A ) G l e A c G r ¯ A c (33)with e = ∂ ¯ A E , e = ∂ ¯ A E . (34)Here, the diﬀerentials with respect to the environmenttensors all vanish as we have showed above. This meansthat the gradient g is the exact gradient for the objectivefunction evaluated at a given value of the boundary MPS χχ