High-Precision Thermodynamic and Critical Properties from Tensor Renormalization-Group Flows
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p High-Precision Thermodynamic and Critical Propertiesfrom Tensor Renormalization-Group Flows
Michael Hinczewski and A. Nihat Berker − Feza G¨ursey Research Institute, T ¨UBITAK - Bosphorus University, C¸ engelk¨oy 34684, Istanbul, Turkey, Department of Physics, Ko¸c University, Sarıyer 34450, Istanbul, Turkey, and Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, U.S.A.
The recently developed tensor renormalization-group (TRG) method provides a highly precisetechnique for deriving thermodynamic and critical properties of lattice Hamiltonians. The TRG isa local coarse-graining transformation, with the elements of the tensor at each lattice site playingthe part of the interactions that undergo the renormalization-group flows. These tensor flows aredirectly related to the phase diagram structure of the infinite system, with each phase flowing to adistinct surface of fixed points. Fixed-point analysis and summation along the flows give the criticalexponents, as well as thermodynamic functions along the entire temperature range. Thus, for theferromagnetic triangular lattice Ising model, the free energy is calculated to better than 10 − alongthe entire temperature range. Unlike previous position-space renormalization-group methods, thetruncation (of the tensor index range D ) in this general method converges under straightforwardand systematic improvements. Our best results are easily obtained with D = 24, corresponding to4624-dimensional renormalization-group flows. PACS numbers: 64.60.Ak, 05.10.Cc, 05.70.Jk, 75.10.Hk
I. INTRODUCTION
The tensor renormalization-group (TRG) method, re-cently introduced by Levin and Nave [1] is versatile, ac-curate, and conceptually interesting. The versatility ofTRG comes from being a genuinely local renormalization-group transformation, that is a mapping between lo-cal Hamiltonians on the original and coarse-grained lat-tices, expressed in terms of tensors at sites of each lat-tice. Thus, although TRG has been demonstrated byobtaining very accurate phase transition temperaturesand magnetization curves for frustrated and unfrustratedclassical Ising models on the triangular lattice, we be-lieve that it can be extended to yield the critical prop-erties and entire thermodynamics of complicated multi-critical systems, including quenched random systems, byrenormalization-group flow analysis.In the present work, we show that the renormalization-group flows of the tensor elements can be analyzed inthe same manner as the flows of Hamiltonian interactionparameters, which was not done in the previous work.Each thermodynamic phase region corresponds to a basinof attraction in the space of tensor elements, and criti-cal temperatures and exponents (not obtained before)are deduced from characteristics of the boundaries be-tween the different basins. Practically, this means thatthe phase diagram, thermodynamic, and critical proper-ties of the infinite system can be derived directly fromthe flows. Moreover, the flows themselves have an in-teresting and unconventional nature, since a phase re-gion does not flow to single isolated sink, but rather toa continuous surface of fixed points. Using the ferro-magnetic triangular lattice Ising model as an example,we show that the accuracy of the calculated free energy,critical temperature, and thermal critical exponent canbe easily and systematically improved on by increasing D , the cutoff on the index range of the tensors. In thisability to converge toward exact values with larger cut-offs, TRG is more general and dependable than position-space renormalization-group techniques [2, 3, 4] whosesuccesses have been based on system-specific heuristics[5, 6, 7]. The TRG approach combines the straightfor-ward interpretative framework of traditional renormal-ization group—analysis of flows in a parameter space—with the accuracy of techniques depending on finite-sizescaling of large systems. II. TENSOR RENORMALIZATION-GROUPTRANSFORMATION
The TRG method [1] can be applied to any classicallattice Hamiltonian satisfying the following conditions:(1) it can be expressed in terms of degrees of freedom onthe bonds of the lattice; (2) the Boltzmann weight of aconfiguration can be written as the product of individualBoltzmann weights for each lattice site, that only dependon the bond variables adjoining the site. This is a broadcategory including many common statistical physics sys-tems like the Ising and Potts models, both of which canbe mapped onto this form through duality transforma-tions, as well as all vertex models. For such Hamiltoni-ans, the partition function can be written as a tensor net-work [8, 9]: We start with a lattice of N sites where eachsite has coordination number q and each bond is in one of d possible states. The Boltzmann weight of an individualsite depends on the configuration of the q bonds meetingat the site and can be written as a tensor T i i ··· i q , whereeach index i α runs from 1 to d . For configurations of bondvariables that are not allowed, the corresponding elementof T is zero. The tensor is real-valued and cyclicallysymmetric. The transformation described below works (a) T ijk = ijk S ijν = ijν T ′ νγδ = νγδ (b) kij ml = νmi lj (c) m ihγ νδ = νγδ (d) FIG. 1: (a) Graphical representation of the tensors. (b)Rewiring. (c) Decimation. (d) Renormalization group trans-formation applied to the entire hexagonal lattice, with thefirst arrow showing the rewiring step, and the second arrowshowing the decimation. also for the more general case of complex-valued, cycli-cally symmetric tensors, onto which our original tensornetwork system is mapped after a single renormalizationstep. The partition function is a product over the N sitetensors, with each index contracted between two differenttensors (since each bond is shared between two sites), Z = d X i ,...,i M =1 T i i ··· i q T i i r ··· i s T i i t ··· i u · · · . (1)The position-space renormalization-group transforma-tion of this system allows us to express the partitionfunction equivalently as a tensor network over a coarse-grained lattice with N ′ = N/b vertices, in two spatialdimensions, with a length rescaling factor b >
1. It isaccomplished in two steps, which we call rewiring anddecimation. For simplicity, we shall focus in this workon applying the transformation to the hexagonal lattice,though this two-step procedure can easily be adaptedto a variety of two-dimensional lattices, including thesquare [1] and Kagome lattices.
Rewiring:
Graphically, let us represent the tensor T ijk at each point of the hexagonal lattice as a three-leggedvertex as shown on the left of Fig. 1(a), and a contractionof an index between two tensors as a connection betweentwo vertex legs. We use the convention that the orderof indices in the tensor matches the counterclockwise or-dering of the labels on the vertex legs. The rewiring stepfor a pair of neighboring tensors consists of reconnectingthe bonds in the manner of Fig. 1(b), rewriting them asa contraction of two new tensors S , d X k =1 T ijk T klm = d X ν =1 S miν S jlν . (2) The first two indices in the tensor S miν run up to d , butthe third index ν runs up to d . We shall distinguishthese types of indices by a Greek letter label, and graph-ically depict them as thick legs [Fig. 1(a) center]. To seethat such a rewiring is possible, let us introduce compos-ite indices α ≡ ( m, i ) and β ≡ ( j, l ), and write the tensorcontraction on the left-hand side of Eq. (2) as a d × d matrix M αβ ≡ P k T ijk T klm . Then Eq. (2) becomes M αβ = X ν S αν S βν , (3)or M = SS T . We can find S using the fact that M = M T , as can be checked using the cyclical sym-metry of the T tensors. Any symmetric matrix M ad-mits a variant of singular value decomposition known asTakagi factorization [10]: M = U Σ U T , where U is a d × d unitary square matrix and Σ is a d × d diag-onal matrix containing the singular values of M . If Σ νν is the ν th singular value (assumed ordered from largestto smallest with increasing ν ), then the elements of S are given by S αν = √ Σ νν U αν . The resulting decomposi-tion M = SS T is not unique, because we can replace S by SO , where O is a complex orthogonal matrix ( OO T = I ).However, we shall always use the S given directly by thethe Takagi factorization, so if the ν th singular value Σ νν is nondegenerate, the ν th column of S is uniquely deter-mined up to a factor of ±
1. The rewiring procedure de-scribed above is applied globally to the entire hexagonallattice as shown by the first arrow in Fig. 1(d), groupingthe T tensors into pairs, and replacing each pair by S tensors. Decimation:
The second step in the renormalizationprocedure consists of tracing over the degrees of free-dom in the triangular clusters that are formed after therewiring. As illustrated in Fig. 1(c), each such clustercan be replaced by a single renormalized tensor T ′ , d X m,i,h =1 S hmγ S miν S ihδ = T ′ γδν . (4)Introducing the notation S ( γ ) for the d × d matrix withelements S hmγ at fixed γ , then the above equation hasthe form Tr[ S ( γ ) S ( ν ) S ( δ ) ] = T ′ γδν . (5)It is clear from Eq. (5) that the renormalized tensor willbe cyclically symmetric. Since the S tensor is constructedfrom the unitary matrix U , the elements of T ′ will ingeneral be complex. The second arrow of Fig. 1(d) showsthe decimation applied to every triangular cluster, andthe result is a hexagonal lattice of T ′ tensors with thelattice spacing larger by a factor of b = √ Z . However, therenormalized tensors have a more complicated structurethan the original ones, since the indices of T ′ γδν run up to d . If this procedure is repeated, the index range growsexponentially with iteration number, making numericalimplementation impractical. This is analogous to the dif-ficulty encountered when applying naive position-spacerenormalization to Hamiltonians on lattices, where thenumber of couplings in the Hamiltonian grows with eachiteration. Some kind of approximate truncation is re-quired to keep the complexity of the renormalized modelbounded. In our case this truncation can be done ina straightforward and systematic fashion by setting anupper bound D for the index range. Rather than usethe full d × d matrix S αν to calculate T ′ , we use onlythe first ¯ d columns, where ¯ d ≡ min( d , D ). Using thistruncated d × ¯ d matrix, which we call ¯ S , means thatthe rewiring step is implemented only approximately, M ≈ ¯ S ¯ S T . However, this is the best approximation pos-sible, since the first ¯ d columns of S correspond to the ¯ d largest singular values of M . The decimation step of Eq.(4) is still carried out exactly, but with S replaced by ¯ S .The resulting tensor T ′ γδν has indices that run up to ¯ d . III. RENORMALIZATION-GROUP FLOWSAND THERMODYNAMIC BEHAVIOR OF THETRIANGULAR LATTICE ISING MODEL
In order to illustrate the nature of the renormalization-group flows resulting from the transformation outlinedabove and the methods by which thermodynamic infor-mation can be extracted from them, we turn now to a spe-cific example: the triangular lattice ferromagnetic Isingmodel. This model is mapped onto a tensor networkthrough a duality transformation, giving an equivalentpartition function in terms of bond variables on a hexag-onal lattice.[1] Each bond in this dual lattice has twostates, σ = 1 or σ = −
1, corresponding to a bond be-tween parallel or antiparallel Ising spins. These statesrespectively contribute energies − J < J to the to-tal Hamiltonian H . The tensor for each lattice site hasthe form T ijk = e βJ ( σ i + σ j + σ k )
12 ( σ i σ j σ k + 1) , (6)where β = 1 /k B T , and the factor multiplying the ex-ponential is a projection operator that is equal to 1 forallowed configurations of the bonds and 0 otherwise. Forsimplicity we set J/k B = 1, effectively measuring tem-peratures in units of J/k B . If the triangular lattice has N sites, the dual hexagonal lattice has 2 N sites, and thepartition function is a contraction of the 2 N site tensors.Iterating the renormalization-group transformation,we consider renormalization-group flows in the space oftensor element amplitudes | T αβγ | . For a cyclically sym-metric tensor T ijk with index range d , the maximumnumber of distinct elements is d (2 + d ) /
3. After thefirst few iterations, the index range of the renormalized tensors reaches the cutoff D and remains there, so thatthe flows in the subsequent steps are in a space with di-mension D (2 + D ) /
3. In this work we look at cutoffs inthe range D = 4 to 24, so the most complex flow spacethat we investigate is 4624-dimensional. A large numberof the D (2 + D ) / T = | T | ˜ T , yielding a reduced tensor ˜ T .The T element is a convenient choice, since it remainsnonzero throughout the renormalization-group flow. Therenormalization-group transformation is then applied tothe reduced tensor ˜ T , giving a renormalized tensor T ′ .Let us denote T as T (0) and T ′ as T (1) . The factorizationand renormalization-group transformation are iterated,so at the n th step we have a tensor T ( n ) = | T ( n )111 | ˜ T ( n ) . Aswe shall see below, with this added factorization step theelements of T ( n ) tend to finite limiting values as n → ∞ .Additionally, keeping track of the factors | T ( n )111 | duringthe flow, as in standard position-space renormalization-group theory, allows us to easily determine the free en-ergy of the system in the thermodynamic limit.[12] Defin-ing G ( n ) ≡ ln | T ( n )111 | , the partition function is a con-traction of 2 N tensors T (0) , which we write schemat-ically as Z = ( T (0) ) N = e NG (0) ( ˜ T (0) ) N . After asingle renormalization-group step, this becomes Z = e NG (0) ( T (1) ) N/ = e N ( G (0) + G (1) / ( ˜ T (1) ) N/ . Thus af-ter n steps we find Z = e N P ni =0 − i G ( i ) (cid:16) ˜ T ( n ) (cid:17) N/ n . (7)The free energy per original site βf = − N − ln Z is then βf = − n X i =0 − i G ( i ) − N ln (cid:16) ˜ T ( n ) (cid:17) N/ n . (8)As n → ∞ the elements of ˜ T ( n ) go to finite fixed val-ues, and thus for large N the second term on the rightin Eq. (8) becomes negligible. So the final result forthe free energy is βf = − P ∞ i =0 − i G ( i ) . In practicethe transformation is iterated for n steps until the value T − − − − − − − | f ( D ) ( T ) − f e x a c t ( T ) | / f e x a c t ( T ) D T c FIG. 2: The relative error of f ( D ) ( T ), the free energy per siteof the triangular lattice Ising model calculated using TRGat various cutoffs D = 4 to 24, compared to the exact freeenergy f exact ( T ), as a function of temperature T . The dottedline shows the exact critical temperature T c = 4 / ln 3. of G ( n ) has converged to the limit G ( ∞ ) within the de-sired numerical precision. Then βf ≈ − P n − i =0 − i G ( i ) − G ( n ) P ∞ i = n − i = − P n − i =0 − i G ( i ) − − n G ( n ) .Even for small cutoffs D , the results of this free energycalculation can be remarkably accurate, as seen in Fig.2, which shows the relative error of f ( D ) ( T ), the freeenergy per site using cutoff D , compared to the exactvalue f exact ( T ). The latter is calculated by numericalevaluation of the integral solution in Ref. [13]. Alreadyat D = 4 the TRG free energy is within 0.09% of theexact value at all temperatures, and is considerably moreaccurate than this away from the critical region near T c =4 / ln 3. As noted in Ref. [1], the TRG method is expectedto behave worst at criticality, and this is indeed whatwe see for all D , with the relative error curves peakednear T c . Yet even here there is significant improvementas we go to larger cutoffs. At D = 24, the largest cutoffexamined, the maximum error is 0.0007%. Overall, goingfrom D = 4 to D = 24 we get an improvement betweentwo and three orders of magnitude in the precision of thefree energy result.Further thermodynamic information can be gleaned bylooking in detail at the behavior of the renormalization-group flows. The elements | T ( n ) ijk | go to finite fixedvalues | T ∗ ijk | as n → ∞ , but unlike typical position-space renormalization-group transformations, there areno unique fixed points which act as a basins of attractionfor the low- and high-temperature phases. Instead, ateach different temperature T the system flows to a dif-ferent fixed point | T ∗ ijk ( T ) | . However there is a way todistinguish the low- and high-temperature phases, sinceeach flows to a different continuous surface of fixed pointsin the space of tensor elements. There is a boundary tem-perature T ( D ) c , depending on the cutoff D , such that for T > T ( D ) c the flows tend to one fixed surface, while for D − − − − | T ( D ) c − T c | / T c FIG. 3: The relative error of the critical temperature T ( D ) c for the triangular lattice Ising model calculated from TRGflows at various cutoffs D = 4 to 24, compared to the exact T c = 4 / ln 3. T < T ( D ) c they tend to the other. As D becomes larger, T ( D ) c gives a rapidly converging estimate of the exact crit-ical temperature T c = 4 / ln 3. Fig. 3 shows the relativeerror of T ( D ) c compared to T c for various D . While thedecrease in error is not monotonic with D , there is anoverall trend which appears to be roughly exponential inthe cutoff, so that the error at D = 24, 0 . D = 4.This segregation of the flows between two distinct fixedsurfaces can be seen directly in Fig. 4, which plots onetensor element | T ( n )111 | as a function of iteration number n for D = 8 and D = 12. The flows are shown for differentvalues of | ∆ T | = | T − T ( D ) c | , from 10 − to 10 − . Thereare two curves for each | ∆ T | , one corresponding to T >T ( D ) c , and the other to T < T ( D ) c . The two curves stayclose to one another for a number of iterations, but thenveer off in opposite directions and reach different fixedvalues | T ( ∗ )111 | . These fixed values change continuously as | ∆ T | is varied, as they map out a slice of the two fixedsurfaces corresponding to the low- and high-temperaturephases.As | ∆ T | gets smaller, there is an interesting differencein the flow behaviors of the D = 8 and D = 12 cases.For D = 8, as seen in Fig. 4(a) and in more detail inFig. 4(b), the ∆ T >
T < | T αβγ | elements,since the flows are attracted to a unique critical fixedpoint, which we denote | T c ∗ αβγ | . The smaller the valueof | ∆ T | , the larger the number of iterations which arespent in the vicinity of | T c ∗ αβγ | before flowing to one ofthe fixed surfaces. Once we are in the vicinity of the crit-ical fixed point, we can isolate it to high precision usinga Newton-Raphson procedure. The analysis of the fixedpoint proceeds just as in a standard renormalization-group approach: we calculate a recursion matrix, whoseeigenvalues can be related to the critical exponents of thesystem. To do this, let us denote the nonzero elements of Iteration number n . . . . | T ( n ) | (c) D = 12 . . . | T ( n ) | (b) D = 8 . . . . | T ( n ) | (a) D = 8 10 − − − − − − − − − FIG. 4: The behavior of tensor element | T ( n )111 | as a functionof iteration number n . The different data sets correspond toflows at different temperatures T near T ( D ) c , with each symbolin the legend corresponding to a value of | ∆ T | = | T − T ( D ) c | between 10 − and 10 − . For each | ∆ T | there are two datasets, one for ∆ T >
T < | T ( n )111 | flows with cutoff D = 8.(b) Same as in (a), but zoomed in to see more clearly theflows for ∆ T <
0. (c) | T ( n )111 | flows with cutoff D = 12. | T c ∗ αβγ | , not related by cyclical symmetry, as K through K m . In the case of D = 8, m = 80. For small pertur-bations away from the critical fixed point, the numberand locations of these nonzero elements stay the sameafter a renormalization-group transformation, which al-lows us to numerically evaluate the m × m recursion ma-trix R ij ≡ ∂K ′ i /∂K j . Writing the eigenvalues of R inthe form b y i , i = 1 , . . . , m , we find only one eigenvaluewhere y i >
0, as expected at a critical fixed point. Thisrelevant eigenvalue, which we denote y T , is related to thespecific heat critical exponent α by α = (2 y T − /y T .For D = 8, y T = 1 . α = 0 . − − − − T − T ( D ) c Sp ec i fi c h e a t C D FIG. 5: Specific heat per site C of the triangular lattice Isingmodel, calculated using TRG for various cutoffs D = 4 to24 as a function of T − T ( D ) c in the high-temperature phasenear criticality. The gray curves superimposed on the datapoints are best-fit curves of the form A + B ( T − T ( D ) c ) − α ,with parameters A , B , and α . The thick black line is theexact specific heat as a function of T − T c . D T ( D ) c T ( D ) c − T c y T . × − − . × − . × − . ± . − . × − . ± . . × − . ± . . × − . ± . D , compared to the exact values inthe last row. The thermal eigenvalues y T for D = 4 and 8 arecalculated from the recursion matrix evaluated at the criticalfixed point. For D > y T is from the best-fitresult to the specific heat near T ( D ) c , as plotted in Fig. 5. compare well to the exact values of y T = 1 and α = 0.We can check this result using an alternative approach,by calculating the specific heat from derivatives of thecalculated free energy for small T − T ( D ) c , and we findthat the singularity in the specific heat agrees with the α derived from the thermal eigenvalue y T of the recursionmatrix.For the six values of the cutoff D where we investigatedthe near-critical flows in detail, D = 4 , , , , , D = 4and 8, namely in 24- and 196-dimensional flow spaces.Despite the different dimensionalities of the flow spaces,both cases yielded the same eigenvalue y T , within theprecision of 5 decimal places. Higher values of D showedvery different flow behaviors, as exemplified in Fig. 4(c)for D = 12. Here the ∆ T >
T < | T ( n )111 | do not stay nearly horizontal before divergingto the fixed surfaces: regardless of how small | ∆ T | ismade, the flows do not gravitate toward a unique criticalfixed point, but map out a continuous spectrum of pointswhich attract the flows before the two curves spread out.Using an arbitrary precision version of the TRG algo-rithm implemented in Mathematica , we checked | ∆ T | values as small as 10 − without finding convergence to-ward a critical fixed point. Nevertheless, in these casesfor D > C per site near T ( D ) c as obtained from the numerical derivative of the calcu-lated free energy. We plot C within the high-temperaturephase for various D in Fig. 5, for T − T ( D ) c between 10 − and 10 − . For each D the data points are fit to a function A + B ( T − T ( D ) c ) − α , which provides an accurate descrip-tion of the singularity, and the best-fit value of α is usedto determine the thermal eigenvalue y T . The results arelisted in Table I. For D = 12, y T = 1 . ± . D = 4or 8, but this value improves as D is increased, reaching y T = 1 . ± . D = 24. Moreover, as seen inFig. 5, our calculated C curve for D = 24 nearly overlapsthe exact C as a function of T − T c , which diverges witha logarithmic singularity at the critical point. IV. CONCLUSION
In summary, we have seen that the flows of the tensorelements in the TRG transformation can be used to ex-tract the phase diagram structure and critical behaviorof a classical two-dimensional lattice Hamiltonian. For the triangular lattice Ising model, the low- and high-temperature phase regions are basins of attraction fortwo distinct surfaces of fixed points. The boundary be-tween these basins defines a critical temperature T ( D ) c ,dependent on the TRG cutoff D . At small cutoffs suchas D = 4 and 8, the flows near the boundary betweenthe basins are controlled by a critical fixed point, whileat higher D the flows show more complicated behavior,never converging at a unique point. In the former casethe thermal exponent y T is found from the eigenvalues ofthe recursion matrix at the critical fixed point, while inthe latter we can deduce the exponent from the scaling ofthe calculated specific heat near T ( D ) c . The free energyat all temperatures systematically converges to the ex-act Ising result with increasing D , particularly fast awayfrom the critical region. For the critical properties theimprovement is not monotonic in D , but both the criti-cal temperature T ( D ) c and the exponent y T tend towardthe exact Ising values at larger cutoffs. With very modestcomputational effort, the TRG method provides an accu-rate portrait of global phase diagram characteristics. Itthus warrants further study, both on applications to othertwo-dimensional lattice models, and possible generaliza-tion to systems with quenched randomness [14] and/orhigher spatial dimensions. Acknowledgments
This research was supported by the Scientific and Tech-nical Research Council (T ¨UB˙ITAK) and by the Academyof Sciences of Turkey. [1] M. Levin and C.P. Nave, arXiv:cond-mat/0611687.[2] Th. Niemeijer and J.M.J. van Leeuwen, Phys. Rev. Lett. , 1411 (1973).[3] A.A. Migdal, Zh. Eksp. Teor. Fiz. , 1457 (1975) [Sov.Phys. JETP , 743 (1976)].[4] L.P. Kadanoff, Ann. Phys. (N.Y.) , 359 (1976).[5] L.P. Kadanoff, Phys. Rev. Lett. , 1005 (1975).[6] L.P. Kadanoff, A. Houghton, and M.C. Yalabık, J. Stat.Phys. , 171 (1976).[7] A.N. Berker and M. Wortis, Phys. Rev. B , 4946(1976).[8] I. Markov and Y. Shi, arXiv:quant-ph/0511069. [9] Y. Shi, L. Duan, and G. Vidal, Phys. Rev. A , 022320(2006).[10] T. Takagi, Japan. J. Math. , 83 (1925).[11] E. Anderson et. al. , LAPACK Users’ Guide (SIAM,Philadelphia, PA, 1999), 3rd ed.[12] M. Nauenberg, J. Math. Phys. , 703 (1975).[13] G.H. Wannier, Phys. Rev. , 357 (1950); erratum : Phys.Rev. B , 5017 (1973).[14] A. Falicov, A.N. Berker, and S.R. McKay, Phys. Rev. B51