Sci. Rep. Kanazawa Univ.
Vol. 60, pp. 11–25, 2016

Computation of correlation functions by tensor
renormalization group method
Norihiro NAKAMOTO and Shinji TAKEDA⇤
Institute of Theoretical Physics, Kanazawa University, Kanazawa 920-1192, Japan(Received June 28, 2016 and accepted in revised form August 26, 2016)

Abstract We demonstrate how to compute the correlation function by using the tensor
renormalization group for 2D Ising model. From the correlation function around the criticaltemperature, we extract the correlation length and the spontaneous magnetization and thendetermine some critical exponents. Attractive feature of our method to extract the criticalexponents is that numerical derivative, that causes a loss of significant digits, and additional
simulation parameters, that demand more computational cost and complicated analysis, are
not required. On the other hand, our approach requires an additional treatment of impurity
tensors, but the additional cost is still tolerable.
Key words.
exponents

1

Tensor renormalization group, Correlation function, Phase transition, Critical

Introduction

The original idea of the tensor renormalization group (TRG) method was proposed by Levin andNave [1] for the Ising model on the triangular lattice. And then it has been applied to manyother models in two dimensional square lattice: f 4 theory [2], Schwinger model [3, 4] and finite
fermion density system [5]. A semi-analytic study for Ising model on square lattice was doneby Aoki et al., [6] who showed that the critical temperature and the critical exponent n can be
reproduced with a relatively high accuracy even though the bond dimension is restricted to only
two. Although the original idea was limited to the two dimensional system, a new coarse graining
method suited for any higher dimensional system was invented in Ref. [7]. The new method uses
the higher order singular value decomposition, thus it is called higher order tensor renormalizationgroup (HOTRG). The authors in [8] seriously analyzed the 3D Ising model by using the HOTRG
and obtained the critical temperature and the exponents precisely. The HOTRG was also appliedin two dimensional systems; XY model [9], O(3) model [10] and CP(1) model [11].
A key point of TRG method is to use the singular value decomposition (SVD) in the coarse
graining procedure, and then one naturally sets a bond dimension which controls the truncationin the decomposition. Thanks to this, one can reduce the size of tensor which would grow exponentially along the renormalization steps unless a truncation is applied, and can make the computational cost manageable while maintaining the important information of the system. It is in⇤ Corresponding author

Email: takeda@hep.kanazawa-u.ac.jp

12Computation of correlation functions by tensor renormalization group method

fact known that the TRG and HOTRG show very precisely results in the off-critical region while
quite large bond dimension is required near critical point to maintain the accuracy. To overcome
this difficulty, the tensor network renormalization (TNR) [12] and the loop optimization for tensor
network renormalization (Loop-TNR) [13] have been recently developed. These methods first recast the tensor network such that an entanglement is reduced and then a successive coarse graining
can be done very efficiently, namely the hierarchy of SVD becomes sharp and the required bonddimension can be reduced.
In this way, the efficient coarse graining methods even around critical point have been developed successfully. Next important issue is how to obtain the information of dynamics, for instancespectrum of theory, matrix elements, amplitudes by the TRG-related method. For that purpose,
one needs to compute the correlation function encoding the above information. The idea of how tocompute the correlation function using TRG was originally introduced in Ref. [14], but an actual
implementation of the algorithm has not been done yet for the classical 2D Ising model to ourbest knowledge1 . Therefore, in this paper we demonstrate the computation of spin-spin correlation function by using the algorithm for the model in details. And then we will extract the criticaltemperature and the critical exponents as a bench mark test. Actually, the extraction of the criticalexponents is sometimes a nontrivial task for TRG-related method since usually one needs to introduce numerical derivative or additional simulation parameters, say the external magnetic field,
in order to compute physical quantities. Since numerical results at some finite bond dimension isalready affected by systematic error, a further loss of significant digits due to the numerical derivative makes the situation much worse. An additional parameter simply requires a wider parameterregion to survey and then the cost and analysis tend to be hard. On the other hand, as we will see
later soon the computation of the correlation function requires only an additional treatment of an
impurity tensor and the total cost is only a few times larger than that of the partition function. Actually our method to extract the critical exponents is in sharp contract to the Monte Carlo method
in which the finite size scaling is often used. In fact, the TRG method, whose cost is proportional
to the logarithm of a system size, can easily deal with extremely large volume. Therefore one cansafely compute the correlation length even around critical point and the critical exponent n , which
dictates the scaling behavior, can be directly obtained from the correlation length. In this way, we
believe that this method is useful and enhances the utility of the TRG method.
The organization of this paper is as follows. In Sec. 2, a setup of the model, notations and tensor
network representation are provided. The computation of the correlation function is explained inSec. 3. In Sec. 4, we show numerical results. Our concluding remarks are summarized in Sec. 5.

