Skip to main content

Full text of "Classical and quantum two-dimensional anisotropic Heisenberg antiferromagnets"

See other formats

Classical and quantum two-dimensional anisotropic Heisenberg antiferromagnets 





M. Holtschneider, 1 S. Wessel, 2 and W. Selke 1 

1 Institut fur Theoretische Physik, RWTH Aachen, 52056 Aachen, Germany 
2 Institut fur Theoretische Physik III, Universitat Stuttgart, 70550 Stuttgart, Germany 

The classical and the quantum, spin S = 4, versions of the uniaxially anisotropic Heisenberg 
antiferromagnet on a square lattice in a field parallel to the easy axis are studied using Monte Carlo 
techniques. For the classical version, attention is drawn to biconical structures and fluctuations 
at low temperatures in the transition region between the antiferromagnetic and spin-flop phases. 
For the quantum version, the previously proposed scenario of a first-order transition between the 
antiferromagnetic and spin-flop phases with a critical endpoint and a tricritical point is scrutinized. 

PACS numbers: 75.10.Hk, 75.10.Jm, 75.40.Mg, 05.10.Ln 















Uniaxially anisotropic Heisenberg antiferromagnets in 
an external field along the easy axis have attracted much 
interest, both theoretically and experimentally, due to 
their interesting structural and critical properties. In 
particular, they display a spin-flop phase, and multi- 
critical behavior occurs at the triple point of the an- 
tiferromagnetic (AF), spin- flop (SF) and paramagnetic 

A prototypical model for such antiferromagnets is the 
XXZ model, with the Hamiltonian 

H = J £ [ A(5f Sf + SfS*) + SfS?] - H S? , 


where the sum runs over neighboring spins of a cubic, 
dimension d = 3, or square lattice, d = 2. The coupling 
constant J and the field H are positive; the anisotropy 
parameter A may range from zero to one. Furthermore, 
Sf, Sf , and Sf denote the spin components at lattice 
site i. 

For the three-dimensional case, early renormalization 
group arguments^ and Monte Carlo simulations^ sug- 
gested that the triple point is a bicritical point with 
0(3) symmetry. Only a few years ago, this scenario has 
been questioned, based on high-order perturbative renor- 
malization group calculationSfii It has been predicted 
that there may be either a first order transition, or that 
the 'tetracritical biconical' - fixed point, due to an inter- 
vening 'mixed' or 'biconical' phase in between the AF 
and SF phase s 12 ' 13 ' 14 , may be stable. 

In two dimensions, conflicting predictions on the 
nature of the triple point have been put forward 
recentl y 15 ' 16 i 17 i 18 , when analyzing the classical version of 
the above model with spin vectors of unit length, and the 
quantum version with spin S = h ■ 

Indeed, in the classical case, simulational evidence for a 
narrow (disordered) phase between the AF and SF phases 
has been presented^ 5 -, extending presumably down to zero 
temperature^ On the other hand, in the quantum case, 
based on simulations as well, a direct transition of first 
order between the AF and SF phases has been argued to 
occur at low temperatures ! 18 ! 19 ' 20 

FIG. 1: Ground state configurations of the classical model 
sketched by the directions of spins on the two sublattices (i. e. 
at neighboring sites), from left to right: AF, SF, and biconical 
state. The circles denote the trivial degeneracy in the xy- 

Obviously, experimental data have to be viewed with 
care because deviations from the XXZ Hamiltonian, 
Eq. |l| , such as crystal field anisotropics or longer-range 
interactions, may affect relevantly the critical behavior 
of the triple point ^^14,21,22 

In the following, we present results from large-scale 
Monte Carlo simulations of the XXZ model on a square 
lattice for both the classical and the quantum variant. 
In the quantum Monte Carlo simulations, the method of 
the stochastic series expansion (SSE}22 is used, and the 
standard Metropolis algorithm is applied for the classical 
case. The simulations are augmented by a ground-state 
analysis of the classical model, showing the significance 
of biconical structures. The outline of the paper is as 
follows: First we shall discuss our findings on the classical 
model, followed by a section on the quantum version of 
the XXZ model. A summary concludes the paper. 


The ground states of the classical model on a square 
lattice, see Hamiltonian ([T]), can be determined exactly. 
The AF structure is stable for magnetic fields below the 
critical value 

H cl = 4JV1 - A 2 , 


