Papers in Physics, vol. 8, art. 080006 (2016)

Received: 2 August 2016, Accepted: 12 October 2016
Edited by: K. Hallberg
Licence: Creative Commons Attribution 3.0
DOI: http://dx.doi.org/10.4279/PIP.080006

www.papersinphysics.org

ISSN 1852-4249

An efficient density matrix renormalization group algorithm for chains
with periodic boundary condition

Dayasindhu Dey,1 Debasmita Maiti,1 Manoranjan Kumar1∗

The Density Matrix Renormalization Group (DMRG) is a state-of-the-art numerical tech-
nique for a one dimensional quantum many-body system; but calculating accurate results
for a system with Periodic Boundary Condition (PBC) from the conventional DMRG has
been a challenging job from the inception of DMRG. The recent development of the Matrix
Product State (MPS) algorithm gives a new approach to find accurate results for the one
dimensional PBC system. The most efficient implementation of the MPS algorithm can
scale as O(p × m3), where p can vary from 4 to m2. In this paper, we propose a new
DMRG algorithm, which is very similar to the conventional DMRG and gives comparable
accuracy to that of MPS. The computation effort of the new algorithm goes as O(m3) and
the conventional DMRG code can be easily modified for the new algorithm.

I. Introduction

The quantum many-body effect in the condensed
matter gives rise to many exotic states such as su-
perconductivity [1], multipolar phases [2,3], valence
bond state [4], vector chiral phase [2,5] and topolog-
ical superconductivity [6]. These effects are promi-
nent in the one dimensional (1D) electronic systems
due to the confinement of electrons. The confine-
ment of electrons and the competition between the
electron-electron repulsion and the kinetic energies
of electrons produce many interesting phases like
Spin Density Wave (SDW), dimer or the bond or-
der wave phase and Charge Density Wave (CDW)
phase in one dimensional systems [7–9]. Although
these quantum many-body effects in the system are
crucial for exotic phases, dealing with these systems
is a challenging job because of the large degrees of

∗E-mail: manoranjan.kumar@bose.res.in

1 S. N. Bose National Centre for Basic Sciences, Block JD,
Sector III, Salt Lake, Kolkata 700106, India.

freedom. The degrees of freedom increase as 2N

or 4N for a spin-1/2 system or a fermionic system,
respectively.

In most cases, the exact solutions for these sys-
tems with large degrees of freedom are almost im-
possible to reach. Therefore, during the last three
decades many numerical techniques have been de-
veloped, e.g., Quantum Monte Carlo (QMC) [10],
Density Functional Theory (DFT) [11], Renormal-
ization Group (RG) [12] and Density Matrix Renor-
malization Group (DMRG) method [13, 14]. The
DMRG is a state-of-the-art numerical technique
for 1D systems with Open Boundary Condition
(OBC). However, the numerical effort to maintain
the accuracy for PBC systems becomes exponen-
tial [15, 16]. It is well known that the Periodic
Boundary Condition (PBC) is essential to get rid
of the boundary effect of a finite open chain and
also to preserve the inversion symmetry in the sys-
tems [17].

The DMRG technique is based on the system-
atic truncation of irrelevant degrees of freedom and
has been reviewed extensively in Ref. [15, 16]. In a

080006-1



Papers in Physics, vol. 8, art. 080006 (2016) / D. Deyet al.

1D system with OBC, the number of relevant de-
grees of freedom is small [15, 16]. Let us consider
that for a given accuracy of the OBC system, m

obc

number of eigenvectors of the density matrix is re-
quired, then the conventional DMRG for the PBC
system requires O(m2obc) [18]. In the conventional
DMRG, computational effort for the OBC systems
with sparse matrices goes as O(m3obc), whereas, it
goes as O(m6obc) for the PBC system [19]. The
accuracy of energies for the PBC systems calcu-
lated from the conventional DMRG decreases sig-
nificantly, and there is a long bond in the system
which connects both ends.

The accuracy of operators decreases with the
number of renormalization, especially the rais-
ing/creation S+/a+ and lowering/annihilation
S−/a− operator of spin/fermionic systems. The
conventional DMRG is solved in a Sz basis, there-
fore the exact Sz operator remains diagonal and,
multiple times, renormalization deteriorates the ac-
curacy slowly; but S+/a+ and S−/a− are off diag-
onal in this exact basis and therefore, the multiple
time renormalization of these operators decreases
the accuracy of the operators. A similar type of ac-
curacy problems occurs for multiple time renormal-
ized a+ and a− in the fermionic systems. In fact, it
has been noted that accuracy of energy of a system
with PBC significantly increases if the superblock
is constructed with very few times renormalized op-
erators [9]. To avoid multiple renormalization, new
sites are added at both ends of the chain in such a
way that only second time renormalized operators
are used to construct the superblock. In this algo-
rithm, there is a connectivity between the old-old
sites and their operators are renormalized; and this
connectivity spoils the sparsity of the superblock
Hamiltonian [9].