2

Tensor network representation

The Hamiltonian for 2D Ising model is given byH=Â sa sb ,

(2.1)

ha,bi

where sa 2 {+1, 1} is a spin variable at a site a on 2D square lattice and the summation is takenover all nearest-neighbor lattice sites. In the following, we exclusively consider the vanishing
1 There were some calculations of the correlation function in Ref. [10, 15] with the TRG-related methods.

In
Ref. [10], spin correlation function is computed in HOTRG method. The authors in Ref. [15] studied the phase diagram of the quenched bond-diluted Ising model on the triangular lattice by using the TRG method. The long-distancebehavior of the correlation functions plays an important role in the study, but their algorithm to compute the correlation
function is not explicitly written and they did not refer to [14] in their paper.

13

Norihiro Nakamoto and Shinji Takeda

external magnetic field (h = 0). The partition function of the model is given byZ = Â e b H = Â eb Âha,bi sa sb ,
{s}

(2.2)

{s}

where b denotes the inverse temperature (T = 1/b ). To obtain the tensor network representationof the partition function, first we use an expansion for the Boltzmann weight1

eb sa sb = cosh b + sa sb sinh b = cosh b Â (sa sb tanh b )x .

(2.3)

x=0

After substituting the above expansion to the all Boltzmann weights in the partition function, one
can carry out the summation of spin variables exactly and then the partition function can be writtenasZ = 2V (cosh b )2V

Â ’ Ax y x y ,{x,y,x0 ,y0 } a

(2.4)

0 0
a a a awhere V is the lattice volume and the explicit form of the tensor is given byp
Ai jkl = ( tanh b )i+ j+k+l dmod(i+ j+k+l,2),0 .

(2.5)

Next, let us consider the one-point function for a single spin variable at a site a, namely themagnetization,m = hsa i = Â sa e b H /Z.

(2.6)

{s}

In a similar way to the partition function, one can consider a tensor network representation of the
one-point function by using the “impurity” tensor [14],pǍi jkl = ( tanh b )i+ j+k+l dmod(i+ j+k+l+1,2),0 ,(2.7)
where a slight difference appears only in the index of Kronecker delta. The numerator of the
one-point function is now expressed in terms of the impurity tensor,"#

Â sa e b H = 2V (cosh b )2V Â{x,y,x0 ,y0 }{s}Ǎxa ya xa0 y0a ’ Axc yc xc0 y0c .

(2.8)

c6=a

The two-point function of spin variables (a 6= b) is given by,
hsa sb i = Â sa sb e b H /Z = hs0 sr i,

(2.9)

{s}

where r = |a b| is a distance between two sites a and b. Similarly, the numerator of the two-point
function is also expressed in terms of the impurity tensors"#

Â sa sb e b H = 2V (cosh b )2V Â{s}{x,y,x0 ,y0 }Ǎxa ya xa0 y0a Ǎxb yb xb0 y0b ’ Axc yc xc0 y0c .
c6=a,b

(2.10)

14

Computation of correlation functions by tensor renormalization group method

Figure 1: Basic steps of coarse graining: (a) The decomposition of the normal tensor and (b)Producing the new tensor by contraction of the original indices.

3

Coarse graining for correlation function

Before explaining a coarse graining of the tensor network for the correlation function, first of all,we summarize the coarse graining for the partition function in eq.(2.4) which contains only thenormal tensor A. The normal coarse graining consists of two steps: the decomposition and the contraction [1]. Figure 1-(a) shows the decomposition step which uses singular value decomposition(SVD). There are two ways of decomposition:
3
Ai jkl = M(i⌦ j)(k⌦l) = Â U(i⌦ j)n Sn (V † )n(k⌦l) = Â Si1jn Skln,n

