Simulation of two-dimensional quantum systems on an infinite lattice revisited:
Corner transfer matrix for tensor contraction
Román Orús
*
and Guifré Vidal
School of Mathematics and Physics, The University of Queensland, Queensland 4072, Australia
共Received 20 May 2009; revised manuscript received 9 August 2009; published 8 September 2009
兲
An extension of the projected entangled-pair states 共PEPS兲 algorithm to infinite systems, known as the
iPEPS algorithm, was recently proposed to compute the ground state of quantum systems on an infinite two
dimensional lattice. Here we investigate a modification of the iPEPS algorithm, where the environment is
computed using the corner transfer matrix renormalization group 共CTMRG兲 method, instead of using one-
dimensional transfer matrix methods as in the original proposal. We describe a variant of the CTMRG that
addresses different directions of the lattice independently, and use it combined with imaginary time evolution
to compute the ground state of the two dimensional quantum Ising model. Near criticality, the modified iPEPS
algorithm is seen to provide a better estimation of the order parameter and correlators.
DOI: 10.1103/PhysRevB.80.094403 PACS number共s兲: 05.50.⫹q, 03.67.Mn, 03.65.Ud, 03.67.Hk
Understanding the emergent properties of many-body sys-
tems is one of the main goals of modern Physics. Projected
entangled-pair states 共PEPS兲
1
were recently proposed by Ver-
straete and Cirac to describe the ground state of finite, inho-
mogeneous quantum systems on a two dimensional 共2D兲 lat-
tice. A PEPS consists of a two dimensional network of
tensors whose coefficients are optimized so as to approxi-
mate the ground state of a local Hamiltonian. An extension
of the PEPS algorithm to infinite, homogeneous 2D lattices,
known as infinite PEPS 共iPEPS兲 algorithm,
2
was also subse-
quently developed. As a many-body ansatz, the infinite PEPS
had already been discussed by Sierra and Martín-Delgado
3
under the name of vertex matrix product ansatz, and had
been successfully used by Nishino and Okunishi,
4
under the
name of tensor product variational state, to evaluate the par-
tition function of a three-dimensional 共3D兲 classical system.
In the context of 2D quantum systems, a simplified version
of the ansatz, with only three free parameters, had also been
used in Ref. 5 prior to the iPEPS algorithm.
2
By considering
an evolution in imaginary time, the iPEPS algorithm over-
came the stability problems of previous proposals while al-
lowing for the optimization of all the coefficients in the an-
satz. Subsequently, ingenious approaches to optimize
homogeneous PEPS in finite systems with periodic boundary
conditions have also been proposed, such as those in Refs. 6
and 7, where the ansatz is often called tensor product state.
So far, algorithms based on the PEPS formalism have pro-
vided accurate ground-state properties of several models of
spins and hard-core bosons.
1,2,5–13
Importantly, they can ad-
dress systems beyond the reach of quantum Monte Carlo,
such as frustrated antiferromagnets.
12,13
An approximation 兩⌿ 典 to the ground state of a local
Hamiltonian H with an infinite PEPS is typically obtained
either by minimizing the expected value of the energy
具⌿兩H兩⌿典 or by simulating an evolution in imaginary time
兩⌿典⬇e
−H
兩⌿
0
典, where 兩⌿
0
典 is some initial state. In either
case, the tensors that define the infinite PEPS are optimized
iteratively and, in order to properly update a given tensor,
one needs to compute its environment: a 2D tensor network
that accounts for the rest of the tensors in the ansatz. Unfor-
tunately, computing the environment is hard and an approxi-
mation scheme is required. In the case of a homogeneous
system, several approximation schemes have been proposed:
共i兲 One-dimensional 共1D) transfer matrix approaches,
such as the infinite time-evolving block decimation 共iTEBD兲
approach
14
used in the original iPEPS algorithm.
2
共ii兲 2D coarse graining approaches based on the tensor
entanglement renormalization group 共TERG兲.
6,7
共iii兲 Corner transfer matrix 共CTM兲 approaches, such as
the corner transfer matrix renormalization group 共CTMRG兲
algorithm.
15
In this work we explore a modification of the iPEPS al-
gorithm where instead of computing the environment using
the iTEBD approach, as originally proposed in Ref. 2,we
use the CTMRG.
15
The CTM formalism was originally derived by Baxter
16
and later adapted by Nishino and Okunishi in the CTMRG
15
to numerically compute environments. The goal of this paper
is to investigate the performance of CTM approaches to
compute environments within the context of the iPEPS
algorithm.
2
Specifically, we first describe a versatile variant
of the CTMRG, the directional CTM approach, which ad-
dresses different directions of the lattice separately, and use it
in conjunction with imaginary time evolution to study the 2D
quantum Ising model near criticality. Accurate estimates of
the critical magnetic field and