while for larger fields the SF state is energetically favor- 
able. At H c t, the tilt angle &sf of the SF structures, see 








J para- 
t magnetic 



0.24 0.28 0.32 0.36 0.40 
k B T/J 

FIG. 2: Detail of the phase diagram of the XXZ model on a 
square lattice with A = | , see Ref . [3 Squares refer to the 
boundary of the SF, circles to that of the AF phase. The solid 
line refers to the magnetic field H/ J — 2.41, where the prob- 
ability distribution P(6 m ,8 n ), depicted in Fig. [3] has been 
obtained. Here and in the following figures error bars are 
shown only if the errors are larger than the symbol size and 
dotted lines are guides to the eye. 

FIG. 3: Joint probability distribution P(8 m ,8 n ) showing the 
correlations between the tilt angles 8 m and 8 n on neighboring 
sites m and n for a system of size L — 80 at H/J = 2.41, 
k B T/J = 0.255, and A = §. P(8 m ,8 n ) is proportional to 
the gray scale. The superimposed black line depicts the rela- 
tion between the two angles in the biconical ground state, see 
Eq. gjl. 

Fig. [JJ is given by 

9 SF = arccos J . (3) 

Increasing the field beyond H C 2 = 4J(1 + A), all spins 
perfectly align in the z-direction. 

At the critical field H&, see Eq. the ground state 
is degenerate in the AF, the SF, and biconical structures, 
as illustrated in Fig. [TJ This degeneracy in the biconi- 
cal configurations, following from straightforward energy 
considerations, seems to have been overlooked in the pre- 
vious analyses. The structures may be described by the 
tilt angles, 0\ and 82, formed between the directions of 
the spins on the two sublattices of the antiferromagnet 
and the easy axis. For a given value of 9%, the other 
angle 62 is fixed by 

/ VT -- A~2 - COS0! \ 

6 2 = arccos . . (4) 

\1 - VT^A^cosfli/ V ' 

Obviously, the biconical configurations transform the AF 
into the SF state: The spins on the "up-sublattice" of the 
AF structure, with the spins pointing into the direction of 
the field, may be thought of to vary from 9\ = to ±#sf, 
while the spins on the " down-sublattice" vary simultane- 
ously from 9% = 7t to T^SF- Accordingly, 9\ determines 

uniquely 62 and vice versa. Apart from this continu- 
ous degeneracy in 9\ (or 92 ), there is an additional rota- 
tional degeneracy of the biconical configurations in the 
spin components perpendicular to the easy axis, the xy- 
components, as for the SF structure, see Fig. UJ. These 
components are, of course, antiferromagnetically aligned 
at neighboring sites. 

To study the possible thermal relevance of the bi- 
conical structures at T > 0, we performed Monte 
Carlo simulations analyzing the joint probability distri- 
bution P(9 m ,9 n ) for having tilt angles 9 m and 9 n at 
neighboring sites, m and n. For comparison with the 
previous studies ^ 15 ' 16 we set A = |, leading to the 
phase diagram depicted in Fig. [21 For example, fixing 
the field at H = 2.41J, we observed at k B T/J m 0.255 
an Ising-type transition on approach from higher tem- 
peratures and a Kosterlitz-Thouless-type transition on 
approach from the low-temperature side, extending our 
corresponding previous findings^ to even lower temper- 
atures, and in agreement with recent results^ Indeed, 
as depicted in Fig. [31 in that part of the phase diagram, 
being in the vicinity of the very narrow intervening, sup- 
posedly disordered phase, the joint probability P(9 m , 9 n ) 
exhibits a line of local maxima following closely Eq. (J4|, 
obtained for the ground state. That behavior is largely 
independent of the size of the lattices we studied. Simi- 
lar signatures of the biconical structures are observed in 


1.32 - 




I I 


h-BH . 








tx. ; 


, , , i , , , i , 

, i i i i i 

0.08 0.10 0.12 0.14 0.16 

k B T/J 

FIG. 4: Phase diagram of the XXZ Heisenberg antiferromag- 
net with spin-i and A = |. The straight solid lines denote 
the choices of parameters where our very extensive simula- 
tions, discussed in the text, have been performed. The arrows 
mark the previously^ suggested locations of the tricritical 
point (Tt) and the critical endpoint (T ce ). 

the simulations at nearby temperatures, when fixing the 
field at H = 2.41J, as well as in the vicinity of the en- 
tire transition region between the AF and SF phases, see 
Fig. [2j at higher fields and temperatures. 

Accordingly, we tend to conclude that biconical fluctu- 
ations are dominating in the narrow intervening phase. 
Whether that phase exists as a disordered phase down 
to the ground state or whether there is a stable biconical 
phase in two dimensions, remain open questions, being 
beyond the scope of this article. 

Note that our additional Monte Carlo simulations for 
the anisotropic XY antiferromagnet in a field on a square 
lattice show that the analogues of 'biconical' structures 
(the orientation of the spins being now given by the two 
tilt angles only) and fluctuations play an important role 
near the transition regime between the AF and SF phases 
in that case as well. In fact, Eq. Q provides an excellent 
guidance for interpreting our simulational data similar to 
the ones presented in Fig. 3. 


The aim of the study on the quantum version, S = |, 
of the XXZ model, Eq. ([I]) , has been to check the previ- 
ously suggested scenario of a first-order phase transition 
between the AF and SF phases extending up to a critical 
endpoint and with a tricritical point on the AF phase 

boundary, see Fig. [4] 

We performed quantum Monte Carlo (QMC) simula- 
tions in the framework of the stochastic series expan- 
sion (SSE) 23 using directed loop updates 2 ^. We con- 
sider square lattices of L x L sites with the linear di- 
mension L ranging from 2 to 150, employing full periodic 
boundary conditions. Defining, as usual^ a single QMC 
step as one diagonal update followed by the construction 
of several operator-loops, each individual run typically 
consists of 10 6 steps and is preceded by at least 2 • 10 5 
steps for thermal equilibration. Averages and error bars 
are obtained by taking into account results of several, 
ranging from 8 to 32, Monte Carlo runs, choosing dif- 
ferent initial configurations and random numbers. Es- 
pecially for large systems and low temperatures we ad- 
ditionally utilize the technique of parallel tempering (or 
exchange Monte Carlo) 25 ' 26 to enable the simulated sys- 
tems to overcome the large energy barriers between con- 
figurations related to different phases more frequently. 
We typically work with a chain of 16 to 32 configura- 
tions in parallel which are simulated at different equally 
spaced temperature or magnetic field values allowing for 
an exchange of neighboring configurations after a con- 
stant number of QMC steps. The achieved reduction of 
the autocorrelation times, e.g. of the different magneti- 
zations discussed below, amounts up to several orders of 
magnitude and therefore results in significantly smaller 
correlations between subsequent measurements which, in 
turn, allows for shorter simulation times. 

To determine the phase diagram and to check against 
previous worki&, we calculated various physical quanti- 
ties including the z-component of the total magnetiza- 

M z 

72 X^ 5 * 


and the square of the ^-component of the staggered mag- 

E< 5 - 


summing over all sites, i a and z&, of the two sublattices 
of the antiferromagnet. A useful quantity in studying 
the SF phase is the spin-stiffness p s which is related to 
the change of the free-energy on imposing an infinitesi- 
mal twist on all bonds in one direction of the lattice. In 
QMC simulations the spin-stiffness can conveniently be 
measured by the fluctuations of the winding numbers W x 

The winding numbers themselves are given by 


W a = - {N+ - N~) 








n1> "3.0 


FIG. 5: Positions of the maxima of the magnetization his- 
tograms as a function of the inverse system size. The inset 
exemplifies two histograms for systems of size L — 32 (circles) 
and L = 150 (squares) at ksT / J = 0.13 and the coexistence 
fields H/J = 1.23075 and H/J = 1.232245. 

FIG. 6: Doubly logarithmic plot of the staggered magnetiza- 
tion (M^) 2 vs. the system size L at H/J = 1.225 for the 
temperatures k B T / J = 0.095, 0.0955, 0.096, 0.0965, 0.097, 
and 0.0975 (from bottom to top). The straight line propor- 
tional to Li illustrates the expected finite-size behavior close 
to a continuous transition of Ising type. 

where and N~ denote the number of opera- 

tors S^~S~ and S~ Sj~ in the SSE operator sequence with 
a bond in the ct-direction, a € {x,y}. 

All data for the quantum model presented here are 
obtained at an anisotropy parameter of A = | to allow 
for comparison with previous finding s 15 ' 18 . The phase 
diagram in the region of interest, where all three phases, 
the AF, the SF, and the paramagnetic phase occur, is 
displayed in Fig. [4j 

The earlier studj*^ asserted a phase diagram with a tri- 
critical point at ksTt/J ~ 0.141 and a direct first-order 
transition between the SF and AF phases below the criti- 
cal endpoint at ksT cc / J w 0.118, see Fig.0J In detail the 
authors identified a first-order AF to paramagnetic tran- 
sition at ksT/J = 0.13 by means of an analysis of the 
magnetization histograms p(M z ). We studied that case, 
improving the statistics and considering even larger lat- 
tice sizes. Indeed, as expected for a discontinuous change 
of the magnetization, the histograms of finite systems are 
confirmed to display two distinct maxima corresponding 
to the ordered and the disordered phase in the vicinity of 
the AF phase boundary (see inset of Fig. [5]) . Note how- 
ever, that such a two-peak structure can also be found 
for small systems at a continuous transition, with a sin- 
gle peak in the thermodynamic limit. Thence, a careful 
finite-size analysis is needed to, possibly, discriminate the 
two different scenarios. We simulated lattice sizes with 
up to 150 x 150 spins adjusting the magnetic field such 
that coexistence of the phases, i.e. equal weight of the 

two peaks, is provided. As depicted in Fig. [5] the posi- 
tions of the maxima as a function of the inverse system 
size exhibit a curvature, which becomes more pronounced 
for larger lattices. In contrast, in the previous analysis^ 
at the same temperature, linear dependences of the peak 
positions as a function of 1/L had been presumed, lead- 
ing to distinct two peaks in the thermodynamic limit. We 
conclude, that the previous claims of a first-order transi- 
tion at ksT I J = 0.13 and of the existence of a tricritical 
point at ksT t / J ss 0.141 needs to be viewed with care. 
Indeed, the tricritical point seems, if it exists at all, to 
be shifted towards lower temperatures. 

In the previous worki£ a direct transition of first order 
between the AF and SF phases has been suggested to take 
place at lower temperatures, fc^T/J < kBT ce /J « 0.118. 
To check this suggestion we studied the system at con- 
stant field H/J = 1.225, where such a direct tran- 
sition would occur, see Fig. |U Calculating the ex- 
pectation values of the different magnetizations as well 
as the corresponding histograms we obtain an esti- 
mate of the critical temperature of the AF phase, 
k B T AF = 0.09625 ± 0.0005. 

Surprisingly, approaching the transition from the AF 
phase, the finite-size behavior of the squared staggered 
magnetization (M^) 2 , being the AF order parameter, is 
still consistent with a continuous transition in the Ising 
universality class: As illustrated in Fig.[6]the asymptotic 
region is very narrow, similar to the observations in the 


classical modeL 15 ' 16 The dependence on the system size 
seems to obey (M s z t ) 2 oc L 1 / 4 right at the transition, as 
expected for the Ising universality class.— 

Furthermore, approaching the transition from the SF 
phase, an analysis of the spin-stiffness p s at the same field 
value of H / J = 1.225 results in about the same transition 
temperature, k B T SF /J = 0.09625 ± 0.001. Thence, there 
may be either a unique transition between the SF and AF 
phases, or, as observed in the classical case, an extremely 
narrow intervening phase, with phase boundaries of Ising 
and Kosterlitz-Thouless (KT) type. 

To determine, whether a KT transition describes the 
disordering of the SF phase, we check the theoreti- 
cal predictio n 28 ! 29 that for the infinite system the spin- 
stiffness is finite within the SF phase, takes on a universal 
value at the KT transition related to Tkt by 


' KT 


fc B Ti 



and discontinuously vanishes in the disordered phase. As 
depicted in Fig. [3 the spin-stiffness p s at T = Tgp seems 
to be, at first sight, significantly larger than the KT- 
critical value given by Eq. ©. Indeed, in the earlier 
studyi^ it has been argued, based on similar observations, 
that there is a direct first order AF to SF transition. 
However, the finite-size effects close to the transition de- 
serve a careful analysis: For the KT scenario, renormal- 
ization group calculation a 30 ' 31 predict the asymptotic size 
dependence at T = Tkt to obey 

p s {T = T KTl L) = 
p s (T = T KT ,L = oo) (l + 


21nL - Co 


where Co denotes an apriorily unknown, non-universal, 
parameter. By studying the quantity^ 2 - 

C{L) = -2 

k R T 

- 2 

- In I 


which, according to Eqs. [9land[T0l converges for L — ► oo 
and T = Tkt to the value Co at a KT transition, we 
obtain a rough estimate of Co ~ 5. A prediction of the 
finite-size behavior at Tsf is obtained by inserting this 
value, Co = 5, into Eq. (Tl0|) . Comparing the data of the 
spin-stiffness p s in the direct vicinity of the boundary of 
the SF phase with the prediction according to the KT 
theory, see Fig. [7] b) , one may conclude that the lattice 
sizes accessible by simulations, L < 64, seem to be too 
small to capture the asymptotic finite-size behavior. In 
any event, in case of a KT transition, the spin-stiffness p s 
drops asymptotically very rapidly to its universal criti- 
cal value as a function of system size, being consistent 
with the relatively large values for the simulated finite 
lattices. Thus, a scenario with a KT transition between 
the SF and a narrow intervening disordered phase cannot 
be ruled out by the present large-scale simulations down 
to temperatures as low as ksT/J = 0.09625 ± 0.001. 

3 ) 0.10 







k B T/J 






FIG. 7: a) Spin stiffness p s /J vs. temperature ksT/J for the 
different system sizes L — 8 (circles), 16 (squares), 32 (di- 
amonds), and 64 (triangles). The straight line denotes the 
critical value of the spin-stiffness according to the formula of 
Nelson and Kosterlita^, see Eq. ©. 

b) Finite-size behavior of the spin-stiffness p s /J at H/ J — 
1.225 as a function of the inverse system size, 1/L for 
the temperatures k B T/J = 0.0955, 0.096, 0.0965, 0.097, 
and 0.0975 (from top to bottom). The dashed curve illustrates 
the estimated asymptotic behavior according to Eq. (|10p 
with fcsTkT / J = 0.09625 and Co = 5, the corresponding crit- 
ical value p s (Tkt,L — oo) w 0.0613 is marked by the filled 

Of course, it is desirable to quantify the role of bicon- 
ical fluctuations in the quantum case as well. However, 
accessing the probability distributions of the tilt angles 
studied in Sect. |TT] for the quantum case is beyond the 
scope of the present numerical analysis. 


We studied the classical and quantum, S = |, versions 
of the XXZ Heisenberg antiferromagnet on the square lat- 
tice in an external field along the easy axis. The model 


is known to display ordered AF and SF as well as dis- 
ordered, paramagnetic phases. Here we focused atten- 
tion to the region of the phase diagram near and below 
the temperature where the two boundary lines between 
the AF and the SF phases and the disordered phase ap- 
proach each other, meeting eventually at a triple point. 
We performed Monte Carlo simulations, augmented, in 
the classical case, by a ground state analysis. 

In the classical version, we presented first direct ev- 
idence for the importance of biconical structures in the 
XXZ model. Indeed, such configurations do exist already 
as ground states at the critical field H c \, separating the 
AF and SF phases. The interdependence of the two tilt 
angles, characterizing the biconical ground states, per- 
sists at finite temperatures, in the region where the nar- 
row phase between the AF and SF phases is expected to 
occur. Indeed, the joint probability distribution of the 
tilt angles at neighboring sites demonstrates the thermal 
significance of those configurations. Previous arguments 
on 0(3) symmetry in that region and down to zero tem- 
perature thus have to be viewed with care. The results 
of the present simulations suggest that, if the biconical 
configurations do not lead to a stable biconical phase in 
two dimensions, the narrow intervening phase is a disor- 
dered phase characterized by biconical fluctuations. In 
this sense the "hidden bicritical point" at T = may 
then be coined into a "hidden tetracritical point." 

In the quantum version, previous simulations sug- 

gested, on lowering the temperature, the existence of a 
tricritical point on the boundary line between the AF and 
disordered phases, followed by a critical endpoint being 
the triple point of the AF, SF and disordered phases, and 
eventually by a first-order transition between the AF and 
SF phases at sufficiently low temperatures. The present 
simulations, considering larger system sizes and improved 
statistics, provide evidence that this scenario, if it exists 
at all, has to be shifted to lower temperatures than pro- 
posed before. Of course, simulations on even larger lat- 
tices and lower temperatures would be desirable, but are 
extremely time consuming. 

A clue on possible distinct phase diagrams for the clas- 
sical and quantum versions may be obtained from an 
analysis of biconical fluctuations in the quantum case. 
Experimental studies on signatures of those fluctuations 
are also encouraged. 


Financial support by the Deutsche Forschungsgemein- 
schaft under grant No. SE 324/4 is gratefully acknowl- 
edged. We thank A. Honecker, B. Kastening, R. Leidl, 
A. Pelissetto, M. Troyer, and E. Vicari for useful discus- 
sions and information. 

1 M. E. Fisher and D. R. Nelson, Phys. Rev. Lett. 32, 1350 
(1974); D. R. Nelson, J. M. Kosterlitz, and M. E. Fisher, 
Phys. Rev. Lett. 33, 813 (1974); J. M. Kosterlitz, D. R. 
Nelson, and M. E. Fisher, Phys. Rev. B 13, 412 (1976). 

2 H. Rohrer, Phys. Rev. Lett. 34, 1638 (1975). 

3 K. Binder and D. P. Landau, Phys. Rev. B 13, 1140 (1976); 
D. P. Landau and K. Binder, Phys. Rev. B 24, 1391 (1981). 

4 K. Takeda and K. Koyama, J. Phys. Soc. Jpn. 52, 648 
(1983); J. Phys. Soc. Jpn. 52, 656 (1983). 

5 H. J. M. de Groot and L. J. de Jongh, Physica B 141, 1 

6 R. A. Cowley, A. Aharony, R. J. Birgeneau, R. A. Pel- 
covits, G. Shirane, and T. R. Thurston, Z. Phys. B 93, 5 

7 R. van de Kamp, M. Steiner, and H. Tietze-Jaensch, Phys- 
ica B 241-243, 570 (1997). 

8 R. J. Christianson, R. L. Leheny, R. J. Birgeneau, and R. 
W. Erwin, Phys. Rev. B 63, 140401(R) (2001). 

9 K. Ohgushi and Y. Ueda, Phys. Rev. Lett. 95, 217202 

10 A. Cuccoli, G. Gori, R. Vaia, and P. Verrucchi, J. Appl. 
Phys. 99, 08H503 (2006). 

11 P. Calabrese, A. Pelissetto, and E. Vicari, Phys. Rev. B 
67, 054505 (2003). 

12 C. J. Gorter and T. Van Peski-Tinbergen, Physica 22, 273 

13 H. Matsuda and T. Tsuneto, Prog. Theor. Phys., Supp. 
46, 411 (1970). 

14 K.-S. Liu and M. E. Fisher, J. Low. Temp. Phys. 10, 655 


15 M. Holtschneider, W. Selke, and R. Leidl, Phys. Rev. B 
72, 064443 (2005). 

16 C. Zhou, D. P. Landau, and T. C. Schulthess, Phys. Rev. 
B 74, 064407 (2006). 

17 A. Pelissetto and E. Vicari, |cond-mat/0702273"1 (2007). 

18 G. Schmid, S. Todo, M. Troyer, and A. Dorneich, Phys. 
Rev. Lett. 88, 167208 (2002). 

19 M. Kohno and M. Takahashi, Phys. Rev. B 56, 3212 

20 S. Yunoki, Phys. Rev. B 65, 092402 (2002). 

21 A. D. Bruce and A. Aharony, Phys. Rev. B 11, 478 (1975); 
D. Mukamel, Phys. Rev. B 14, 1303 (1976); E. Domany 
and M. E. Fisher, Phys. Rev. B 15, 3510 (1977). 

22 R. Leidl and W. Selke, Phys. Rev. B 70, 174425 (2004); 
R. Leidl, R. Klingeler, B. Biichner, M. Holtschneider, and 
W. Selke, Phys. Rev. B 73, 224415 (2006). 

23 A. W. Sandvik and J. Kurkijarvi, Phys. Rev. B 43, 5950 
(1991); A. W. Sandvik, J. Phys. A 25, 3667 (1992). 

24 O. F. Syljuasen and A. W. Sandvik, Phys. Rev. E 66, 
046701 (2002). 

25 K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 

26 P. Sengupta, A. W. Sandvik, and D. K. Campbell, Phys. 
Rev. B 65, 155113 (2002). 

27 L. Onsager, Phys. Rev. 65, 117 (1944). 

28 J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 

29 D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett. 39, 


1201 (1977). 

P. Minnhagen and H. Weber, Physica B 152, 50 (1988). 
H. Weber and P. Minnhagen, Phys. Rev. B 37, 5986 (1988). 

K. Harada and N. Kawashima, Phys. Rev. B 55, R11949