Si1jn = U(i⌦ j)n
and

p
Sn ,

3
⇤
Skln
= V(k⌦l)n

p
Sn ,

4
Ai jkl = M̃( j⌦k)(l⌦i) = Â Ũ( j⌦k)n S̃n (Ṽ † )n(l⌦i) = Â S2jkn Slin,n

S2jkn = Ũ( j⌦k)n

p
S̃n ,

(3.1)

n

(3.2)

(3.3)

n

4
⇤
Slin
= Ṽ(l⌦i)n

p
S̃n ,

(3.4)

where M and M̃ are matrix representation of the normal tensor, S and S̃ are diagonal matrix whose
diagonal elements are singular values, U, Ũ, V and Ṽ are unitary matrix and S1,2,3,4 are considered
as a rank-3 tensor. The new integer n is regarded as a new variable. In an actual implementation,
the summation of n is truncated and the bond dimension Dcut is introduced. Due to the truncation,
the rotational symmetry is broken thus S1,2,3,4 are considered as independent, although they shouldbe the same if the truncation is not applied. Figure 1-(b) shows the contraction step where thecontractions of original indices are carried out and the new tensor is produced
Anew
i jkl =

1
2
3 4
Sbal
Scbi
Sdc j .
Â Sadk(3.5)

a,b,c,dBy repeating the decomposition and the contraction steps iteratively, the coarse graining of thetensor network can be done as shown in Fig. 2, and then finally one can obtain the partition
function by contracting a small number of tensors.

Norihiro Nakamoto and Shinji Takeda

15

Figure 2: Coarse graining of tensor network. Step from 1 to 2 (3 to 4) is the decomposition. Step
from 2 to 3 (4 to 5) is the contraction.

Let us explain the coarse graining for the correlation functions [14], namely a coarse grainingof tensor network including the impurity tensor. As a concrete example, we address the one- andtwo-point function namely we consider two types of renormalization procedures in the following
First, we consider how to renormalize the impurity tensor for one spin variable whose procedure is shown in Fig. 3. We define “first sweep” which consists of four steps from 1-1 to 1-5 in
Fig. 3 where the linear lattice extent shrinks a half size. At the end of the first sweep, one obtains
three new impurity tensors. Similarly, the second sweep is composed of steps from 2-1 to 2-5 in
Fig. 3 and there are four new impurity tensors at the end of this sweep. On the other hand, in thethird sweep, the number of impurity tensors remains intact thus after the third sweep the impurity does not diffuse more than four sites and note that there is a cycle composed of steps from
3-1 to 3-5 in Fig. 3. This is a practically important feature to have a sustainable renormalization
procedure. Finally, one can perform a contraction for a small number of tensors and obtain the
numerator of the one-point function.
As the second example, Figure 4 shows the renormalization of the impurity tensor for two spin
variables. As seen in the previous paragraph, the number of impurity tensor for one spin variable
turns out to four after several sweeps. The initial configuration in Fig. 4-1-1 represents that the“two” impurity tensors for one spin variable after several sweeps are placed next to each other,thus in total there are eight impurity tensors in the network. At the end of the first sweep, thenumber of impurity tensor is reduced to six (see Fig. 4-1-5). In the second sweep, the numberof impurity tensor does not change any more (see Fig. 4-2-5). In contrast to the previous case,the cycle appears in the second sweep. The sustainable renormalization is also attained for thetwo-point function.

16

Computation of correlation functions by tensor renormalization group method

Figure 3: Renormalization of the impurity tensor for one spin variable. The open circle (white
point) represents the normal tensor. The colored point denotes the impurity tensor. In the firstand second sweep, the number of impurity tensor increases while it remains intact from the third
sweep and a cycle (from 3-1 to 3-5) is seen.

Norihiro Nakamoto and Shinji Takeda

17

Figure 4: Renormalization of the impurity tensor for two spin variables. The symbols are thesame as those of Fig. 3. The initial configuration represents the two impurity tensors for one spin
variable after several sweeps. In the first sweep, the number of impurity tensor decreases fromeight to six, while it remains intact from the second sweep and a cycle (from 2-1 to 2-5) is seen.

18

4Computation of correlation functions by tensor renormalization group method

