Neel order in the two-dimensional S=1/2 Heisenberg Model
Néel order in the two-dimensional S = 1
-Heisenberg Model
Ute Löw1
Theoretische Physik, Universität zu Köln, Zülpicher Str.77, 50937 Köln, Germany
(Dated: November 1, 2018)
The existence of Néel order in the S = 1
Heisenberg model on the square lattice at T = 0 is
shown using inequalities set up by Kennedy, Lieb and Shastry in combination with high precision
Quantum Monte Carlo data.
The ground state order of quantum spin systems, in
particular the issue whether the ground state shows long
range magnetic order, has attracted long and continuous
interest. For the prototype of spin models, the antiferro-
magnetic Heisenberg model, the existence of Néel order
at low temperatures was proved in the seminal paper of
Dyson, Lieb and Simon [1] in 1978 for spin S ≥ 1 and
spatial dimension d ≥ 3 and also for S = 1
and d > 3.
Ten years later Kennedy, Lieb and Shastry [2] showed
that also for S = 1
and d = 3 Néel order in the ground
state exists.
The situation in two dimensions is different and more
subtle, since the Mermin-Wagner-Hohenberg theorem
forbids Néel order at finite T , leaving open however the
possibility of Néel order in the ground state. The exis-
tence of Néel order for the two-dimensional model and
S ≥ 1 was shown in [3, 4] and later in [2] by an indepen-
dent derivation of the relevant inequality at T = 0.
However the inequalities sufficient to show Néel order
for S = 1 in the two-dimensional case are not sufficient
to construct an analogous proof for S = 1
. Thus the
case of S = 1
remains an open problem. Still it is pos-
sible to derive inequalities concerning spin-spin correla-
tions at short distances [2] which are violated if Néel order
is present. That is, with a minimum of numerical infor-
mation, the question of Néel order in the ground state
can be decided.
The issue of this paper is to evaluate the spin-spin cor-
relations of the two-dimensional S = 1
antiferromagnetic
Heisenberg model at short distances and demonstrate
that these results combined with the analytic expres-
sions of [2] show the existence of Néel order in the two-
dimensional S = 1
antiferromagnetic Heisenberg model
at T = 0. Such a study has become possible, due to the
developement of high precision Monte-Carlo techniques
over the last decade.
In Ref.[2] Kennedy, Lieb and Shastry used data of
Gross, Sanchez-Velasco and Siggia [5] for a comparison,
however these data clearly deviate from the results pre-
sented here. The authors of [5] used a Quantum Monte
Carlo method without loop updates and with discrete
Trotter time (see below). Their data served only as a
crude comparison to extrapolated Lanczos data and data
produced by the Neumann-Ulam method, which were
the best algorithms to study the properties of the two-
dimensional Heisenberg model in 1988. Today modern
loop algortihms by far outreach both methods.
As will be shown in the following an accurate evalua-
tion of correlation functions at short distances is possible
with modern Quantum Monte Carlo methods, which al-
low us to compute expectation values at very low temper-
atures and even though the short distance results have a
certain finite size and finite temperature correction, these
uncertainties are well controlled and allow to draw defi-
nite conclusions.
The approach and intention of this paper is diffrent
from a completely numerical evaluation of e.g. the corre-
lation length, which involves a calculation of correlations
at long distances and an appropriate extrapolation to in-
finite distances, which cannot be used as a proof of long
range order in any rigorous sense.
At first sight a ”Quantum Monte Carlo algorithm”
seems a puzzling concept, since an important step in
any Monte-Carlo-method is the evaluation of Boltzmann
weights for given energies of the system. For quantum
models these energies are hard if not impossible to calcu-
late. A key idea to make Monte Carlo methods applicable
to quantum systems is to map the quantum model onto a
classical model by introducing an extra dimension, usu-
ally referred to as Trotter-time [6].
In the first generation of algorithms this mapping
was straightforwardly applied to the quantum Heisen-
berg model. Though this allowed for a wealth of new
studies of the finite temperature properties in one and in
particular in two-dimensional systems, these algorithms
had two major drawbacks, which became most evident
at low temperatures. Firstly the extra Trotter dimension
was discretized, introducing the number of time slices as
a parameter which had to be eliminated from the final
results by an extrapolation. Secondly the update pro-
cedure, i.e. the construction of new independent config-
urations, was done locally. As a consequence one had
to move through the lattice site by site several times to
obtain a configuration independent of the starting config-
uration and useful for a new evaluation of an observable.
A first improvement was introduced by the so called
loop-algorithms [7], which uses nonlocal updates similar
to the Swendsen-Wang algorithm for classical models. A
second and important step towards high precision Quan-
tum Monte Carlo techniques were algorithms which work
directly in the Euclidian time continuum [8] and require
no extrapolation in Trotter time. For the algorithm [9]
used for the analysis presented here no approximations
enter, and statistical errors are the only source of inac-
curacy.
Since this work intends to produce highly accurate
http://arxiv.org/abs/0704.0713v1
data it seems appropriate to assess the precision of the
method by a comparison with exact results. The best
candidate for such a comparison are the correlations of
one-dimensional systems evaluated by the Bethe-ansatz
with almost arbitrary precision up to distance seven [10]
and with results for finite chains from Ref.[11]. This is
done in the Appendix for chains of 400 sites at T=0.005.
After these introductory remarks we now return to our
actual goal, which is the two-dimensional system. Our
starting point is a S = 1
Heisenberg model
x,yεΛ
~Sx~Sy (1)
with nearest neighbour interaction on a finite square lat-
tice Λ with an even number of sites in every direction
and periodic boundary conditions.
The Fourier transform of the spin-spin correlation
function at T = 0 is given by
gq = 〈S−qSq〉 =
e−iqx〈S30S
x〉 (2)
where
e−iqxS3x. (3)
For the corresponding finite temperature expectation
value of gq an upper bound fq was derived in [1]. The
T = 0 limit of this bound was obtained in Ref.[4] and
a direct proof of the bound at T = 0 was given in [2].
Following the notation and arguments of Kennedy, Lieb
and Shastry [2] the inequality for d = 2 reads
gq ≤ fq for q 6= Q (4)
where fq =
12Eq−Q
, Eq = 2 − cos q1 − cos q2, Q =
(π, π) and −e0 is the ground state energy per site of the
Heisenberg model Eq.1 on the lattice Λ.
The fundamental idea is, that the existence of Néel
order in the limit of infinite system size corresponds to a
delta-function in the Fourier transform of the spin-spin
correlation gq at Q. This means, if Eq. 4 is integrated
over the whole Brillouin zone one finds in the case of Néel
order
d2q fq ≥
d2q gq = S(S + 1)/3 (5)
where m2 is the coefficient of the delta-function at Q.
If there is no Néel order m2 is zero. By numerically
evaluating the integral over fq, and by using exact varia-
tional upper and lower bounds on the ground state energy
−e0 one sees, that the above inequality and its analogon
for d ≥ 3 cannot be fulfilled with m2 = 0 and S ≥ 1,
which proves Néel order.
Inequalities of type Eq. 5 are not sufficient to prove
the existence of a nonzero m2 for d = 2, 3 and S = 1
but a new relation is obtained by multiplying gq by cos qi
and again integrating over the Brillouin zone:
ddq gq cos qi = 〈S
〉 = −e0/3d (6)
with i=1,2 for d=2 and i=1,2,3 for d=3 and δi the unit
vector in i-direction and the value of the ground state
energy form Ref.[12] is e0 = 0.669437(5).
Carrying out an analogous integral over fq and using
again Eq.4 one finds:
d2Eq−Q
cos qi
were the f+ means the positive part of a function, which
equals f, when f is positive and is zero otherwise.
Again Eq.7, which is valid if no Néel order exists, was
shown to be violated for d = 3 and S = 1
in Ref.[2] by
using bounds on e0 and thus the existence of Néel order
was proved also for d = 3 and S = 1
For S = 1
and d = 2 one cannot construct a contra-
diction by using only the ground state energy. Here more
input from numerical data is needed. This can be incor-
porated by multiplying gq by cos(mqi) with m = 2, 3...
and again integrating over the whole Brillouin zone:
d2q gq cos(mqi) = 〈S
〉 (8)
with i=1,2.
Next, defining ḡ(n) as
ḡ(n) =
(−1)m〈S30S
〉 (9)
and using again inequality 4 one constructs the follow-
ing relations involving the correlation functions:
ḡ(n) =
2n+ 2
(−1)m{cos(mq1) + cos(mq2)} gq
2n+ 2
(−1)m{cos(mq1) + cos(mq2)}+ fq.
Whenever the inequality Eq. 10 is violated for a cer-
tain n, a nonzero m2 multiplying the delta-function at
Q is needed and therefore the existence of Néel order is
proved.
The ḡ(n) as defined in Eq.9 were calculated by the
QuantumMonte Carlo method [9]. The results, displayed
in table I, show that the ḡ(n) calculated by Quantum
Monte Carlo cross the bound obtained by integrating
over fq at n = 8. This is also depicted in Fig. 1. Thus
inequality Eq.10 is violated and Néel order must exists in
n Bound T = 0.005 T = 0.025 T = 0.075
1 2.297e-01 1.80799e-01 ± 3.63e-06 1.80794e-01 1.80792e-01
2 1.714e-01 1.40308e-01 ± 5.63e-06 1.40302e-01 1.40298e-01
3 1.383e-01 1.17686e-01 ± 6.84e-06 1.17678e-01 1.17670e-01
4 1.166e-01 1.03005e-01 ± 7.67e-06 1.02997e-01 1.02985e-01
5 1.013e-01 9.27815e-02 ± 8.27e-06 9.27743e-02 9.27544e-02
6 8.990e-02 8.52115e-02 ± 8.73e-06 8.52048e-02 8.51770e-02
7 8.107e-02 7.93875e-02 ± 9.10e-06 7.93811e-02 7.93436e-02
8 7.400e-02 7.47551e-02 ± 9.40e-06 7.47496e-02 7.47012e-02
9 6.820e-02 7.09844e-02 ± 9.64e-06 7.09795e-02 7.09191e-02
10 6.334e-02 6.78504e-02 ± 9.85e-06 6.78464e-02 6.77734e-02
11 5.921e-02 6.52055e-02 ± 1.00e-05 6.52021e-02 6.51163e-02
12 5.563e-02 6.29418e-02 ± 1.02e-05 6.29389e-02 6.28404e-02
13 5.252e-02 6.09835e-02 ± 1.03e-05 6.09806e-02 6.08695e-02
14 4.976e-02 5.92718e-02 ± 1.04e-05 5.92692e-02 5.91456e-02
15 4.732e-02 5.77638e-02 ± 1.06e-05 5.77617e-02 5.76255e-02
TABLE I: Bound obtained by integrating numerically over
the right hand side of Eq.10 compared with ḡ(n) for a 40×40
lattice and different temperatures.
the two-dimensional antiferromagnetic Heisenberg model
with S = 1
at T = 0.
There are three type of corrections to the data of table
I, which need to be taken into account, but which, as
we shall show in the following, do not change the above
conclusion of a crossing of the curves at n = 8:
(i) effects of finite temperature,
(ii) effects of the finiteness of the system,
(iii) statistical errors.
In the following we comment on how these corrections
modify the data.
(i) The Quantum Monte Carlo data presented are at
T ≥ 0.005. The overall effect of finite temperature is to
lower the absolute value of the correlations and therefore
also the value of the ḡ(n). The effect of finite temperature
is to shift the crossing of the bound and ḡ(n) to larger n,
or eventually to destroy a crossing completely.
The functional dependence of the internal energy
U(T ), which up to an overall factor 3z (z = 2 is the
coordination number of the two-dimensional square lat-
tice) equals the correlation-function at distance one, has
been determined for low T by spin wave theory [13, 14]
U(T ) = −e0 + bT
3. (11)
The coefficient is given in [15] as b =
≈ 0.2853626,
so the correction for distance one is ≈ b
10−7, which is
two orders of magnitude smaller than the statistical error,
(see point (iii)).
For distances larger than one, we fitted the data as
a function of temperature (taking the exponent of T as
n Bound T = 0.025 T=0.025 extrapolated
1 2.297e-01 1.80794e-01 1.80791e-01 ± 5.09e-06
2 1.714e-01 1.40302e-01 1.40295e-01 ± 7.87e-06
3 1.383e-01 1.17678e-01 1.17668e-01 ± 9.53e-06
4 1.166e-01 1.02997e-01 1.02983e-01 ± 1.07e-05
5 1.013e-01 9.27743e-02 9.27534e-02 ± 1.15e-05
6 8.990e-02 8.52048e-02 8.51762e-02 ± 1.21e-05
7 8.107e-02 7.93811e-02 7.93428e-02 ± 1.26e-05
8 7.400e-02 7.47496e-02 7.46995e-02 ± 1.30e-05
9 6.820e-02 7.09795e-02 7.09154e-02 ± 1.34e-05
10 6.334e-02 6.78464e-02 6.77663e-02 ± 1.37e-05
11 5.921e-02 6.52021e-02 6.51035e-02 ± 1.39e-05
12 5.563e-02 6.29389e-02 6.28188e-02 ± 1.41e-05
13 5.252e-02 6.09806e-02 6.08346e-02 ± 1.43e-05
14 4.976e-02 5.92692e-02 5.90923e-02 ± 1.45e-05
15 4.732e-02 5.77617e-02 5.75473e-02 ± 1.46e-05
TABLE II: Bound obtained by integrating numerically over
the right hand side of Eq.10 compared with ḡ(n) extrapolated
for N=40,36,32,24 at T = 0.025.
fit parameter) for T = 0.005, 0.025, 0.05, 0.075 and found
the corrections due to finite temperature all of the order
of 10−5, which is the order of the statistical error. There-
fore we do not give any finite temperature corrections.
(ii) The absolute value of the correlations in the ther-
modynamic limit are smaller than in systems of finite
size. This means that the effect of finite system size is
opposite to the effect of temperature. The finite size be-
haviour of the ground state energy is well studied for
the Heisenberg model on the square lattice. Arguments
originating from the quantum nonlinear sigma model de-
scription [16] of the Heisenberg model to lowest order in
system size give
− e0 = −e0(N) +
, with c > 0 (12)
where −e0(N) is the ground state energy of a system of
size N ×N . Though the corrections are not substantial,
they do effect the results, and taking into account, that
the finite size errors in contrast to the finite temperature
effects, might falsely lead to a crossing, we extrapolated
the data for N = 24...40 using the functional dependence
Eq.12, which we found well satisfied also for larger dis-
tances. The results are shown in table II. One sees that
the numeric values are changed but the crossing point is
still at n = 8.
(iii) We compute ∆x = 1√
〈x2〉 − 〈x〉2 (where the
observable x stands for the value of a correlation at a
given distance, temperature and system size and NMC is
the number of Monte Carlo iterations), which is a reliable
estimate for the statistical error of the mean value 〈x〉,
since for the algorithm of Ref.[9] the autocorrelation time
is of order one and the Monte Carlo configurations are al-
most independent. To assess the quality of our error anal-
0 5 10 15 20
bound 
 24x24  
 36x36 
 40x40 