In this paper, a new DMRG algorithm is pro-
posed, which can be implemented upon the existing
conventional DMRG code in a few hours and gives
accurate results which are comparable to those of
MPS algorithm. In fact, this algorithm can be im-
plemented for two-legged ladders without much ef-
fort [20]. We have studied the spectrum of den-
sity matrix of the system block, ground state en-
ergy and correlation functions of a Heisenberg anti-
ferromagnetic Hamiltonian for spin-1/2 and spin-1
on a 1D chain with PBC.

This paper is divided into five sections. The
model Hamiltonian is discussed in section II. The

new algorithm and the comparative studies of al-
gorithms are done in section III. The accuracy of
various quantities is studied in the section IV. In
section V, results and algorithm are discussed.

II. Model Hamiltonian

Let us consider a strongly correlated electronic sys-
tem where Coulomb repulsion is dominant, there-
fore the charge degree of freedom gets localized, for
example, Hubbard model in large U limit in a half
filled band. In this limit, the system becomes in-
sulating, but the electrons can still exchange their
spin. The Heisenberg model is one of the most
celebrated models in this limit, and only the spin
degrees of freedom are active in the model. The
Heisenberg model Hamiltonian can be written as

H =
∑
〈ij〉

Jij ~Si · ~Sj (1)

where, Jij is the anti-ferromagnetic exchange in-
teraction between the nearest neighbor spin. In the
rest of the paper, Jij is set to be 1.

III. Comparison of algorithms

A ground state wavefunction calculated from the
conventional DMRG can be represented in terms
of the Matrix Product State (MPS), as shown by
Ostlund and Rommer [21]. The wavefunction can
be written as

|ψ〉 =
∑

n1,n2,...,nL

tr(A1n1A
2
n2
. . .ALnL )|n1,n2, . . . ,nL〉,

(2)

where Aknk are a set of matrices of dimension m×m
for site k and with nk degrees of freedom. The wave
function |ψ〉 can be accurately found if m is suffi-
ciently large. The expectation value of an operator
Ok in the gs [19, 22] can be written as

〈Ok〉 = tr


 ∑
nk,n

′
k

〈nk|Ok|n
′

k〉A
k
nk
⊗Ak

n
′
k


 , (3)

080006-2



Papers in Physics, vol. 8, art. 080006 (2016) / D. Deyet al.

where nk is the local degrees of freedom of site k.
The matrix A can be evaluated by using the follow-
ing equations

Hk|ψk〉 = λNk|ψk〉, (4)
Nk = Ek+11 E

k
1 . . .E

L
1 E

1
1 . . .E

k−1
1 , (5)

where

Ek1 =
∑
nk,n

′
k

〈nk|1|n
′

k〉A
k
nk
⊗Ak

n
′
k

. (6)

Here, Hk is the effective Hamiltonian of k
th site

and λ is the expectation value of energy. The A
matrices are evaluated at this point and the ma-
trices are rearranged to keep the algorithm stable.
The Hk and N can be calculated recursively while
evaluating A, one site at a time [19]. Here, N ma-
trices are ill-conditioned and require storing, ap-
proximately L2 matrices as well as multiplication
of L2 matrices of m×m size [19] at each step. The
evaluation of A and N are done for all the sites and
backward and forward sweeps for all the sites are
executed similarly to the finite system DMRG. The
mathematical operations of matrices of dimension
m2 ×m2 represent the Hamiltonian cost ∼ (o)m6,
but the special form of these matrices reduces the
cost by a factor of m. Therefore, this algorithm
scales as ∼ (o)m5 [19].