ResultIn this section we show the numerical results of the correlation function in 2D Ising model andthen present the analysis to obtain the critical exponents.
1

TRG T=2.26700
Fitting point T=2.26700
Fitting function T=2.26700
TRG T=2.26933
Fitting point T=2.26933
Fitting function T=2.26933
TRG T=2.27600
Fitting point T=2.27600
Fitting function T=2.27600

0.9
0.8
0.7

G(r)

0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
20

25

210

215

220

r=2n
Figure 5: The correlation function as a function of distance r with fixed Dcut = 32 for 2D Isingmodel together with the fitting curve in eq.(4.4). Selected three temperatures are shown, T =2.26700 (the low temperature phase), 2.26933 (around exact Tc = 2.269185 · · · ), 2.27600 (thehigh temperature phase). The error bar is just the difference given in eq.(4.2).

We calculate the two-point function
G(r) = hs0 sr i,(4.1)

with the distance r = 20 , 21 , 22 , · · · , 220 on a fixed lattice volume V = L2 = 221 ⇥ 221 . The periodic
boundary condition is imposed on the system. The numerical results of the G(r) at fixed Dcut = 32
are shown in Fig. 5 for some selected temperatures. The error bar shown in the figure is just thedifference of the correlation function at Dcut = 24 and 32,

d = |hOiDcut =32

hOiDcut =24 |.

(4.2)

At higher temperature than the critical temperature Tc , the correlation function shows exponential
decay clearly. Such an exponential decay is also observed in the low temperature side T < Tc ,but it stays a nonzero constant for larger distance. This indicates a development of the nontrivial
spontaneous magnetization and we will address this issue later on. The error bars at the temperature in the off-critical region are very small, this means that the Dcut dependence is negligible.
On the other hand, around the critical point T = 2.26933, G(r) shows a long range correlation
as expected and Dcut dependence is clearly visible especially around r ⇡ 210 at this temperature.
Therefore one should be careful that the temperature should not be too close to the critical otherwise a large contamination due to the truncation error ruins the result.

Norihiro Nakamoto and Shinji Takeda

19

Let us see how to extract the magnetization and the correlation length from G(r). Around thecritical temperature T ⇡ Tc , the “connected” two-point function Gc (r) should behave as
Gc (r) = hsr s0 ic = hsr s0 i

hsr ihs0 i =A r/xe, for r ! •,
rE

(4.3)

where x is the correlation length and A is just a proportional constant. The exponent E is also afitting parameter which is supposed to be the critical exponent h = 1/4 at the critical temperature(x 1 = 0) for 2D Ising model. To compute the connected correlation function, the magnetizationm = hs0 i is required. Instead of using Gc (r), however, we fit G(r) with the following functional
form at each temperature,A
G(r) = C2 + E e r/x ,
r

(4.4)

where a constant term C2 , a fitting parameter, represents the tadpole contribution in eq.(4.3) andC is regarded as the magnetization2 . In the fitting procedure, we take eq.(4.2) as an error but notethat this does not mean that the systematic error originated from the finite Dcut is incorporated. The
fitting range of r is determined by the following way. Initially, a fit is done including all available
r data points. When the resulting chi square is large, say c 2 /d.o.f. > 10, we discard data point at
the smallest r and then perform a next fit. This procedure is repeated until we obtain a reasonablechi square of O(1) and we also check a stability of fitting parameters when changing the fit range.
Furthermore, we also monitor that the chosen lattice volume can afford a long range correlationx ⌧ L for all simulated temperatures to avoid the finite size effect on the correlation function.
Finally, let us show how to extract the critical exponents b , n and the critical temperature. As
mentioned, the fitting parameter C in eq.(4.4) is the spontaneous magnetization m, thus one canextract the b and Tc by fitting with the following temperature dependence3 ,C = m(T ) = D|T

Tc |b , for T < Tc .

(4.5)

The critical exponent n and the another estimate of Tc can be obtained from the correlation lengththat diverges at the critical point4 ,x (T ) = B|T

Tc | n .

(4.6)