exponent are obtained. Then
we compare the results obtained using the original iPEPS
algorithm 共where the environment was computed using the
iTEBD兲
14
with this new version of the algorithm. The modi-
fied iPEPS algorithm is seen to converge significantly faster
to the ground state and provide a better characterization of
the critical point.
17
Let us consider an infinite PEPS for the state 兩⌿典 of an
infinite 2D lattice L , which for concreteness we take to be a
square lattice, with each site labeled by two integers r
ជ
=共x , y兲 and represented by a complex vector space V of di-
mension d. In the simplest scenario, the infinite PEPS is
characterized by a single tensor A that is repeated on all
lattice sites. It has components A
s udlr
, where s labels a local
共physical兲 basis of V 共s =1,¯ ,d兲 and u ,d ,l ,r are bond indi-
ces ranging from 1 to D, with D the bond dimension of the
infinite PEPS, see Fig. 1共a兲. Let a denote the reduced tensor
PHYSICAL REVIEW B 80, 094403 共2009兲
1098-0121/2009/80共9兲/094403共4兲 ©2009 The American Physical Society094403-1
a⬅兺
s=1
d
A
s
丢 A
s
ⴱ
, with double bond indices such as u
¯
=共u ,u
⬘
兲关see Fig. 1共b兲兴. Then the scalar product 具⌿ 兩⌿典 can
be expressed as a two dimensional network E made of infi-
nitely many copies of a, Fig. 1共c兲. 共Notice that with a proper
choice of tensor a, E can also represent the partition function
of a 2D classical statistical model兲. The environment of site
r
ជ
, E
关r
ជ
兴
=
E /
a
关r
ជ
兴
, is obtained from E by removing the tensor a
on site r
ជ
, see Fig. 1共d兲. The goal of the CTMRG algorithm is
to compute an approximation G
关r
ជ
兴
to E
关r
ជ
兴
by finding the fixed
point of four CTMs. This effective environment is given
in terms of a small tensor network, G
关r
ជ
兴
=兵C
1
,T
1
,C
2
,T
2
,
C
3
,T
3
,C
4
,T
4
其, where tensors C
1
,C
2
,C
3
,C
4
represent four
CTMs 共one for each corner兲, and tensors T
1
,T
2
,T
3
,T
4
repre-
sent two half-column and two half-row transfer matrices,
Fig. 1共e兲.
In the directional variant of the CTMRG used in this
work, the eight tensors of G
关r兴
are updated according to four
directional coarse-graining moves, namely left, right, up,
and down moves, which are iterated until the environment
converges. Given an effective environment G
关r
ជ
兴
=兵C
1
,T
1
,C
2
,T
2
,C
3
,T
3
,C
4
,T
4
其, a move, e.g., to the left, con-
sists of the following three main steps, Fig. 2:
共1兲 Insertion: insert a new column made of tensors T
1
, a,
and T
3
as in Fig. 2共b兲.
共2兲 Absorption: contract tensors C
1
and T
1
, tensors C
4
and
T
3
, and also tensors T
4
and a, resulting in two new CTMs C
˜
1
and C
˜
4
, and a new half-row transfer matrix T
˜
4
, see Fig. 2共c兲.
共3兲 Renormalization: truncate the vertical indices of C
˜
1
,
T
˜
4
, and C
˜
4
by inserting the isometry Z, Z
†
Z=I. This produces
renormalized CTM’s C
1
⬘
=Z
†
C
˜
1
, C
4
⬘
=C
˜
4
Z and half-row trans-
fer matrix T
4
⬘
, Figs. 2共d兲 and 2共e兲.
A proper choice of isometry Z in the renormalization step
is of great importance. One possibility is to use the eigen-
value decomposition of the product of the four CTMs C
˜
1
, C
2
,
C
3
, and C
˜
4
as in Ref. 15. Here we consider instead the ei-
genvalue decomposition of C
˜
1
C
˜
1
†
+C
˜
4
†
C
˜
4
=Z
˜
D
Z
Z
˜
†
, Fig. 2共f兲,
and use the isometry Z that results from keeping the entries
of Z
˜
corresponding to the
largest eigenvalues of D
Z
. Simi-
larly to Ref. 15, this isometry targets the CTMs of the effec-
tive environment instead of the wave function itself. The cost
of implementing these steps scales with D and
as O共D
6
3
兲.
The net result is a new effective environment G
⬘
关r
ជ
兴
for site r
ជ
given by tensors 兵C
1
⬘
,T
1
,C
2
,T
2
,C
3
,T
3
,C
4
⬘
,T
4
⬘
其, see Fig. 2共d兲.
By composing the four moves of the directional CTM we
recover one iteration of CTMRG.
15
The additional flexibility
provided by individual moves can be used to accelerate con-
vergence in a specific direction, e.g., in highly anisotropic
systems. In addition, the prescription used to compute the
isometry Z is still valid—and produces stable results
17
—in
the context of simulating imaginary time evolution described
in this work. As with other similar methods,
14,15,19
an imme-
diate application of the directional CTM is to compute ex-
pected values from 2D classical partition functions 共results
not shown兲.
In order to compute the ground state of 2D quantum mod-
els by simulating imaginary time evolution, we consider an
infinite PEPS characterized by two tensors A and B. The
optimization of tensors A and B proceeds in the same way as
we proposed as part of the iPEPS algorithm,
2
but with the
crucial difference that here the required environment for two
contiguous sites is computed with the directional CTM in-
stead of using 1D transfer matrix 共1DTM兲 techniques. The
scalar product 具⌿兩⌿典 consists of an infinite 2D tensor
network made of copies of the reduced tensors a and b.We
first consider the environment E
关r
ជ
1
,r
ជ
2
,r
ជ
3
,r
ជ
4
兴
of a four-site unit
cell, see Fig. 3共a兲, and approximate it with an effective
environment G
关r
ជ
1
,r
ជ
2
,r
ជ
3
,r
ជ
4
兴
=兵C
1
,T
b1
,T
a1
,C
2
,T
a2
,T
b2
,C
3
,T
b3
,
T
a3
,C
4
,T
a4
,T
b4
其 made of twelve tensors, which are com-
puted by iterating left, right, up and down moves. These
directional moves are a natural adaptation to a four-tensor
unit cell of those in Fig. 2. As shown in Fig. 3共b兲, the half-
row and half-column transfer matrices T
i
共i =1,2,3,4兲 are re-
placed with pairs of half-row and half-column transfer ma-
FIG. 1. 共Color online兲 Diagrammatic representation of 共a兲 infi-
nite PEPS tensor A with physical index s and bond indices u ,r,d
and l; 共b兲 reduced tensor a; 共c兲 infinite 2D tensor network E ; 共d兲
environment E
关r
ជ
兴
for site r
ជ
; 共e兲 eight-tensor effective environment
G
关r
ជ
兴
.
FIG. 2. 共Color online兲共a兲–共d兲 Main steps of a left move: inser-
tion, absorption and renormalization; 共e兲 the CTMs C
˜
1
, C
˜
4
and the
half-row transfer matrix T
˜
4
are renormalized with isommetry Z; 共f兲
eigenvalue decomposition for the sum of the squares of CTMs C
˜
1
and C
˜
4
.
ROMÁN ORÚS AND GUIFRÉ VIDAL PHYSICAL REVIEW B 80, 094403 共2009兲
094403-2
trices T
ai
,T
bi
. This time, in order to implement, e.g., a left
move, two new columns are inserted in the system in step 共1兲,
see Fig. 3共c兲. For each of the inserted columns, we perform
the absorption and renormalization steps 共2兲 and 共3兲. The
renormalization step requires introducing an additional isom-
etry W, which we compute in an analogous way as isometry
Z, see Fig. 3共e兲. As before, the cost of a move scales as
O共D
6
3
兲. Finally, from a converged environment for the
four-site unit cell, an effective environment for any pair of
nearest-neighbor sites is easily obtained with an additional
directional move.
To demonstrate the performance of the approach, we have
computed an infinite PEPS approximation to the ground state
of the spin-1/2 quantum Ising model on a transverse mag-
netic field, H
I
共兲=−兺
具r
ជ
,r
ជ
⬘
典
z
关r
ជ
⬘
兴
z
关r
ជ
⬘
兴
−兺
r
ជ
x
关r
ជ
兴
, by simulating
an imaginary time evolution. The simulation proceeds as in
Ref. 2, but we use the directional CTM to obtain the effec-
tive environment at each step of the imaginary time evolu-
tion. This evolution is performed with decreasing time steps
␦
ranging from 10
−1
to 10
−5
, and until convergence of local
observables and two-point correlators is attained.
Figure 4 shows the order-parameter m
z
⬅具⌿兩
z
兩⌿典 as a
function of the transverse magnetic field , for 共D,
兲 equal
to 共2,20兲 and 共3,30兲, where the value of
is chosen so that
the results are converged with respect to this parameter. Re-
markably, an infinite PEPS with bond dimension D =3 al-
ready produces results within less than a percent from the
best quantum Monte Carlo estimates for the critical magnetic
field
c
MC
⬇3.044 and critical exponent for the order-
parameter