The above algorithm is extended by Verstraete
et. al., for translational invariant systems [23].
Only two types of matrices, A1 and A2, are consid-
ered [23]. Product of the two matrices can be re-
peated to compute N. In this algorithm, only two
matrices, A1 and A2, are updated and optimized
to get the gs properties. This algorithm scales as
(o)m3, although it does not work for a finite sys-
tem, or systems with impurity, etc. Pippan et. al.
introduced another MPS based efficient algorithm
for translational invariant PBC systems [19]. In the
old version of MPS, most of the computation cost
goes to constructing the product of m2×m2 matri-
ces E defined in Eq. (6). The new MPS algorithm
overcomes this problem by performing a SVD of
the product of sufficiently large (l � 1) number
of m2 × m2 transfer matrices[19, 28]. The singu-
lar values, in general, decay very fast; therefore,
only p (� m2) among m2 singular values are sig-
nificant [28]. Thus, the computational cost now is

reduced to (o) p×m3 [19]. However, one requires
p ∼ m to reach adequate numerical accuracy of
physical measures, as pointed out in Ref. [28].

Although the above technique is efficient and ac-
curate, there are various reasons for developing the
new algorithm. First, the modified MPS works ef-
ficiently for a system where the singular values of
products of matrices decay exponentially and this
algorithm scales as (o) pm3, where p can vary lin-
early with m. Second, the implementation of the
MPS based numerical technique is quite different
from the conventional DMRG, and the MPS algo-
rithm should be written from scratch. Third, many
conventional numerical techniques like the dynami-
cal correction vector [24] or continued fraction [25],
and implementation of symmetries like parity or in-
version symmetries are difficult. In this paper, we
will explain a new algorithm which is very similar to
the conventional DMRG technique, and also show
that the new algorithm can give accuracy compa-
rable to that of MPS based techniques. This al-
gorithm is applied for S = 1/2 and S = 1 chains
with PBC. But first, let us try to understand the
algorithm before discussing the results.

In this algorithm, we will try to avoid the multi-
ple renormalization of operators, whereas the other
parts of the algorithm remain the same as the con-
ventional DMRG. Before going to the new algo-
rithm, let us recap the conventional DMRG.

1. Start with a superblock of four sites consisting
of one site for both the left and the right block
and two new sites.

2. Get desired eigenvalues and eigenvectors of the
superblock and construct the density matrix ρ
of the system which consists of the left or the
right block and a new site.

3. Now, construct an effective ρ̃ with m num-
ber of eigenvectors of ρ, corresponding to the
m largest eigenvalues. The effective system
Hamiltonian and all operators in the truncated
basis are constructed using the following equa-
tions:

H = ρ̃†Hρ̃; O = ρ̃†Hρ̃ (7)

4. Superblock is constructed using the effective
Hamiltonian and operators of the system block
and two new sites.

080006-3



Papers in Physics, vol. 8, art. 080006 (2016) / D. Deyet al.

(a)  (b) 

(c) 

Figure 1: Pictorial representation of the new
DMRG algorithm with only one site in the new
block. (a) One starts with two blocks left and right
represented by filled circles and two new sites blocks
represented as open circles. The dotted box repre-
sents the system block for the next step. (b) Su-
perblock of the next DMRG step is shown. (c) The
final step of infinite DMRG of N = 4N system size
is shown with 2N−1 number of sites in each of the
left and the right blocks and two new sites.

5. Repeat all the steps from 2 until the desired
system size is reached. The full process is
called infinite DMRG.

As mentioned earlier, the conventional algorithm
is excellent for a 1D open chain as superblock is
constructed with only one time renormalized op-
erators. However, for a PBC system, one needs a
long bond; therefore, at least two operators of su-
perblocks are renormalized multiple times. In the
new algorithm, the multiple time renormalization
of operator is avoided and the algorithm goes as:

1. Start with a superblock with four blocks con-
sisting of a left and a right block and two new
site blocks. The blocks are shown in Fig. 1
as filled circles and may have more than one
site. Here, we have considered only one site in
each block. New blocks may also have more
than one site and are shown as open circles.
In this paper, new blocks have one site in a
chain or two for a ladder, like a structure with
PBC [17].

2. Get the eigenvalues and eigenvectors of the su-

perblock and construct the density matrix ρ of
the system which consists of the left or the
right block and two new blocks. The left sys-
tem block is shown inside the box in Fig. 1a.

3. Now, construct an effective ρ̃ with m num-
ber of eigenvectors of ρ corresponding to the
m largest eigenvalues of the density matrix.
The effective system Hamiltonian and opera-
tors in the truncated basis are calculated using
Eq. (7).

4. The superblock is constructed using the ef-
fective Hamiltonian and operators of system
blocks and two new sites.

5. Now, go to step 2 and the processes 2 – 5 are
repeated till the desired system size is reached.