FIG. 1: Bound on ḡ(n) obtained from Eq.9 and ḡ(n) for 24×
24, 36× 36 and 40× 40 at T = 0.025.
Distance Quantum Monte Carlo Bethe-Ansatz
0 0.25000000 ( 0)
1 -0.14771586 (198) -0.1477157268
2 0.06067787 (324) 0.0606797699
3 -0.05024194 (282) -0.0502486272
4 0.03464515 (281) 0.0346527769
5 -0.03088096 (260) -0.0308903666
6 0.02443619 (255) 0.0244467383
7 -0.02248413 (242) -0.0224982227
8 0.01895736 (236)
TABLE III: Correlations for a chain with N = 400 sites at
T/J=0.005 compared with results from Ref. [10].
ysis we also returned to the case of the one-dimensional
antiferromagnetic Heisenberg model ( see Appendix ) and
compared results with independent streams of random
numbers.
To calculate an upper limit to the errors of ḡ(n), the er-
rors of the correlations where added up ( being evaluated
with the same configurations, they are not independent).
To conclude, the error analysis shows that the short
range correlations entering Eq.9 were determined with
sufficiently high accuracy to prove the existence of a
crossing of the bound and the Quantum Monte Carlo
data for ḡ(n) at n = 8 and therefore to show the
existence of long range order.
Appendix
(1) In this Appendix we list the correlations of a one-
dimensional Heisenberg model with periodic boundary
conditions and chain length N = 400 at T=0.005 com-
pared with results of Ref.[10] for infinite chain length and
T = 0.
For the internal energy U(T ) of the Heisenberg chain
the temperature dependence for low T is U(T ) = −e10 +
aT 2 with the ground-state energy e10 = 0.4431471804
for 400 sites and e10 = −
+ ln 2 for the infinite size
system[17]. and the coefficient a = 1
given in Ref.
[18, 19]. This means that the correction for the corre-
lations in tableIII due to finite temperatures are of the
order of 10−5.
(2) The exact values of the correlation functions
[11, 20] for distance one and two at T = 0 for a chain
with 400 sites are 〈S30S
1〉400 = −0.147717441765735 and
〈S30S
2〉400 = 0.0606813790491800. The above data show
that the error analysis concerning statistical errors and
finite temperature effects is consistent.
Acknowledgement
I am indebted to Prof. E.H. Lieb for bringing the prob-
lem of longrange order to my attention and for his interest
in this work.
[1] F.J. Dyson, E.H. Lieb and B. Simon, J.Stat.Phys. 18
335-383 (1978).
[2] T. Kennedy, E. H.Lieb and S. Shastry, J.Stat.Phys. 53,
1019-1030,(1988).
[3] I. Affleck, T. Kennedy, E.H. Lieb and H. Tasaki, Comm.
Math. Phys. 115:477-528.
[4] E. Jordão Neves and J. Fernando Perez, Phys.Lett.114A
331-333 (1986).
[5] M. Gross, E. Sanchez-Velasco, and E. Siggia, Phys.Rev.B
39 2484(1989).
[6] M. Suzuki, Commun. Math. Phys. 51, (1976).
[7] H.G. Evertz, G. Lana, M.Marcu, Phys. Rev. Lett 70, 875
(1993).
[8] E. Farhi and S. Gutmann, Ann.Phys. (N.Y.)213, 182
(1992).
[9] B. B. Beard and U. -J. Wiese, Phys. Rev. Lett. 77, 5130
(1996).
[10] J. Sato, M. Shiroishi, M. Takahashi hep-th/0507290.
[11] J.Damerau, F.Göhmann, N.P.Hasenclever, A.Klümper,
cond-mat/0701463.
[12] A. W. Sandvik, Phys.Rev.B56 (1997) 11678.
[13] R. Kubo, Phys.Rev.87, 568 (1952).
[14] T. Oguchi, Phys. Rev.117, 117 (1960).
[15] M. Takahashi, Phys. Rev.B 40, 2494 (1989).
[16] S. Chakravaty, B.I. Halperin, D.R. Nelson,
Phys.Rev.B39,2344(1989).
[17] L. Hulthén, Arkiv Mat.Astron.Fysik 26A,1 (1938).
[18] H.M. Babujian, Nucl.Phys. B215, 317 (1982).
[19] I. Affleck, Phys.Rev.Lett. 56,746 (1986).
[20] J. Damerau, private communication.