MC
⬇0.327,
20
namely with relative errors ⬇0.1%
and ⬇0.3% respectively. Table I contains a comparison with
results obtained with the original version of the iPEPS algo-
rithm and with the TERG algorithm for large systems.
6
It is particularly instructive to compare the performance
of the original and present versions of the iPEPS algorithm,
since they are both based on imaginary time evolution and
only differ in how the two-site environment is computed: by
means of the iTEBD and directional CTM approaches, re-
spectively. One finds that when computing environments us-
ing the directional CTM, a significantly better infinite PEPS
approximation to near-critical ground states is obtained, lead-
ing to a more accurate characterization of the quantum phase
transition. As shown in Fig. 5, the resulting infinite PEPS
also displays stronger correlators S
zz
共l兲⬅具⌿兩
z
关r
ជ
兴
z
关r
ជ
+le
ˆ
x
兴
兩⌿典
−共m
z
兲
2
.
However, further comparison of results involving also
other spin models reveals that, away from the quantum criti-
cal point, both the directional CTM and the iTEBD ap-
proaches yield equivalent accuracies for ground-state prop-
erties. In particular, both versions of the iPEPS algorithm are
equally suited to study first-order phase transitions, a task for
TABLE I. Critical-point
c
and exponent

for the 2D quantum
Ising model as estimated by the new and old versions of the iPEPS
algorithm, as well as the TERG 共for a finite lattice of up to
2
9
⫻2
9
spins兲. For reference, the quantum Monte Carlo estimation
is
c
MC
⬇3.044 and

