Accurate exponents from approximate tensor renormalizations
We explain the recent numerical successes obtained by Tao Xiang’s group, who developed and applied Tensor Renormalization Group methods for the Ising model on square and cubic lattices, by the fact that their new truncation method sharply singles out a surprisingly small subspace of dimension two. We show that in the two-state approximation, their transformation can be handled analytically yielding a value 0.964 for the critical exponent much closer to the exact value 1 than 1.338 obtained in the Migdal-Kadanoff approximation. We propose two alternative blocking procedures that preserve the isotropy and improve the accuracy to and 0.993 respectively. We discuss applications to other classical lattice models, including models with fermions, and suggest that it could become a competitor for Monte Carlo methods suitable to calculate accurately critical exponents, take continuum limits and study near-conformal systems in arbitrarily large volumes.
The Renormalization Group (RG) ideas have triggered considerable conceptual and numerical progress in many branches of physics cardy96 (); Meurice:2011wy (). However, the basic method to thin down the number of degrees of freedom in configuration space wilson74 (), often called “block spinning”, has remained a formidable computational challenge for most classical lattice models (e. g., spin models and lattice gauge theories). A few years ago, inspired by the so-called tensor network states vc (); cv () introduced in the context of the density matrix RG method PhysRevLett.69.2863 (); uli (), a Tensor RG (TRG) approach of two-dimensional (2D) classical lattice models was proposed PhysRevLett.99.120601 (). Successful approximations PhysRevLett.99.120601 (); PhysRevLett.103.160601 (); PhysRevB.81.174411 () were found for the Ising model on honeycomb and triangular lattices.
Very recently, the TRG method was successfully extended to the Ising model on square and cubic lattices by Tao Xiang’s group PhysRevB.86.045139 (). There are two important ingredients in their calculations. First, their formulation allows an exact block spinning procedure which separates neatly the degrees of freedom inside the block, which are integrated over, from those kept to communicate with the neighboring blocks. As explained below, this provides a more systematic way to implement ideas initiated by Migdal Migdal:1975zf () and Kadanoff Kadanoff:1976jb () (abbreviated as MK hereafter). The indices of the tensors run over some finite vector space of “states” associated with finite volume link configurations. Second, they used a new method, based on higher order singular value decomposition, which selects in a very economical way the most important states that insure the communication among the blocks. Calculations using of the order of 20 states can be carried on a laptop computer. The excellent agreement found with the Onsager solution in 2D for arbitrarily large volume suggests that TRG-based methods could become competitors for conventional Monte Carlo methods.
In this Letter, we show that the truncation method of Ref. PhysRevB.86.045139 () for the 2D Ising model sharply singles out a two-dimensional subspace of states. Keeping only these two states, we show that we can construct approximate RG transformations with 3 or 4 parameters, find the nontrivial fixed point and obtain precise estimates of the critical exponent associated with the correlation length from a linear analysis. The accuracy of the estimates is significantly better than for textbook examples such as the MK approximation Migdal:1975zf (); Kadanoff:1976jb (); Martinelli:1980sr (), the so-called approximate recursion formula wilson71b () or other hierarchical approximations baker72 (); hmreview ().
The TRG formulation for the Ising model can be extended to nonlinear sigma models and recent numerical implementations for xiangprogress () indicate an optimistic outlook. It seems also possible to formulate models with local invariance and avoid sign problems. In this context, it is important to understand why the method works so well for the Ising model.
The paper is organized as follows. We review the basic TRG formulation for the Ising model on a square lattice emphasizing the connection with the MK ideas. We then consider the cases of an isotropic blocking (as in the Migdal recursion) and an anisotropic blocking (as in the Kadanoff version and in Ref. PhysRevB.86.045139 ()). We also propose a new type of accurate isotropic projection based on a transfer matrix. We briefly discuss the 3D Ising model and how the method can be applied for lattice fermions. The implications for the study of near-conformal systems and the calculations of critical exponents are discussed in the conclusions.
We consider the nearest neighbor Ising model on a square lattice with an inverse temperature . For easy reference, we stay close to the notations of Ref. PhysRevB.86.045139 () where it is shown that the partition function can be written as
The tensor is attached to each site and Tr is a short notation for contractions over the links joining nearest neighbors on the lattice. The horizontal indices and vertical indices take the values 0 and 1. The tensor is zero for an odd number of 1. For an even number of 1, a factor (with appears for each 1 irrespectively of the direction. The partition function can be interpreted as a sum over intermediate states attached to the links. The reader familiar with the high temperature expansion of the model will recognize that this expression reproduces exactly the proper closed paths with the proper weights.
We now use this reformulation to blockspin. First, we follow Migdal Migdal:1975zf () by using an isotropic procedure. We consider a square block enclosing 4 sites and sum over the states, inside the block, associated with the nearest neighbor links joining these 4 points. This defines a new rank 4 tensor where each index now takes 4 values.
where is a notation for the product states. For reasons that will become clear later, we use the convention: . Later, we also use the ket notation for etc.. This is represented graphically in Fig. 1.
The new tensor can be used to define an exact expression of the partition function of the same form as (1), however the number of states proliferates as after steps and approximations are needed in order to get an expression useful for practical purposes.
which can be put in a canonical form by using a Higher Order Singular Value Decomposition defined by a unitary transformation on each of the four indices (see Ref. PhysRevB.86.045139 () for justifications and refinements). The unitary transformation for each index is the one that diagonalizes the symmetric tensor obtained by summing over all the other indices in the following way:
We now consider approximations where for each index we only keep the two states which correspond to the two largest eigenvalues of . With our convention on the product states, this matrix is block diagonal because it does not connect states with even number of 1’s () to states with odd number of 1’s (). Numerically, we found that the two smallest eigenvalues are always very small compared to the largest one for any value of and that the second largest one is small at small and almost as large as the largest one at large . This situation is illustrated for the initial step in Fig. 2. After iterations the gap sharpens as if going to smaller for or larger for .
The two new states have the form
The angle is obtained by diagonalizing the even-even block. The symmetric form of is a consequence of , itself due to reflection symmetry.
With this rather drastic projection, we obtain a new rank 4 tensor with indices taking again two values and the same parity selection rule. We divide the new tensor by a constant in such a way that, dropping all the primes herafter, . In the isotropic case, the reflection symmetry imposes that , , . For the initial tensor, we have what we call later the “Ising condition”:
but this property is not preserved by the blocking procedure which can be expressed as a mapping of the three dimensional parameter space into itself. This can be expressed as rational expressions involving the values of the tensor and trigonometric functions of . They are analytical expressions and derivatives can be taken to calculate the linearized map. They are not written explicitly here but can be obtained easily with symbolic manipulation programs.
A fixed point is found for not far from the exact value 0.44068. The fixed point is approximately and shows significant departure from the Ising condition (6). The eigenvalues of the linearized RG transformation are =2.0177, =0.2022 and =0.09967. The scaling factor is 2, and this implies the values for the critical exponents and which can be compared with the exact values 1 and 2 respectively.
For comparison, the also isotropic Migdal recursion for the Ising model and can be written as
The fixed point is at and corresponds to and = 1.338. This approximation can be seen as a two-state approximation which in addition requires the Ising condition. Improvements of MK discussed in Ref. Martinelli:1980sr () lead to a value of = 0.796 at first order in their expansion parameter and 0.93(1) at second order, the error bar coming from the use of different Padé approximants.
This suggests that the Ising condition is too restrictive. It is nevertheless possible to modify the angle in Eq. (5) in such a way that Ising condition remains valid. This is nontrivial because the Ising condition amounts to two equations. However, an explicit calculation shows that if the two conditions are satisfied and the mapping takes the form
which unfortunately has only the low temperature fixed point.
The above results can be compared with the two-state truncation for the anisotropic blocking used in Ref. PhysRevB.86.045139 () where a first projection occurs after blocking pairs of vertical sites as already introduced in Eq. (3). The complete transformation is then obtained by repeating the procedure with a horizontal blocking as in Kadanoff’s proposal Kadanoff:1976jb (). We use our previous notations but with and subscripts denoting the horizontal and vertical couplings respectively: and . We keep the same notations for and which are invariant under the interchange of the vertical and horizontal directions. We have now a map with four parameters. In the Ising condition, we need to replace by . We take initial values that satisfy the Ising condition and are isotropic (). The critical value is then with a nontrivial fixed point for approximately =0.21358, =0.37924, =0.41998 and = 0.27177, which is clearly anisotropic and violates the Ising condition. We have an additional, intermediate, eigenvalue which is approximately 0.657 and specific to the anisotropic case. This value is not very far from unity which is why it requires more fine-tuning of to get rid of the irrelevant directions. The first and third eigenvalues are 2.052 and 0.1934 respectively which implies = 0.964 and = 2.37.
This anisotropic version can be compared with the Kadanoff recursion Kadanoff:1976jb () for =2, where first, the horizontal bonds are slided vertically with doubled and squared. This corresponds to doubling the vertical lattice spacing first as we did above. After repeating the procedure with horizontal moves, we obtain the Migdal recursion of Eq. (7) for while for we obtain
which corresponds to reversing the order of the two operations. The fixed point is and =0.295598. The eigenvalue is the same in both directions and the value of identical to the Migdal case.
We have also considered the isotropic map with a different truncation. Instead of Eq. (4), we use
The trace of this matrix is the partition function for a model with periodic boundary conditions. This gives a slightly displaced fixed point at . =0.42229, =0.28637 and =0.27466. The values of the exponents are and . It should be noted that in this case the two small eigenvalues at criticality (0.00128 and 0.0000698) are much smaller than in the first calculation (0.118 and 0.0525) which may explain the improved accuracy on .
Extensions of these methods for more states, more components and more dimensions are in progress. For practical purposes, the analytical methods discussed above need to be implemented numerically. We have succeeded to reproduce all the results obtained so far with adequate accuracy using numerical procedures which can be implemented using the most common programming languages. The fixed point can be found by monitoring successive bifurcations in tensor values. The high and low temperature phases are characterized by the fact that some tensor values go to zero (when is too small) or one (when is too large) if we iterate enough times. We can then fine-tune and observe the stabilization of tensors at some nontrivial values. The eigenvalues can be found by taking numerical derivatives of the one step transformation with respect to the initial values as close as possible to the nontrivial fixed point. This requires variations small enough but not too small since we have only limited accuracy on , however accurate results can be obtained using linear extrapolations to zero variation.
The formulation can be extended in dimension using tensors with with indices. The reason the TRG blocking works well in any dimension is that the links are orthogonal (dual) to domain boundaries. We keep the link variables across the boundaries of the block fixed and sum unrestrictedly over all the states inside the block. We have extended the third method described above to the 3D Ising model. The block is then a cube with four external legs coming out of each of the six faces. The blocking and the projection can be built out of the loop made by the four edges of a face with 4 external legs attached to each of the 4 corners. The initial transfer matrix can be obtained by tracing the external legs in the plane of the loop with their opposite leg. We then obtain a matrix corresponding to the 4 legs coming out of the plane of the loop in each direction. This matrix splits into two blocks. Numerically, we found for a two-state projection. The initial eigenvalues of at are 1.2325, 0.5082 and a pair with value 0.1682 which cannot be considered as small. For this reason, the value of for the 2D Ising model in the two-states approximation is not very close to the accurate value 0.630(2) pelissetto00b () but nevertheless more accurate than the MK approximation value 1.055.
There is good empirical evidence PhysRevB.86.045139 (); xiangprogress () supporting the idea that as we increase the number of states in TRG calculations, the numerical estimates of energy and entropy get closer to values obtained by exact methods (for the 2D Ising model) or Monte Carlo simulations (for the 2D model). However, in order to reach a reasonably large number of states (20-30) on a laptop, one needs to use efficient methods. The anisotropic TRG methods discussed above involve less contractions or external legs and the computational cost for states can be limited PhysRevB.86.045139 () to a growth in two dimensions and a growth in three dimensions. Numerical implementations of our method for up to 12 states are now in progress judahprogress (). In all cases, we observe a sharp split between a few large eigenvalues, their number being characteristic of the phase, and the small ones as in the two states case. The interesting behavior of these new maps will be discussed elsewhere judahprogress ().
Because of the binary nature of Grassmann numbers, the techniques developed for Ising models can also be used for lattice models with fermions. In the case where we have Grassmann variables with at every site and nearest neighbor interactions specified by a matrix at every link , we can write
The are linear combinations of the at the same site obtained from the decomposition with and unitary and diagonal with elements . The terms can then be factorized at every site and the local integrations performed. The states are now parametrized by and there are of them at every link. Translational invariance is essential to perform large volume calculations. For this reason, possible gauge interations would need to be averaged inside the blocks.
In summary, we have shown that two-state approximations of the TRG capture the universal behavior of Ising models much better than the MK approximation. Building on the numerical success of Ref. PhysRevB.86.045139 () which uses more states, we expect to be able to use the methods presented here to calculate the exponents of the 3D Ising model with unprecedented accuracy and study the analytical picture of the critical behavior provided in Ref. PhysRevD.86.025022 (); ElShowk:2012hu (). Recent numerical results for the TRG method applied to the model xiangprogress () suggest that improvements of the MK approximation could be applied to abelian models and resolve the controversy regarding the confinement in 4D gauge theory discussed in Ref. Tomboulis:2009zz (). We are hoping to be able to extend the TRG method for lattice gauge theories with fermions. Being able to block spin accurately would provide an efficient tool to study the continuum limits of asymptotically free models and the conformal window of models that could provide alternatives to the fundamental Higgs mechanism DeGrand:2010ba ().
Acknowledgments. Our work on the subject started while attending the KITPC workshop “Critical Properties of Lattice Models” in summer 2012. We thank M.C. Banuls, S. Chandrasekharan, A. Denbleyker, A. Li, Y. Liu, Y, M. Ogilvie, P. Orland, M. Qin, T. Tomboulis, J. Unmuth-Yockey, X-G Wen, T. Xiang, Z. Xie, J. Yu, and H. Zou for valuable conversations and suggestions. This research was supported in part by the Department of Energy under Contract No. FG02-91ER40664
- (1) J. L. Cardy Scaling and Renormalization in Statistical Physics, (Cambridge Univ. Pr., Cambridge, 1996)
- (2) Y. Meurice, R. Perry, and S.-W. Tsai, Phil. Trans. Roy. Soc. A369, 2602 (2011) and Refs. therein
- (3) K. Wilson and J. Kogut, Phys. Rep. 12, 75 (1974)
- (4) F. Verstraete and J. I. Cirac, cond-mat/0407066 (2004)
- (5) J. I. Cirac and F. Verstraete, J. of Phys. A 42, 4004 (2009)
- (6) S. R. White, Phys. Rev. Lett. 69, 2863 (1992)
- (7) U. Schollwoeck, Phil. Trans. R. Soc. A 369, 2643 (2011)
- (8) M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 (2007)
- (9) Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 103, 160601 (2009)
- (10) H. H. Zhao, Z. Y. Xie, Q. N. Chen, Z. C. Wei, J. W. Cai, and T. Xiang, Phys. Rev. B 81, 174411 (2010)
- (11) Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Phys. Rev. B 86, 045139 (2012)
- (12) A. A. Migdal, Sov.Phys.JETP 42, 743 (1975)
- (13) L. Kadanoff, Annals Phys. 100, 359 (1976)
- (14) G. Martinelli and G. Parisi, Nucl.Phys. B180, 201 (1981)
- (15) K. Wilson, Phys. Rev. B. 4, 3184 (1971)
- (16) G. Baker, Phys. Rev. B 5, 2622 (1972)
- (17) For a review see: Y. Meurice, J. Phys. A40, R39 (2007)
- (18) A. Denbleyker, Y. Liu, Y. Meurice, M. Qin, T. Xiang, Z. Xie, and J. Yu, preprint in progress
- (19) A. Pelissetto and E. Vicari, Phys. Rept. 368, 549 (2002) cond-mat/0012164
- (20) Y. Meurice, J. Unmuth-Yockey, and H. Zou, work in progress
- (21) S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin, and A. Vichi, Phys. Rev. D 86, 025022 (Jul 2012)
- (22) S. El-Showk and M. F. Paulos (2012), arXiv:1211.2810
- (23) E. T. Tomboulis, Mod. Phys. Lett. A 24, 2717 (2009)
- (24) T. DeGrand, Phil. Trans. R. Soc. A 369, 2701 (2011)