The fitting should be done at higher and lower temperature regions separately. Figure 6 shows
the m, x and E as a function of the temperature together with the fitting line at Dcut = 32. The
magnetization m stays zero in the symmetric phase (the high temperature phase) while it has anonzero value in the broken phase (the low temperature phase). The correlation length x diverges
at the critical point. We observe that E converges to 1/4 near the critical temperature as expected.
The fitting range of temperature for m and x to obtain Tc and the critical exponents should be
carefully chosen such that it is not too close to the critical point, around which Dcut truncation error
is large, while the temperature should be in the scaling region. To find a proper region, we plot the
Tc and n in the high temperature phase as a function of the fitting range in Fig. 7 for Dcut = 32.
2 In fact, we directly compute the one-point function hs i by following the procedure given in the previous section0

and then it turns out that the value is exactly zero at a small number of the iteration but it takes significantly nonzero
value after several times of the iteration. Actually, the value of hs0 i is consistent with C in eq.(4.4) within errors foreach volume.
3 Practically, it is preferred to transfer to m1/bexp , where b
exp is a proper constant number (say, bexp = 1/8 for 2D
Ising model), and then fit with the form, m1/bexp µ |T Tc |b /bexp .
4 In the actual fitting procedure, it is useful to fit the inverse correlation length, x 1 µ |T
Tc |n .

20

Computation of correlation functions by tensor renormalization group method

Dcut
26
28
30
32
34
36
38
40
exact

Tc from x (T > Tc )2.26877(3)2.269030(7)2.269290(3)2.269330(8)2.26920(5)2.26935(1)2.269270(5)2.26917(4)2.2691853 · · ·n from x (T > Tc )0.904(11)0.977(4)0.945(6)0.931(6)1.035(21)0.997(11)0.967(7)1.013(16)1

Tc from m (T < Tc )2.26882(1)2.26903(1)2.269320(6)2.269330(6)2.26938(2)2.269440(6)2.26928(2)2.26924(1)2.2691853 · · ·b from m (T < Tc )0.1176(5)0.1179(8)0.1244(2)0.1253(3)0.1286(12)0.1270(2)0.1244(11)0.1249(4)0.125 = 1/8

Table 1: The resulting Tc and the critical exponents determined from the magnetization m and the
2pcorrelation length x for 2D Ising model. The exact critical temperature is Tc = log(1+.
2)

The x-axis of the figure is the value of the lowest temperature of the fitting range in which the hightemperature end is fixed. We take a representative value of n and Tc from the saturation region. We
repeat the same procedure for different Dcut = 26, 28, 30,..., 40. The resulting Tc and n for each
Dcut are summarized in Table. 1 and plotted in Fig. 9. We roughly observe that they converge tothe exact values for larger Dcut . A similar analysis is attempted for the low temperature phase. The
window analysis, however, turns out to be unstable since x (and also E) in the low temperaturephase is flattering as seen in Fig. 6, thus we abandon to present the result.
We also carry out the similar analysis to obtain b and Tc from the magnetization m in the
low temperature phase. We first check the finite size effect on the magnetization by performing
additional calculations on the larger lattice 231 ⇥ 231 , and then it turns out that it is completely
negligible in the temperature region we investigated. The window analysis for b and Tc is shownin Fig. 8 for Dcut = 32. The resulting values of Tc and b for various Dcut are given in Table. 1.
They tend to converge to the exact values for large Dcut as shown in Fig. 9.

5

Summary

In this paper, we have demonstrated how to compute the correlation function in the framework of
TRG. This technique opens up a possibility to gain some dynamical information of a system fromthe correlation function. Furthermore we have extracted the critical temperature and the universal
critical exponents from the correlation function for 2D Ising model, and have shown that they
are correctly reproduced. The method demonstrated here requires no numerical derivative andno additional parameters, thus it will be a powerful tool to determine the critical exponents formore complicated systems, where using a small Dcut is unavoidable due to a computational cost.
When using the method, however, one should be careful to ensure a proper window where thetemperature should not be too close to the critical point, but being in the scaling region.
So far, we have restricted ourself within the framework of TRG. As future work, it would beinteresting to establish a method to compute the correlation function in the framework of TNR and
Loop-TNR. This will open up further possibilities to capture more insights of dynamics in more
precise and less costly way.

Norihiro Nakamoto and Shinji Takeda

21