We notice that the superblock Hamiltonian is
constructed using the effective Hamiltonian of
blocks and operator which are renormalized once.
Therefore, the massive truncation because of long
bond is avoided in this algorithm. Bonds in the
superblock are only between new-new sites or new-
old sites. For the construction of a Hamiltonian
matrix of old-new site bond, a new site operator
is highly sparse. However, old sites renormalized
operators are highly dense. The old-old sites in-
teraction in the conventional algorithm generates
a large number of non-zero matrix elements in the
superblock Hamiltonian matrix and the diagonal-
ization of dense matrix goes as m4. But, in the
new algorithm, old-new sites interaction in the su-
perblock generates only a sparse Hamiltonian ma-
trix, and its digonalization scales as m3.

IV. Results

We consider spin-1/2 and spin-1 chains with PBC
of length up to N = 500 to check the accuracy of re-
sults for the Heisenberg Hamiltonian. In this part,
we study the truncation error of density matrix T,
error in relative ground state energy ∆E|E0| and de-

pendence of correlation function C(r) on m. The
correlation function C(r) is defined as

C(r) = ~S0 · ~Sr, (8)

where ~S0 corresponds to the reference spin and ~Sr
is the spin at a distance r from the reference spin.

080006-4



Papers in Physics, vol. 8, art. 080006 (2016) / D. Deyet al.

16 32 64 128 256 512
m

i

10
-12

10
-9

10
-6

10
-3

10
0

T

N = 102

    = 502

100 200 300 400 500
m

i

10
-12

10
-9

10
-6

10
-3

10
0

T

S = 1

S = 1/2

N = 502

Figure 2: Truncation error T of the density ma-
trix for spin-1/2 chain (main). The inset shows
the truncation error for the S = 1 chain. For
S = 1/2, the truncation error follows power law
decay whereas it follows exponential decay for a
S = 1 system.

The relative ground state energy can be defined
as ∆E/|E0|, where ∆E = E(m) − E0 with E0 is
the most accurate value for S = 1 chain [19] and
E0 = E(m = 1200) for S = 1/2 chain.

As discussed earlier, the DMRG is based on the
systematic truncation of the irrelevant degrees of
freedom and the eigenvalues of the density matrix
represent the importance of the respective states.
In the DMRG, only selected states corresponding
to the highest eigenvalues are kept and the rest of
the other states are thrown. We define truncation
error T as

T = 1 −
m∑
i=1

λi (9)

where λi are eigenvalues of density matrix of the
system block arranged in descending order. In the
main Fig. 2, T is shown as function of m for S = 1/2
and the inset shows the same for a S = 1 system
with N=102 and 502. We notice that m ∼ 350 for
S = 1/2 and m ∼ 300 for S = 1 are sufficient to
achieve T = 10−9. In the main Fig. 2, T vs m in
a log-log plot shows a linear behavior, i.e., T for
both system sizes of N = 102 and 502 for S = 1/2
follows a power law T ∝ mαi with α = 4.0 and 3.4,
respectively. The m dependence of T for S = 1
ring is shown in the inset of Fig. 2. The truncation

16 32 64 128 256 512

m

10
-8

10
-6

10
-4

10
-2

10
0

∆
E

/|
E

0
| 0 200 400 600

m

10
-8

10
-6

10
-4

10
-2

10
0

∆
E

/|
E

0
|S=1/2 S=1

Figure 3: Energy accuracy ∆E/|E0| for spin half
chain with PBC (main) which shows a power law
behavior with m. The inset shows the energy ac-
curacy for spin one chain with PBC which shows
exponential behavior with m.

error T in this case decays exponentially, i.e., T ∝
exp(−βmi) with β = 0.03 for both N = 102 and
502.

The relative error in energies ∆E/|E0| for S =
1/2 and 1 with N = 102 spins are shown in Fig. 3,
main and inset respectively. The exact energies
of a spin-1 system is E0/N ∼ 1.4014840386 and
this value is obtained by using conventional DMRG
with m = 2000 and up to N = 100 site chain [19].
Extrapolation of energies with m is done to ob-
tain the above value [19]. We notice that ∆E

E0
for

S = 1/2 systems goes to 10−6 for m = 256 whereas
it goes to 10−8 for m ≈ 500, as shown in the main
Fig. 3. Although similar accuracy of the energy
can be achieved with m = 100 in MPS approach,
the scaling is ∼ m4, rather than ∼ m3 in our algo-
rithm. For S = 1 accuracy of 10−8 can be reached
with m ∼ 450, as shown in the inset of Fig. 3.

