Table Of ContentDo dipole correlations exist in two-dimensional spin ice?
M. I. Ryzhkin∗
Institute of Solid State Physics of Russian Academy of Sciences
2 Academician Ossipyan str., Chernogolovka Moscow District, Russia, 142432
(Dated: January 29, 2014)
Here we study the statistical properties of two-dimensional spin ice in its ground state by the
Monte Carlo simulation method. Using a new sampling algorithm, we show that the short-range
ice rule in two dimensions gives rise to long-range but not dipole-like correlations, to non-Gaussian
probability density function for magnetization and to non-extensive conditional entropy (entropy
with given value of magnetization).
4
PACS numbers: 75.10.Hk, 75.40.Mg, 05.10.Ln
1
0
2
Introduction. - Originally, the name ”spin ice” was
n
given to compounds like Ho Ti O ,Dy Ti O which
2 2 7 2 2 7
a
demonstrate unusual magnetic ordering [1]. The mag-
J
neticionsinthesecompoundssitattheverticesofregular
8
tetrahedrons linked into a three dimensional pyrochlore
2
lattice. Owing to strong anisotropy, the spins of mag-
] netic ions can be directed along local anisotropic axes
h
connecting centers of nearest tetrahedrons. The ground
c
state is characterized by the ice rule: two spins of each
e
m tetrahedronaredirectedtowarditscenter,andtwoother
spins, from its center. Subsequently, the term spin ice
-
t was applied to any physical system with spins oriented
a
according to the ice rule (two spins in and two spins out
t
s rule). The spin ice shows such amazing properties as
.
t an exponentially degenerate ground state [2], long-range
a
m dipole correlations [3–6], self-screening of dipole interac-
tion [7, 8] and emergent magnetic monopoles [9, 10]. It
-
d isalsoconsideredasaclassicalanalogueofquantumspin
n liquid [11].
o Thetwo-dimensionalimplementationofspinicecalled FIG. 1. Ising spins in two-dimensional spin ice sit on each
c an artificial spin ice creates new possibilities for infor- bond of square lattice and are constrained to point in only
[ towardoroutofthenearestneighborvertices. Intheground
mation processing technology [12, 13]. The correspond-
state configurations, all vertices satisfying the ice rule (two
1 ing theoretical model belongs to the class of six-vertex
spins in and two spins out of each vertex) have the same
v models allowing exact solutions (for some specific val-
energy. Theloopsofspinsalignedalongloopsarehighlighted
6
ues of parameters) [14–16]. The availability of analytical in bold line.
7
and experimental results makes two-dimensional spin ice
1
a unique testing area for debugging of numerical algo-
7
. rithms. That explains our initial purpose: to develop
1 state,thereareonlycorrectverticeswhichsatisfyicerule
a new sampling algorithm for the Monte Carlo simula-
0 (see Fig.1). It is easy to see that in a ground state, the
tionofspiniceusingthetwo-dimensionalspinicemodel.
4 flip of a spin results in violations of ice rule for two ver-
However, the obtained results were so unexpected that
1 tices: with three (one) spins in and one (three) spin out
: thestudyofstatisticalpropertiesoftwo-dimensionalspin
v from the vertex. These vertices carry an effective mag-
ice became the main substance of the paper.
i neticchargeandhavefinitecreationenergy. Bythisrea-
X Let us describe the initial and modified purpose in de-
son, a usual sampling algorithm (random choice of spin
r tail. The model of two-dimensional spin ice can be de- and its flip with the probability defined by Metropolis
a scribed as follows. Ising spins sit at bonds of square lat-
function) becomes ineffective at low temperatures since
tice and are only allowed to be directed along the bonds
the simulation time grows exponentially as the tempera-
(see Fig. 1). All vertices satisfying ice rule (two spins in
ture decreases [17]. Moreover, it is not applicable for the
andtwospinsoutofeachvertex)havethesamestatistical
studies of the ground state at all.
weight or energy (taken as a zero), violations of ice rule
give some finite positive energies. Therefore, in ground To overcome this problem, special algorithms were de-
veloped. They are all based on searching for a closed
pathof spins alignedalong thepath (see, asan example,
theloopinfigure1)andonthefollowingcollectiveflipof
∗ [email protected] all spins in the loop [18–21]. Such collective flip do not
2
violate ice rule and strongly reduce the simulation time.
But they do not solve the problem completely. In fact,
the configuration space of spin ices has a very involved
structuresplittingintoseparatedtopologicalsectors[22].
Thetransitionsbetweenthesectorsarekineticallyslowed
down, in particular, they cannot be realizable by loop
flips. Indeed,loopflipsdonotchangemagnetizationand,
therefore, they cannot lead to transitions between topo-
logical sectors with different magnetization. This is the
main reason why we planned to explore a new sampling
algorithm which could provide spin configurations from
different topological sectors.
It is generally believed, that the simple two-
dimensional model described above keeps all essential
features of three-dimensional spin ice, for example such
properties as (1) dipole correlations, (2) extensive condi-
tional (an entropy at given magnetization) entropy and
FIG. 2. Two-dimensional histogram G(M ,M ) for the clus-
(3) Gaussian density distribution function of magnetiza- x y
ter of 1986024 spins as a function of magnetization. Total
tion. However, our results reveal that in two dimensions number of configurations is about 107.
conditional entropy is non-extensive, probability density
functionofmagnetizationisnon-Gaussian,andthesespe-
cial features are consequences of long-range non-dipole of the cluster. Note that N < N(cid:48) because we exclude
correlations. two outer layers of vertices. The numbers of spins in the
Method. - A specific feature of our approach is a new interior domain are taken from N = 40 to N = 77224,
sampling algorithm consisting of two stages. At the first andthenumberofconfigurationistakenaboutT =105.
stage, we define absolutely by chance all spins of cluster. In individual cases, we take N = 1986024 and T = 107.
As a result, we get an absolutely random configuration From the results of simulation, we calculate one and two
with a great number of incorrect vertices, where ice rule dimensional histograms (or sampling probability density
is broken. At the second stage, we use the Metropolis functions)g(M ),G(M ,M ),andthefirstsamplingmo-
x x y
annealing algorithm to reduce the number of incorrect ments
vertices to zero. This means that we take a random spin
and try to turn it with the probability defined by the (cid:88)T
E(Mm )=T−1 Mm (t), (2)
Metropolisfunction[17]. Thisprocedureisrepeateduntil x,y x,y
a ground state configuration is obtained. t=1
After the second stage, we get a ground state configu- (cid:88)T
E(Mm Mn )=T−1 Mm (t)Mn (t), (3)
ration which can be used for calculation of physical vari- x,y x,y x,y x,y
ables. We repeat the two stages procedure many times t=1
togetsettledandconsistentresults. Thissamplingalgo-
form,n=1,2,3,4. Allquantitiesarecalculatedasfunc-
rithmcreatestheconfigurationsdistributedoverconfigu-
tionsofN,andforfutureanalysiswealsocalculatethem
rational space more uniformly than the usual ”loop pro-
for the non-interacting model whose properties are well
cedure” [18–21]. The computer implementation of our
known.
algorithm is simple, rather fast and it can be directly
Results - First, we have found that simple generators
modified for many processor systems by a simple accu-
of pseudo-random numbers result in artificial effects and
mulation of results from individual processor.
correlations. To test the employed pseudo-random num-
With this algorithm, we generate T ground state con-
ber generators, we compared the simulation performed
figurationsforsquareclustersoftwo-dimensionalspinice
with them and the one made with files of genuine ran-
with N(cid:48) spins. We use free boundary conditions and ex-
dom numbers [23]. Unfortunately, one could not use
clude several outer layers from further computations to
only these files to perform long Monte Carlo simula-
avoid boundary effects. For each configuration, we cal-
tion, because they are rather short. We have found that
culate magnetization according to the formula
the pseudo-random number generator called Mersenne
twister gives a sufficiently good results, moreover it has
N
(cid:88) an enormous period of 219937 that allows very long sim-
M(t)= σ (t)e (1)
k k
ulations [24].
k=1
Our results are given in Figs.2-4. In Fig.2, we give a
Here σk(t) = ±1 is the spin on the bond k at the con- two-dimensional histogram G(Mx,My) as a function of
figuration with integer index t, e are unit vectors equal magnetization. Thecalculationsofthefirstsamplingmo-
k
to e ,e for the horizontal and vertical bonds respec- ments show that odd moments equal zero E(M1,3) = 0,
x y x,y
tively, and summation is over all N < N(cid:48) interior spins and the even ones are isotropic E(M2,4) = E(M2,4), in
x y
3
important result implies violation of the central limit
7
theorem for spins interacting according to the ice rule.
Therefore, it is natural to assume that the probability
6
density function for magnetization can differ from the
5 (2) Gaussian one.
Thisassumptionwillbethoroughlyexaminedandcon-
N4 firmed below. However, it is first reasonable to con-
/
Mx sidertherelationbetweenEq.(4)andassertionaboutthe
D3 dipolar character of the correlation function [3–6]. The
latter, strictly speaking, concerns the spin ice mode in
2
three dimensions. But its extension for two dimensions
is straightforward and quite common. In our case, the
1
(1)
above statement can be written in the form:
0
0 2 4 6 8 (cid:68) (cid:69) δ −2n n
N x104 Sαβ(r)= Mα(r)Mβ(0) ∝ αβ r2 α β (5)
FIG.3. Dispersionsofmagnetizationasfunctionofthenum- where nα =rα/r. Let us show that the result expressed
berofspinsforthenon-interactingspins(1)andforthespin by Eq. (4) disagrees with that of Eq.(5). To this end we
ice (2). estimate the variance (below summation is restricted to
the horizontal bonds):
addition the mixed moment E(MxMy) = 0. The accu- (cid:34)(cid:32)N/2 (cid:33)(cid:32)N/2 (cid:33)(cid:35)
(cid:88) (cid:88)
racyofthesestatementsbecomeshigherasthesimulation DM =E σ σ = (6)
x i k
time T grows (that has been checked up to T =107). In
i=1 k=1
Fig.3,wegivethedispersionsofmagnetizationasafunc-
N/2 N/2 N/2
tion of the number of spins N for both the spin ice and E(cid:88)σ2+E(cid:88)σ σ ≤ N(cid:104)1+(cid:88)E(σ σ )(cid:105) (7)
forthesystemwithoutanyinteractionbetweenspins. In i i k 2 i 0
i=1 i(cid:54)=k i
Fig.4, the dependence of ∂lnG(M)/∂M2 upon N is also
shown for both the spin ice and for a paramagnetic sys- The quantity E(σ σ ) is the spin-spin correlation func-
i 0
tem. Now let us discuss the obtained results. tion S (r ,0). The use of Eq.(5) yield the following es-
xx i
Discussion - Fig.2 reveals that a density of configura- timate:
tion is isotropic and has a peak at zero magnetization
√
M = 0. These results were expected from general con- (cid:32) N (cid:33)
N (cid:90) dr DM
siderations, and can be regarded as a proof of the ab- DM ≤ 1+a ⇒ x ≤a +b lnN (8)
sence of essential errors in the simulation procedure. In- x 2 1 r N 2 2
deed, the plot in Fig.2 is very sensitive to any attempts b1
to accelerate the configuration generation algorithm via
where a and b are positive constants. The right
1,2, 1,2
insertion of some faster but determinate elements. For
hand side of inequality (8) increases slower than the
example, replacement of the random choice of spins for
power function (4). To reproduce the dependence (4),
possible flip by sequential scan of all spins results in a
one has to assume that at large distances the correlation
shift of the maximum from zero. Further, it becomes function behaves as S(r) ∝ r−2+β with the parameter
clear from the histogram shape that E(M1,3) = 0 and
x,y β = 0.74, i.e., correlations decay much slower than dic-
E(Mx2,4)=E(My2,4),whichwasalsotestedbydirectcal- tatedbydipolarlaw(5), namelyasS(r)∝r−1.26. Thus,
culations. Accordingly, the histogram peak for the non- Eq.(4) is at variance with the assertion on the dipolar
interacting system also locates at M=0, possesses axial character of the correlation function, at least in two di-
symmetry, but exhibits a more narrow peak. mensions.
From Fig.3 one can see that for a noninteracting This variance should be discussed in more detail. The
system the quantity D(Mx)/N = const as should be result expressed by Eq.(5) is based on theoretical inves-
expected for noninteracting and, therefore, for non- tigations of two types. First, it is often stated that such
correlatedspins. Thisresultisinagreementwithcentral a behavior corresponds to the exact solution obtained in
limiting theorem for non-correlated spins. But for spin [15]. This is just an extension of the results of [15] be-
ice with a short-range interaction determined by the ice yond its applicability domain. In [15], the exact solution
rule, the quantity D(Mx)/N (cid:54)= const. We have found for the correlation functions was derived for the six ver-
that the dependence of D(Mx)/N upon N can be ap- tex model with the parameter ∆ = 0, whereas for the
proximated by power-like function two-dimensional spin ice ∆ = 1/2. For exactly solvable
D(M )/N =a+bNα (4) models such a difference can be critically important and
x
can change the result qualitatively. Second, references
wherea,baresomepositiveconstantsandα=0.37. This aremadetoworks[3–6]whichinvolvetherepresentation
4
is proportional to N, the entropy s(M) also proves to be
4
proportionaltoN;i.e.,itisanextensivequantity. Onthe
contrary, for system with the ice rule, K ∼ N−1.23 (see
6 line2)andtheconditionalentropyequalss(M)∼N0.77,
i.e., non-extensive. The nature of this non-extensivity as
wellastheresultexpressedbyEq.(3)areexplainedbythe
8 (1) long-rangecharacterofcorrelationsinspinice. Itshould
K) be emphasized that the non-extensivity is due not to the
n(
l long-range character of the interaction (it is extremely
10
short range) but exactly to the long-range spin-spin cor-
relations.
(2)
12 Ourresultsalsoshowsthattheiceruleleadstobroad-
ening of the probability density function for magnetiza-
tioncomparedtotheGaussianoneforthenoninteracting
14
spins. The test of the hypothesis about Gaussian char-
4 5 6 7 8 9 10
ln(N) acterofdensitydistributionfunctionbyPearsonmethod
yieldsthePearsonparameterχ2 =3·103. Thatisgreater
FIG. 4. Dependence of K = ∂lnG(M)/∂M2 on the number than the quantile χ2 for all confidence probabilties γ
γ,97
of spins. This quantity defines the entropy as a function of (weused100intervals,thenumberofthedegreesoffree-
magnetization (see Eqs.9). dom was k = 100−2−1 = 97). That means that the
hypothesisoftheGaussiancharacteroftheexperimental
density cannot be accepted at any confidence probabil-
ity. We found the experimental distribution function for
of the magnetization distribution function in the form
the spin ice to be closest to the self-similar distribution
f(x,t) = t−Hp(xt−H) with the parameter H ≈ 0.68.
K
P{M}∝exp[s(M)], s(M)=s(0)− M2 (9) In this case, the role of the discrete time is played by
2 the number of spins in the system, t ∼ N. Distribution
functions of this kind appear in the theory of anomalous
wheres(M)istheentropyofthesystematgivenmagne- diffusion (normal diffusion corresponds to H ≈0.5).
tization. However, the use of Eq.(9) actually implies the In conclusion, we summarize the basic results of this
use of the Gaussian distribution of magnetization from work. The behavior of the variance indicates correla-
the very outset. And as was mentioned above, the result tions which decay slower than dipolar correlations (see
expressed by Eq.(3) forces us to question the validity of Fig.3 and the estimates). The conventional entropy is
thisassumption(atleastintwodimensions). Formally,a a nonextensive quantity (see Fig.4). The magnetization
non-Gaussiandistributionimpliesthatthefunctions(M) of two-dimensional spin ice does not obey the central
has a more complicated form than that given in Eq.(9), limit theorem; i.e., the distribution of the magnetization
it can include higher terms of the expansion in powers of isnon-Gaussian(seeFig.4andtheresultsofthehypoth-
the magnetization [25] or even have no expansion at all. esis testing). All these results are the consequences of
Inapplicability of formulas (9) in two dimension is long-range spin-spin correlations.
confirmed by the following result. If one presents the
density of configurations by formula (9), the quantity
∂lnG(M)/∂M2 will play the role of the experimental ACKNOWLEDGMENTS
value of the coefficient K in Eq. (9). The dependence
of this quantity on the number of spins in the cluster is I am grateful to S. N. Molotkov and I. A. Ryzhkin for
showninFig.4. Inthiscase,fornoninteractingspins,one theformulationoftheproblem,supportofthework,and
has K ∼N−1 (see line 1). Since the total magnetization discussion of the results.
[1] M.J.Harris,S.T.Bramwell,D.F.McMorrow,T.Zeiske, Phys. Rev. Lett. 93, 167204 (2004).
and K. W. Godfrey, Phys. Rev. Lett. 79, 2554 (1997). [6] C. L. Henley, Phys. Rev. B 71, 014424 (2005).
[2] A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, [7] M.J.P.GingrasandB.C.denHertog,Can.J.Phys.79,
and B. S. Shastry, Nature 399, 333 (1999). 1351 (2001).
[3] R.Youngblood,J.D.Axe, andB.M.McCoy,Phys.Rev. [8] S.V.Isakov,R.Moessner, andS.L.Sondhi,Phys.Rev.
B 21, 5212 (1980). Lett. 95, 217201 (2005).
[4] D.A.Huse,W.Krauth,R.Moessner, andS.L.Sondhi, [9] I. A. Ryzhkin, JETP 101, 481 (2005).
Phys. Rev. Lett. 91, 167004 (2003). [10] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature
[5] S.V.Isakov,K.Gregor,R.Moessner, andS.L.Sondhi, 451, 42 (2008).
5
[11] L. Balents, Nature 464, 199 (2010). (1979).
[12] R.F.Wang,C.Nisoli,R.S.Freitas,J.Li,W.McConville, [20] D. J. Adams, J. Phys. C: Solid State Phys. 17, 4063
B.J.Cooley,M.S.Lund,N.Samarth,C.Leighton,V.H. (1984).
Crespi, and P. Schiffer, Nature 439, 303 (2006). [21] G.T.BarkemaandM.E.J.Newman,Phys.Rev.E57,
[13] A. Y. Kitaev, Ann. Phys. 303, 2 (2003). 1155 (1998).
[14] E. H. Lieb, Phys. Rev. Lett. 18, 692 (1967). [22] L. D. C. Jaubert, M. J. Harris, T. Fennell, R. G. Melko,
[15] B. Sutherland, Phys. Lett. A 26, 532 (1968). S. T. Bramwell, and P. C. W. Holdsworth, Chem. Phys
[16] R. J. Baxter, Exactly Solved Models in Statistical Me- 43, 329 (1979).
chanics (Academic Press, 1982). [23] G. Marsaglia, http://www.stat.fsu.edu/pub/diehard
[17] K.BinderandD.W.Heerman,Monte Carlo Simulation (1995).
in Statistical Mechanics (Springer, 2010). [24] M. Matsumoto, http://www.math.sci.hiroshima-
[18] A.RahmanandF.H.Stillinger,J.Chem.Phys.57,4009 u.ac.jp/ m-mat/MT/emt.html (2007).
(1972). [25] I. A. Ryzhkin and R. W. Whitworth, J. Phys.: Cond.
[19] A. Yanagawa and J. F. Nagle, Chem. Phys 43, 329 Matter 9, 395 (1997).