Acknowledgments
We would like to thank H. Kawauchi and R. Sakai for helpful advice and useful discussion. N. N.
is grateful to S. Kumamoto and S. Onai for useful conversation. This work was supported by
Kanazawa University SAKIGAKE Project.

References[1] M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 (2007)[2] Y. Shimizu, Mod. Phys. Lett. A 27, 1250035 (2012)[3] Y. Shimizu and Y. Kuramashi, Phys. Rev. D 90, 014508 (2014)[4] Y. Shimizu and Y. Kuramashi, Phys. Rev. D 90, 074503 (2014)[5] S. Takeda and Y. Yoshimura, PTEP 2015 (2015) no.4, 043B01[6] K. Aoki, T. Kobayashi and H. Tomita, Int. Jour. Mod. Phys. B 23, 18 (2009).
[7] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang and T. Xiang, Phys. Rev. B 86, 045139 (2012)[8] S. Wang, Z. Y. Xie, J. Chen, B. Normand and T. Xiang, Chin. Phys. Lett. 31, 070503 (2014)[9] J. F. Yu, Z. Y. Xie, Y. Meurice, Y. Liu, A. Denbleyker, H. Zou, M. P. Qin, J. Chen and T. Xiang, Phys. Rev. E 89,
013308 (2014)[10] J. Unmuth-Yockey, Y. Meurice, J. Osborn, H. Zou, arXiv:1411.4213.
[11] H. Kawauchi and S. Takeda, Phys. Rev. D 93, no. 11, 114503 (2016)[12] G. Evenbly and G. Vidal, Phys. Rev. Lett. 115, 180405 (2015)[13] S. Yang, Z. Gu, and X. Wen, arXiv:1512.04938v1[14] Z. Gu, M. Levin and X. Wen, Phys. Rev. B 78, 205116 (2008)[15] C. Gueven, M. Hinczewski and A. N. Berker, Phys. Rev. E 82, 051110 (2010)

22

Computation of correlation functions by tensor renormalization group method

0.014

Spontaneous magnetization for TRG
Fit [2.2645:2.2680]

0.012
0.01

m8

0.008
0.006
0.004
0.002
0
0.02

correlation length for TRG
Fit [2.2695:2.2780]

0.015

-1

0.01

0.005

0

-0.005

0.25

1.2
1

E

0.8
0.6
0.4
0.2

2.264

2.266

2.268

2.270

2.272

2.274

2.276

2.278

Temperature

Figure 6: The spontaneous magnetization m, the correlation length x and E as a function oftemperature around Tc in the 2D Ising model at Dcut = 32. The vertical yellow line represents
the location of the exact critical temperature Tc . Fitting lines, whose fitting range is given in thelegend, are also shown.

23

Norihiro Nakamoto and Shinji Takeda

Fit result
Exact critical temperature

2.27052.2703Tc

2.27012.26992.26972.26952.26932.2691
Fit result
=1

1.05
1
0.95
0.9
0.85
0.8
0.75
0.7
0.65

2.26952.27002.27052.27102.27152.2720Temperature

Figure 7: Fitting range dependence of the fitting parameters Tc and n in the high temperaturephase at Dcut = 32. The x-axis shows the lowest temperature in a fitting range where the hightemperature end is fixed to 2.278.

24Computation of correlation functions by tensor renormalization group method

2.27

Fit result
Exact critical temperature

2.2695Tc

2.269

2.2685

2.268

Fit result=0.125

0.14
0.13
0.12
0.11
0.1
0.09
0.08
2.26552.26602.26652.26702.26752.2680Temperature

Figure 8: Fitting range dependence of the fitting parameters Tc and b in the low temperature phaseat Dcut = 32. The x-axis for Tc and b shows the highest temperature in a fitting range where thelow temperature end is fixed to 2.2645.

25

Norihiro Nakamoto and Shinji Takeda

2.2694Tc

2.26922.269

2.2688

Tc at each Dcut fromTc at each Dcut from m
Exact critical temperature

1.06

at each Dcut from

=1

1.04
1.02
1
0.98
0.96
0.94
0.92
0.9
0.88
0.86
0.84
0.82
at each Dcut from m=0.125

0.129
0.127
0.125
0.123
0.121
0.119
0.117
20

25

30

35

Dcut

Figure 9: Dcut -dependence of Tc , n and b .

4045