The dependence of accuracy of correlation func-
tion |C(r)| of S = 1/2 for N = 130 as a func-
tion of m is shown in the main Fig. 4. We notice
that m = 256 is sufficient to achieve an accuracy
of ∼ 10−4. We also notice that |C(r)| ∝ r−1 with
the well known logarithmic correction specially for
smaller r [26]. We have fitted the correlations with

the well known form |C(r)| = Ar−1 ln1/2(r/r0)
with A = 0.22 and r0 = 0.08 [29]. Deviation in
function for large r is the artifact of finite sys-

080006-5



Papers in Physics, vol. 8, art. 080006 (2016) / D. Deyet al.

1 2 4 8 16 32 64
r

10
-2

10
-1

10
0

|C
(r

)|

m =  64
    = 132
    = 256
    = 400

C(r) = A r
-1

 ln
1/2 

(r/r
0
)

0 0.01 0.02 0.03

N
-1

0

0.01

0.02

0.03

0.04

0.05

C
(r

=
N

/2
)

m = 132
    = 256
    = 400

Figure 4: The main figure shows the variation of
|C(r)| (defined in Eq. (8)) as a function of r for dif-
ferent m. The solid curve in the main figure is the
logarithmic correction formula: Ar−1 ln1/2(r/r0)
with A = 0.22 and r0 = 0.08. [26, 29] The in-
set shows the variation of C(N/2) with the in-
verse of N for different m. The solid curve is the
form: A(N/2)−1 ln1/2(N/2r0) with A = 0.323 and
r0 = 0.08.

tems. In our algorithm, two sites are added sym-
metrically and the new sites are added N/2 sites
apart, consequently the distance between two new
sites is N/2. Therefore, the new-new sites correla-
tion function is C(N/2) and is plotted with N−1

in the inset of Fig 4. We observe that m = 256
is sufficient for N ∼ 200 to achieve sufficient nu-
merical accuracy. The curve behaves almost lin-
early with the logarithmic correction: |C(N/2)| =
2AN−1 ln1/2(N/2r0) with A = 0.323 and r0 = 0.08.

V. Summary

The DMRG is a state-of-the-art numerical tech-
nique for solving the 1D quantum many-body sys-
tems with open boundary condition. However, the
accuracy of the 1D PBC system is rather poor.
The MPS approach gives very accurate results but
the computational cost goes as (o) m5 [23]. Later,
this algorithm was modified and the computational
cost of the modified algorithm goes as (o) p×m3,
where p in general varies linearly with m [28], but
p can go as m2 in case of long range order in the
system. The computational cost of the algorithm
presented in this paper scales as (O) m3 because
of the sparse superblock Hamiltonian and is very

similar to the conventional DMRG. To achieve this
goal, we avoid the multiple times renormalization
of the operators which are used to construct the
superblock. This algorithm can readily be used to
solve general 1D and quasi-1D systems, e.g., J1–J2
model, two-legged ladder. The new algorithm can
be implemented with ease using the conventional
DMRG code.

Our calculation suggests that most of the quan-
tities, e.g., ground state energies, energy gaps and
correlation function, can accurately be calculated
by keeping m ∼ 400. The superfluity stiffness [27]
and dynamical structure factors using the correc-
tion vector technique [24] or continued fraction
method [25] can be calculated with this algorithm.
The symmetries, e.g., spin parity, electron-hole, in-
version, can easily be implemented in this algo-
rithm [24]. This algorithm is used by us in cal-
culating accurate static structure factors and cor-
relation function for J1 − J2 model for a spin-1/2
ring geometry [17].

Acknowledgements - We thank Z. G. Soos
for his valuable comments. M. K. thanks DST
for Ramanujan fellowship and computation facil-
ity provided under the DST project SNB/MK/14-
15/137.

[1] M Tinkham, Introduction to Superconductiv-
ity: Second Edition (Dover Books on Physics)
(Vol I), Dover Publications INC, New York
(2004).

[2] A V Chubukov, Chiral, nematic, and dimer
states in quantum spin chains, Phys. Rev. B
44, 4693 (1991).

[3] L Kecke, T Momoi, A Furusaki, Multimagnon
bound states in the frustrated ferromagnetic
one-dimensional chain, Phys. Rev. B 76,
060407 (2007).

[4] P W Anderson, The resonating valence bond
state in la2cuo4 and superconductivity, Science
235, 1196 (1987).