MC
⬇0.327 共Ref. 20兲.
iPEPS with
directional CTM
iPEPS with
iTEBD
a
TERG
b
c
D=2 3.08 3.10 3.08
c
D=3 3.04 3.06

D=2 0.333 0.346 0.333

D=3 0.328 0.332
a
Reference 2.
b
Reference 6.
FIG. 3. 共Color online兲共a兲 Environment of the four-site unit cell;
共b兲 twelve-tensor effective environment; 共c兲 two new columns are
inserted, and absorbed toward the left and renormalized individu-
ally. The diagram shows the contraction leading to C
˜
1
, C
˜
4
, T
˜
a4
, and
T
˜
b4
when absorbing the first column, and also to the CTMs Q
˜
1
and
Q
˜
4
; 共d兲 two isometries Z and W are used to obtain the renormalized
half-row transfer matrices T
a4
⬘
and T
b4
⬘
; 共e兲 eigenvalue decomposi-
tion for the sum of squares of CTMs Q
˜
1
and Q
˜
4
.
2.7 2.75 2.8 2.85 2.9 2.95 3 3.05 3.1 3.15 3.2 3.25 3.3
0
0.1
0.2
0.3
0.4
0.5
0.6
m
z
−5 −4 −3 −2 −1
−2
−1.5
−1
−0.5
ln|−
c
|
ln(m
z
)
2.95 2.975 3 3.025 3.05
0
0.1
0.2
0.3
0.4
0.5
m
Z
iTEBD D=3
iTEBD D=2
directional CTM D=3
directional CTM D=2
D=2,
c
= 3.08,
= 0.333 ± 0.002
D=3,
c
= 3.041,
= 0.328 ± 0.005
FIG. 4. 共Color online兲 Order-parameter m
z
as a function of the
transverse field , computed with the directional CTM approach.
Lines are a guide to the eye. The lower-left inset shows a log plot
共in natural logarithms兲 of m
z
versus 兩 −
c
兩, including our estimates
for
c
and