[5] A Parvej, M Kumar, Degeneracies and ex-
otic phases in an isotropic frustrated spin-1/2
chain, J. Magn. Magn. Mater. 401, 96 (2016).

080006-6



Papers in Physics, vol. 8, art. 080006 (2016) / D. Deyet al.

[6] X L Qi, S C Zhang, Topological insulators and
superconductors, Rev. Mod. Phys. 83, 1057
(2011).

[7] M Nakamura, Tricritical behavior in the ex-
tended hubbard chains, Phys. Rev. B 61, 16377
(2000).

[8] P Sengupta, A W Sandvik, D K Campbell,
Bond-order-wave phase and quantum phase
transitions in the one-dimensional extended
hubbard model, Phys. Rev. B 65, 155113
(2002).

[9] M Kumar, S Ramasesha, Z G Soos, Tun-
ing the bond-order wave phase in the half-filled
extended hubbard model, Phys. Rev. B 79,
035102 (2009).

[10] M Suzuki (Ed.), Quantum Monte Carlo Meth-
ods in Condensed Matter Physics, World Sci-
entific, Singapore (1993).

[11] W Kohn, L J Sham, Self-consistent equa-
tions including exchange and correlation ef-
fects, Phys. Rev. 140, A1133 (1965).

[12] K G Wilson, The renormalization group and
critical phenomena, Rev. Mod. Phys. 55, 583
(1983).

[13] S R White, Density matrix formulation for
quantum renormalization groups, Phys. Rev.
Lett. 69, 2863 (1992).

[14] S R White, Density-matrix algorithms for
quantum renormalization groups, Phys. Rev.
B 48, 10345 (1993).

[15] U Schollwöck, The density-matrix renormal-
ization group, Rev. Mod. Phys. 77, 259 (2005).

[16] K A Hallberg, New trends in density matrix
renormalization, Adv. Phys. 55, 477 (2006).

[17] Z G Soos, A Parvej, M Kumar, Numerical
study of incommensurate and decoupled phases
of spin-1/2 chains with isotropic exchange J1 ,
J2 between first and second neighbors, J. Phys.
Condens. Mat. 28, 175603 (2016).

[18] U Schollwöck, The density-matrix renormal-
ization group in the age of matrix product
states, Ann. Phys. 326, 96 (2011).

[19] P Pippan, S R White, H G Evertz, Effi-
cient matrix-product state method for periodic
boundary conditions, Phys. Rev. B 81, 081103
(2011).

[20] E Dagotto, T M Rice, Surprises on the way
from one- to two-dimensional quantum mag-
nets: The ladder materials, Science 271, 618
(1996).

[21] S Östlund, S Rommer, Thermodynamic limit
of density matrix renormalization, Phys. Rev.
Lett. 75, 3537 (1995).

[22] F Verstraete, D Porras, J I Cirac, Density ma-
trix renormalization group and periodic bound-
ary conditions: A quantum information per-
spective, Phys. Rev. Lett. 93, 227205 (2004).

[23] F Verstraete, J I Cirac, J I Latorre, E Rico,
M M Wolf, Renormalization-group transfor-
mations on quantum states, Phys. Rev. Lett.
94, 140601 (2005).

[24] S K Pati, S Ramasesha, Z Shuai, J L
Brédas, Dynamical nonlinear optical coef-
ficients from the symmetrized density-matrix
renormalization-group method, Phys. Rev. B
59, 14827 (1999).

[25] K A Hallberg, Density-matrix algorithm for
the calculation of dynamical properties of low-
dimensional systems, Phys. Rev. B 52, R9827
(1995).

[26] I Affleck, D Gepner, H J Schulz, T Ziman,
Critical behaviour of spin-s Heisenberg antifer-
romagnetic chains: analytic and numerical re-
sults, J. Phys. A - Math. Gen. 22, 511 (1989).

[27] D Rossini, V Giovannetti, R Fazio, Spin-
supersolid phase in Heisenberg chains: A char-
acterization via matrix product states with pe-
riodic boundary conditions, Phys. Rev. B 83,
140411 (2011).

[28] D Rossini, V Giovannetti, R Fazio, Stiffness in
1D matrix product states with periodic bound-
ary conditions, J. Stat. Mech. 2011, P05021
(2011).

[29] A W Sandvik, Computational Studies of
Quantum Spin Systems, AIP Conf. Proc.
1297, 135 (2010).

080006-7