. The continuous lines show the linear fits. The upper-
right inset shows a comparison close to criticality with the results
from Ref. 2 using the original version of the iPEPS algorithm,
which used iTEBD 共dashed lines兲. Results correspond to 共D ,
兲
equal to 共2,20兲 and 共3,30兲.
SIMULATION OF TWO-DIMENSIONAL QUANTUM SYSTEMS… PHYSICAL REVIEW B 80, 094403 共2009兲
094403-3
which they are particularly successful.
9,13
Indistinguishable
results are also obtained in models with long or infinite cor-
relation lengths, such as the Heisenberg antiferromagnet on a
square lattice,
13
when using a bond dimension D that is too
small to offer a proper approximation to the ground state 共in
Ref. 13, a bond dimension D =5 still produces a spontaneous
magnetization that is off by 10%兲. It seems, therefore, that
computing environments with a CTM approach leads to bet-
ter results when the following two requirements are simulta-
neously met: 共i兲 the ground state must have a long correla-
tion length 共e.g., near a quantum critical point兲, and 共ii兲 the
bond dimension D must be sufficiently large that the ansatz
can in principle properly approximate the ground state.
The correlators in Fig. 5 also show that numerical calcu-
lations with infinite PEPS introduce an artificial finite corre-
lation length at criticality. This is a consequence of the trun-
cation in parameter
in the calculations of the effective
environments. This truncation effect is similar to the one
discussed in Ref. 21 and could in principle be analyzed in a
systematic way.
To summarize, PEPS are a valuable ansatz to approximate
the ground state of 2D lattice models, as previously demon-
strated by several authors
1,2,5–13
with a variety of systems of
interacting spins and hard-core bosons, including frustrated
spins
12,13
that cannot be addressed with quantum Monte
Carlo techniques. Several methods have been proposed to
optimize the PEPS. The main factor limiting the accuracy of
the results is the bond dimension D. However, for a fixed
value of the bond dimension D, the quality of the approxi-
mation may also depend on the method used to optimize the
ansatz. In this work we have investigated a modification of
the iPEPS algorithm, where the iTEBD of the original pro-
posal has been replaced with a directional CTM, a variant of
CTMRG.
15
The new version of the algorithm provides a sig-
nificantly better description of the ground state near critical-
ity.
Support from the University of Queensland 共Grant No.
ECR2007002059兲 and the Australian Research Council
共Grants No. FF0668731, and No. DP0878830兲 is acknowl-
edged.
*
orus@physics.uq.edu.au
1
F. Verstraete and J. Cirac, arXiv:cond-mat/0407066 共unpub-
lished兲; V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. A 75,
033605 共2007兲.
2
J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys.
Rev. Lett. 101, 250602 共2008兲.
3
G. Sierra and M. Martín-Delgado, arXiv:cond-mat/9811170 共un-
published兲.
4
T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 67, 3066 共1998兲.
5
Y. Nishio, N. Maeshima, A. Gendiar, and T. Nishino,
arXiv:cond-mat/0401115 共unpublished兲.
6
Z.-Cheng Gu, M. Levin, and X.-Gang Wen, Phys. Rev. B 78,
205116 共2008兲.
7
H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101,
090603 共2008兲.
8
A. Isacsson and O. F. Syljuasen, Phys. Rev. E 74, 026701
共2006兲; H-Q. Zhou, R. Orús, and G. Vidal, Phys. Rev. Lett. 100,
080601 共2008兲.
9
R. Orús, A. C. Doherty, and G. Vidal, Phys. Rev. Lett. 102,
077203 共2009兲.
10
B. Li, S. Li, and H. Zhou, Phys. Rev. E 79, 060101共R兲共2009兲.
11
J. Jordan, R. Orús, and G. Vidal, Phys. Rev. B 79, 174515
共2009兲.
12
V. Murg, F. Verstraete, and J. Cirac, Phys. Rev. B 79, 195119
共2009兲.
13
B. Bauer, G. Vidal, and M. Troyer, arXiv:0905.4880 共unpub-
lished兲.
14
G. Vidal, Phys. Rev. Lett. 98, 070201 共2007兲; R. Orús and G.
Vidal, Phys. Rev. B 78, 155117 共2008兲.
15
T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 65, 891 共1996兲.
16
R. J. Baxter, J. Math. Phys. 9, 650 共1968兲; J. Stat. Phys. 19, 461
共1978兲; Exactly Solved Models in Statistical Mechanics 共Aca-
demic Press, London, 1982兲.
17
The accuracies obtained in this work for the critical point of the
2D quantum Ising model are comparable to those obtained in
Ref. 18 for the 3D classical Ising model. The present approach,
also valid for 2D quantum systems without 3D classical coun-
terpart, does not suffer from the instabilities reported in Ref. 18.
18
T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima, and Y. Akutsu,
Prog. Theor. Phys. 105, 409 共2001兲; A. Gendiar, N. Maeshima,
and T. Nishino, ibid. 110, 691 共2003兲; A. Gendiar, T. Nishino,
and R. Derian, Acta Phys. Slov. 55, 141 共2005兲.
19
M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 共2007兲.
20
H. W. J. Blöte and Y. Deng, Phys. Rev. E 66, 066110 共2002兲.
21
L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I. Latorre,
Phys. Rev. B 78, 024410 共2008兲.
1 2 3 4 5 6 7 8 9 10 20 30
10
−1
10
−2
10
−3
10
−4
10
−5
l
S
zz
(l)
iTEBD D=3
iTEBD D=2
directional CTM D=3
directional CTM D=2
FIG. 5. 共Color online兲 Log plot of the correlator S
zz
共l兲共in base
10兲, as computed with the original and present versions of the iP-
EPS algorithm, namely using the iTEBD 共dashed lines兲 and the
directional CTM 共solid lines兲 approaches. Lines are a guide to the
eye. Our results are for the 共D,
兲 pairs 共2,20兲共for =3.09兲 and
共3,30兲共for =3.04兲.
ROMÁN ORÚS AND GUIFRÉ VIDAL PHYSICAL REVIEW B 80, 094403 共2009兲
094403-4
