7 . ~.,—~. -- “WM-9:1-
. ,..|'. I,"~.- . J , u'n ' ' I"'n1
u. _: .1..‘."‘ ' " - (".II‘I2I"',"‘.~..'5.I‘-'3fi- "352:" "MI 2““.
1" .\'-I1-' ,11'3 ' 11.11.1116“ '1’?» ' 'II 11"1'1‘1." -‘ II
N .4. : 1.1w~.--.wmw ,-

   
     
  
      
   
  

 
 
 
  
 
  
 
   
    
  
 
 
   
 
  
 
   
 
  
  
  
      
  

   

     
   

:" I' Riga. ‘ "II! ‘J L}! El. 2-":‘1 I ' I ‘ ‘~‘ .‘ ‘ u I
111.” "11?". p' .1..#z1:1‘t:11u'\h' $11.17;; 44"”: ‘5 11 11.1.3)“. \ .
.' 1' 1".""I‘“.I' ..-"',,-.-:- 11w- .,- ..-' .- , .-
- . I : "I :11"... 'I ' 'x "J." .\‘1 i 'l? >
" I11'II"‘-‘:1'II:.I111111III'IA1-"11'5'111‘Ig11'II ‘ i" “1'3"“? 1‘11" 5.”.
ilfif'”. .I' "'.A ' ,DI ‘71:}; 1,(.( .I. '. W C". I :.I.I_ ’ : 2'
I~m . 1.1.«r-yafiwymnwc- g.
g1! I11-||"1(.$p‘h&" 2‘31“ “‘11. :13. 'fii-Hfi‘jith 1 1“! fr “1.1, l.: . m-Z-R‘g
VI ”.l‘l 1" 1") "V I...
' u.‘:.. 1'11? 111-] II" ‘1 "1r". 1 ‘1
11 1" 1 " 11111 .

’13

  

1'35. .
\IW I "( I..I.)-‘II- xi. "'2“. .
8'" 1‘1 I 11 I"'*IIII..I";<““‘""II’ ' 1|
-,QM‘1I¢$W11%I'H
., 1" ‘1 11111- 1'11 f. 11131.- W. 1 '
1 II..- .. -.

'.’ ‘ I~ _' II':|.|'
I 1‘ $0 II
. 1‘

  
 

-.-. , 1‘ '51? I11 i'r“ 111:. £40.
.1.- 1.1-.
.11 1 LI“? ' 1...” 'Il1l1l‘\ll1lfil;' . l;‘\'r {S'W
1111.314 .1 '1“i".u"5"1:11'"§l ”I11“ 1 h":

1 'I I (‘N'IRHSP {115.1.” ”1:. 1:13;” I333} .111“
:1 ' "1 "1" "’ 1-1-3.- ’-‘ I13. .&:1.Ii?fis
11. I ‘, 9-"-‘I"'."I".1‘ , ,.
.‘11'I‘1I..1.{‘V4‘I~1I‘ U111: 1'11 (I 9%?" 11:“1‘1510.’ 1.\11‘( '1‘??? J'Izr"§.fi“ '1'. ‘. .‘\‘1‘E.I~'.I “Q.
1 ' .111'25': :JJ‘IJR‘I‘ ”‘ '1‘ h”! \".1 1m ”1““... “3:7?! ”1.. N\'.\’(‘ v: '1'. 1% "1‘”,',-"‘1: 1

1111"”‘1l1 I51" I'II“I .

{13133110601‘ . ."_‘I. H "ILITI .‘;"- ”1.0111313;
"-13 .' .h "'.' 1_ 1.51.0.1 1f 1' I
. It. .- 1. 1-! .-
X1‘L"I\"2‘11 ".1" ' "1:111 "_..V1I , ' III",.*..I’- "' 1 1 1.).--
. ““1

11112-11)" ' J1“.

  

1 I1

  
 
  

   
    

.‘s
.0

     
      
      
 
     
      
  
    
  
   
  

  

1I | 1.,W .gw‘ U ’ng .4 A -.
H“ 1“; ”7‘ I I” .,'.‘.:"‘_
. 11“" $11.45.". ':21'.1.‘:'I'E."‘1‘I}:;IIH" ' ' ‘II'h {- if“: , :1‘.‘.‘ :1‘E‘1g: ‘JKI;I.I‘,%‘%‘§L "' 'g’flr‘fib’
, ." I " "1 11 1‘11 "' "'-'-.IgI1'5}1'-111‘\$‘,1"'1‘HI mh‘fl " 1‘15“."7-"1‘ “1" "fl."
. » I ‘..I|'III' H. h -_I. 1‘ y - .II ’ '1‘ . p I I . .
: Ii..l'r"'I"I')>1 >1” .1”- l' '1.'1I " 11'!“ III"I‘1V‘.‘SII""II' ':;’:"P’€ '1‘“ ‘tfh‘fi'l ', [:1 1.1.7:" gzi'11‘v ”it" I} ,1 I"? ‘
£331.11 7-1-1 3.'I‘.'.3.?1\'I":.. "'7!- '1.I'I":I.-'."711' 1.. ‘1'”"6-‘1-I'31?I'I'-w'.*."- ~" I ‘g,':...-‘.':~;.:.".. .4-
..II.' .1101. 1 .111“ .1‘1','-I..I..:.Ij, -11 1 I \'l‘1".\‘ 1.11. ~1.’ ‘. , -'-‘ ., .- .I'«-. .- "11'“. «V. '4
- . - I - .. 2 1-11-1111: ”11-- - -. '-‘.':~.." i w .-
i' 1111' ' ".‘v- ""11' " “ “'1 "'-1""II""' "‘ 7913116111'11.'3"!£"-.'-1v~-'II‘I
"'3'" II 11 ’.' "'I‘IIII 1 I: '1‘ 1111‘" "1 ‘-".1'.".I-"11"' ‘1‘" "I "111\"""I"""II""-'\' .8I
., '1 I 1| I 1 _. .D.‘ .1 1 .' .- ._I. 'l' I- .1 .l) .' -. n
I"" 1.1.. I 1. 1' H I. ‘11l'llll .1. 1.11 1‘.. '11I ' ." 'J‘Mflo ‘4‘31‘1'Q' ,1?" ‘ch 4,11!
,Iuwm1I1N-nvnwm... I1QU‘ mHuMHw?
"'I' I "'II \" II "'I 11' h ‘ ' ‘ I' ‘1 ‘
I II IA .1 {I I ' 0 ‘ I I
.ll' 1.-l1h- In" - 1. 1' . - ‘9 .

    

 
 

”2;. 1-.1' .131
-., I‘M. .1
. Iv ..\~ I " HI. . '\1 '. " 111" MIN! 1 IIII'T
"-" -' I'I ‘ .I II" II ','1' 'I.'I" II.‘\'1I"I~h'I ": I“. I" 1 I. ‘9‘1'-;-“‘ LYN-1...“; I I" IIII‘". 1. 'I. )1?‘XI1".
“'-' '.;2. - ' ' " '1 1- :I-I. - -‘
. I '11.: .1 1 .

    

I ' 1'.’ "It I II
I III1'1'1I1-"IIII “1'. . II1‘112xf111‘II"II"

 

III” III III III III

THESlg‘ W3 ‘1293101667

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I
LIBRARY l
Michigan State
University
\_-1_, “:1

 

 

f

 

This is to certify that the

thesis entitled
PROBLEMS IN ORDERED AND
DISORDERED ISING SYSTEMS

presented by

JOHN MARSHALL THOMSEN

has been accepted towards fulfillment
of the requirements for

 

Ph.D. degree in PhYSiCS

WW/es

Major pro ssor

M. F. Thorpe

Date I? ”up ({1ng

0-7639 MSU is an Affirmative Action/Equal Opportunity Institution

 

MSU

 

 

 

RETURNING MATERIALS:
Place in book drop to

 

LJBRAfiJES remove this checkout from

.—3—_ your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

_ ._‘__7._____ r_ -,

PROBLEMS IN ORDERED AND DISORDERED ISING SYSTEMS

By

John Marshall Thomsen

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

1984

oi

2/9/37"

hp}

ABSTRACT
PROBLEMS IN ORDERED AND DISORDERED ISING SYSTEMS
By
John Marshall Thomsen

A wide range of Ising systems are studied. In one dimension,
attention is focussed on calculating the transverse susceptibility.
By combining transfer matrix techniques with linear response theory,
it is shown that Xi can be calculated for pure and dilute spin S
Ising chains. The low temperature behavior in the presence of
a parallel field is discussed. The temperature dependent probability
distribution of internal fields, P(h), is defined. and it iS shown
to yield several exact thermodynamic results. This function is
studied on a variety of Ising ferromagnets, with exact results
available in one and two dimensions as well as on the Bethe lattice.
Three and four dimensional lattices are studied with Monte Carlo
simulations. It is shown that P(h) is very sensitive to short
range order. This order is in turn sensitive more to the lattice
dimension and less to the coordination number. P(h) is also studied
for the Sherrington-Kirkpatrick infinite range spin glass model. An

analytic result is available in the paramagnetic phase while Monte

Carlo simulations are employed in the spin glass phase. The de—.
velopment of the dip at h = O is investigated in some detail.
Finally, some speculations are made on the possibility of extending

the P(h) formalism to non-Ising systems.

ACKNOWLEDGMENTS

It is a pleasure to acknowledge valuable discussions with
Drs. T. A. Kaplan, N. Rivier, D. H. Lee, S. de Leeuw, J. Bonner,
and J. S. Thomsen. I am in particular indebted to Drs. J. Parkinson
and S. D. Mahanti each for a series of discussions which included
important insights and probing questions.

E. Garboczi, H. He, D. Heys, Dr. D. Sahu, and P. Murray have
all been good sources of advice during the process of converting
general ideas into concrete results. I

I would also like to thank S. Rosenbrook, J. King, and especially
Dr. J. Kovacs, whose efforts allowed me to concentrate more of my
time on research than might otherwise have been possible.

I am indebted to Drs. D. Stump, J. Bass, S. D. Mahanti, and
M. F. Thorpe who, as members of my guidance committee, must tolerate
my writing for another 150 pages. As my research adviser, Dr. Thorpe
has suggested to me more interesting problems than I could possibly
have had time to solve. His insight helped guide me through this
wide assortment of problems to produce some coherent results. I
should also point out that while the preprint in Appendix A was
a collaborative effort among myself and Drs. Thorpe, Sherrington,
and Choy, most of the writing was done by Dr. Thorpe. I have in
addition benefitted from some correspondence with Drs. Choy and

Sherrington.

ii

Finally, I should like to thank Amelia Hefferlin whose en-
couragement and assistance have helped bring this project to com-
pletion. No less important have been her occasional reminders
that Physics is not all there is to life (which is easy for an
Ecologist to say).

Most of this work was performed under the auspices of an

Exxon Fellowship.

iii

TABLE OF CONTENTS

List of Tables ...... . ..... . ...... . ............... ........ ...... .-V
List of Figures ................................................ ~Vi
1. Introduction ................................................. l

2. Transverse Susceptibility in One Dimensional Chains..........6

3. P(h): The Local Field Probability Distribution Function ..... 26
4. P(h) for Bethe Lattices ...... . .............................. 33
5. P(h) for Two Dimensional Ferromagnets ....................... 41
6. P(h) for Ferromagnetic Systems of Higher Dimension .......... SO
7. The Sherrington-Kirkpatrick Spin Glass Model ..... ...........58
8. P(h) for the Sherrington—Kirkpatrick Model..................68
9. Future Work ................ .............. .......... . ........ 82

Appendix A. Local Magnetic Field Distributions I.

Two-Dimensional Ising Models.................85

Appendix B. Some 2D Correlations .............................. 124

Appendix C. Details of Monte Carlo Simulations................128

List of References ....... ...... ..... . ......... . ............... 135

LIST OF TABLES

Table 1. <0001>T as obtained from Monte Carlo simulations
c
and compared to high temperature series results .......... ...53

Table 2. Nearest neighbor locations for three dimensional lattices..134

Figure

Figure

Figure
Figure
Figure
Figure
Figure

Figure

Figure

Figure

Figure

Figure

Figure

LIST OF FIGURES

1. The inverse of the transverse susceptibility per
spin for a dilute (c=0.25) Ising chain in the pre—
sence of three different parallel fields ............... ....25
2. Section of a Bethe lattice for z=3 (upper) and
=4 (lower) ................................................ 34
3. P(h) for Bethe lattices of coordination 3, 4, 6, and w ..... 40
4. P(h) for the ferromagnetic square net ...................... 43
5. P(h) for the ferromagnetic honeycomb lattice ............... 45
6. P(h) for the ferromagnetic triangular lattice .............. 46
7. Generating the 3-12 lattice from the honeycomb lattice.....47
8. Showing P(h) at TC for the honeycomb (upper) and
3—12 (lower) lattices. ................ .... ................. 49
9. Extrapolation of <0001> to the infinite lattice limit
for a simple cubic at TC................. ....... . .......... 52
10. P(h) for four different three dimensional ferromagnets....55
11. P(h) for simple cubics of dimensions one through
four and for the mean field limit ......................... 56
12. The phase diagram of the Sherrington-Kirkpatrick model
(reproduced from Sherrington and Kirkpatrick 1975) ....... 61
13. P(h) from the analytic result for the Sherrington-

Kirkpatrick model in the paramagnetic phase (JO < J) ...... 72

vi

Figure

Figure

Figure

Figure

14.

15.

16.

17.

A comparison between the Monte Carlo simulation of
the Sherrington-Kirkpatrick model and the analytic
result: T = ZTSG.. ......................... .. ........... ..75
A comparison between the Monte Carlo simulation of
the Sherrington-Kirkpatrick model and the analytic

result: T = T ...; .............................. . ....... .76

86
Results of the Monte Carlo simulations of the
Sherrington-Kirkpatrick model: P(h) at four different
temperatures. The T=O data has been symmetrized .......... 77
P(O) as a function of temperature in the spin glass

phase of the Sherrington-Kirkpatrick model. Open

circles indicate simulation results. The dashed line

‘W

represents the linear hypothesis .......... . ............... /9

vii

1. Introduction

A fundamental and much studied problem in condensed matter
physics is the explanation on a microscopic level of the magnetic
properties of various materials. As might be expected, there is
frequently a competition between building a model Hamiltonian that
describes a system with accuracy, and building a model Hamiltonian
which is tractable. Fortunately, there exist a variety of systems
which are described reasonably well by a simple spin Hamiltonian.

For a system of well localized effective moments, one general

form a spin-spin interaction may take is

H12 = film-'3’2 (1)

where 3+is a 3X3 interaction matrix and S1 and S2 are three compo-
nent operators describing spins located at sites 1 and 2. The matrix
elements of UIare related to the overlap of the electronic wave
functions associated with the two spins. (For a justification of
this Hamiltonian, see White 1983.) In certain cases, the Jzz term

may dominate the matrix, reducing the interaction to the Ising (1925)
form I

z z
12 - -J8182 . (2)

Systems well described by this form over at least a restricted

H

temperature range include CoC12° 2NC5H5, a one dimensional ferro-
magnet (Takeda et. al. 1971), COCs3Br5, a two dimensional antifer-

romagnet (Mess et. al. 1967), and DyPOA, a three dimensional

antiferromagnet (Wright et. al. 1971).

The advantage of working with Hamiltonians built with Ising
interactions in the absence of a field is that every piece of the
Hamiltonian commutes with every other piece. In the presence of
external fields, the most general Ising Hamiltonian considered in
this work includes one non-commuting piece:

H = —%.Z.Jij3:8: - u hZZS: - uthisi (3)

ll .
1"] l l

where the Si are taken to be quantum spin operators. The matrix
of Jij's is assumed to be symmetric with vanishing diagonal elements.
The fact that the Si terms do not commute with the rest of the
Hamiltonian will in general cause problems; however hX will be taken
either to be zero, in which case the problems are eliminated, or
infinitessimal, in which case one can work around the commutation
problems.

The Ising Hamiltonian has been studied in a wide variety of
forms, and some of the more relevant studies are reviewed here.
When first introduced, the partition function of a pure spin % chain
with uniform nearest neighbor interactions was calculated (Ising
1925). Analytic results for the specific heat and parallel sus—
ceptibility are now available for spins up to 3/2 (Suzuki et. al.
1967). The site diluted version has been studied by, among others,
Matsubara et. al. (1973) for S=%, and Matsubara and Yoshimura (1973)
for 8:1 and S=3/2. In both cases, the parallel susceptibility and

the specific heat were calculated.

IH

The zero field transverse susceptibility of a pure spin 2

infinite chain was calculated by Fisher (1963), while Crest and

Rajagopal (1974) repeated the calculation in a finite parallel

field. Chapter 2 focusses on extending the transverse susceptibility
calculations to include higher spins and site dilution, both with

and without external parallel fields. Unusual effects arising

from dilution and from the application of the parallel field will

be discussed.

Traditional approaches to studying Ising models focus on de—
termining the free energy of the system and obtaining other thermo—
dynamic quantities by differentiation. The subsequent chapters of
this work deal with an alternative approach involving the distri—
bution of local fields. This function, P(h), has been much studied
in relation to spin glasses at zero temperature (see for instance
Palmer and Pond 1979 or Kirkpatrick and Sherrington 1978). However,
it will be shown in Chapter 3 that at all temperatures, the thermo-
dynamics of an Ising system may be reformulated exactly in terms of
the temperature dependent P(h). From P(h), one may obtain not only
the magnetization and energy (and hence the free ener8Y). but also
the inelastic neutron scattering cross section. Furthermore, since
P(h) describes what is happening locally, it can give some intuitive
feel for what is happening to the system. It was for these reasons
that a study of P(h) on various Ising systems was undertaken.

While the original motivation to study P(h) lay in spin glass
research, it is useful to first study this function on well under-
stood uniform Ising systems in order to obtain a basic understanding
of how it behaves. In Chapter 4, P(h) is calculated for Bethe

lattices with nearest neighbor ferromagnetic interactions. The

Bethe lattice allows one to make a fairly smooth transition from a
network of small coordination to the mean field limit, and hence
may provide some insight into the effect of coordination number
(or perhaps dimensionality) on this function.

The Bethe lattice provides well behaved analytic results, but
it is not a true lattice and the large fraction of spins residing
on its surface makes it unphysical. For this reason, it is useful
to look at the two dimensional Ising model. This model takes on a
special significance in statistical mechanics as it is one of the
few exactly solvable models exhibiting a phase transition at a
finite temperature. The partition function of a square net with
uniform nearest neighbor interactions was first derived by Onsager
(1944). He extended the transfer matrix technique introduced
by Kramers and Wannier (1941). Since then, a variety of solutions
have been developed. The Pfaffian approach taken by Montroll et.
al. (1963) seems to lend itself most readily to calculating cor-
relation functions, although it may not be the quickest way to
determine the free energy. This approach has since been generalized
to a wide variety of two dimensional lattices (see, e.g., Hurst
1963 and Stephenson 1964). These results are drawn on in Chapter 5
to perform exact calculations of P(h) for several two dimensional
Ising systems.

The general features of P(h) on these two dimensional lattices
are relatively insensitive to coordination number. However, there
is a strong dependence of P(h) on lattice dimension. This depen—

dence is explored in Chapter 6 using Monte Carlo simulations.

Results in three and four dimensions are compared to the mean field
limit. The models chosen, nearest neighbor ferromagnetic inter-
actions on diamond, simple cubic, body centered cubic, face centered
cubic, and hypercubic lattices, have all been studied in detail
with high temperature series. Since these series provide suffi-
ciently accurate estimates of TC for the purposes of this study
(see, e.g., Domb 1974), it is possible to concentrate on P(h) with-
out getting bogged down in details such as locating Tc'

In Chapter 7, attention is turned to a spin glass model, the
Sherrington-Kirkpatrick (infinite range) Ising model (Sherrington
and Kirkpatrick 1975). This model was first introduced as an
approximation to a spin glass system dominated by the long range
RKKY interaction. It was hoped that the infinite range interaction
would lead to an exact solution, although the solution in the spin
glass phase has turned out to be more complex than originally hoped
and is still under scrutiny (Parisi 1983). Nevertheless, this is
one of the more tractable spin glass Hamiltonians and one can in
fact obtain an analytic result for P(h) in its paramagnetic phase.
There are some distinctive features associated with the small h
behavior of P(h) and these will be elaborated on in Chapters 7
and 8.

Finally, in Chapter 9 further studies of Ising systems are
suggested as is the possibility of extending the P(h) formalism

to non-Ising systems.

2. Transverse Susceptibility in One Dimensional Chains

The one dimensional Ising model is attractive in that many
exact results for the system may be obtained through straightfor-
ward calculations. While this model does not exhibit any phase
transitions for T>O, there are systems, such as CoC12°2NC5H5, which
behave at least approximately like one dimensional Ising chains
(Takeda et. al. 1971). In this chapter, a generalized method will
be developed for calculating the transverse susceptibility of Ising
chains by using modified transfer matrices.

First, some background on the transfer matrix technique, de-
veloped by Kramers and Wannier (1941), must be introduced. One

starts with a Hamiltonian of the form

N N
_ _ A z z _ 1 z z z
HN ' Jiil Sisi+1 2““ 1:1(Si+si+1) ’ (1)

describing a one dimensional spin S Ising chain in the presence

of an external parallel field, H = h2 2. Defining S§+1 = S? en-
sures periodic boundary conditions.
The partition function QN is given by
QN = Triexp(-BHN>I (2)
N z z z N z z
_ - i
- éz g2...§z exp{B[JiEISiS1+1 + zuh i§l(Si+Si+l)]} . (3)
1 2 N I -
The

notation g2 implies a sum over all (28+1) possible values of
i

S: and B is as usual 1/kBT. The (28+1)X(28+1) transfer matrix,

T, is defined by its entries

_ z z A z z z
T 2 z — exp{B[JSiSi+l + guh (Si+S.+l)]} . (4)

l
SiSi+1

In terms of this matrix,

QN = g ...z T T ...T T . (5)

Z
I S; S; 1 2 2 3
In going from equation (4) to (5), explicit use was made of the
fact that all the operators in the Hamiltonian commute with each
other.

The operation in equation (5) is now recognized as matrix
multiplication, so that

N

QN = :2 (T )st2 <6)
1 1 1
= TrlTNl <7)
ZS+l V
= 2 Ah , (8)
3:1 3 .

where Aj is the jth eigenvalue of T. In the limit of infinite

N, the maximum eigenvalue, AM, dominates:

N
QN = (AM) (II-m) . <9)
Thermodynamic quantities may be obtained by the usual differentiation
of RnQN.

In the absence of an external field, for spin %,

eXP(BJ/4) €XP(-BJ/4)
T = (10)
exp(-BJ/4) exp(BJ/4)

with

A1

~2 cosh(BJ/4)

(11)

A2

2 sinh(BJ/4)

so that for large N one recovers the result (Kramers and Wannier
1941),

QnQN = N£n[2 cosh(BJ/4)] . (12)

For spin 1,
exp(BJ) 1 exp(-8J)
T = 1 1 1 (13)

exp(-BJ) 1 exp(BJ)

giving
A1 = cosh BJ + % + %q
A2 = cosh BJ + % - %q (14)
A3 = 2 sinh BJ
where
q = [(2 cosh BJ - 1)2 + 8]% . (15)
For large N,
anN = N£n[cosh BJ + % + %q] . (16)

With this background, it is now possible to proceed with the
calculation of the transverse susceptibility of an infinite, pure

Ising chain.

J. Phys. C: Solid State Phys. 16(1983) L237- 24-0. Printed in Great Britain

LETTER TO THE EDITOR

Transverse susceptibility in spin S Ising chainsi‘

M F Thorpe and M Thomsen:
Department of Physics and Astronomy. Michigan State Universuty. East Lansing. MI
18824. USA

Received 31 December 1982

Abstract. It is shown that the transverse susceptibility of a one-dimensmnal spin S Ising
system can be calculated by combining linear response theory With the :ranster matrix
technique. Closed form expressions for infinite chains are obtained for 5 = 5 and 5 = 1. and
numencal results are displayed for S =- l to S = i.

The perpendicular susceptibility of a spin é [sing chain has been calculated by several
methods leg. Fisher 1963. Katsura 1962. and GreSI and Rajagopal 1974). However. it
would appear that these methods have not yet been generalised to higher spin. The free
energy and parallel susceptibility have been found for higher spins using transfer matrix
techniques (Suzuki et al 1967). In this paper the spin dependence of the perpendicular
susceptibility will be studied. We use standard linear response theory (see for example
Grest and Rajagopal 1974) and apply it with transfer matrices to obtain exact closed
form expressions for the zero field perpendicular susceptibility for spin 3: and spin 1 as
well as numerical results for spin 3, 2. and 3.

Our Hamiltonian has the general form

H=HO+HI (1)

and we wish to find the expectation value of an operator .-'1 to first order in H I. where
HI... H I and A do not necessanly commute. It can be shown that If

 

Trfexp(-I13HI))A] == 0 I2)
then
I” "l’ I
-TrI I epr-(B - r)H..]HI~ expl - tH..; dr.~11
' -O
(1'1) = T d f 3
. Trlexm -I3Ho) l x
to first order in H 1.
In the Ising model for a closed cyclic chain we have
.V
Ho: -1 25555-. . (4)
_V
H1 = -_uh 215: I5)
-‘ Work supported in part by the National SCIence Foundation.
: Exxon Fellow.
© 1983 The InSIitute of Physics L237

C8—C2

10

L238 Letter to the Editor

where the S .- are quantum spins of magnitude 5 and site N + 1 is equivalent to site 1.
From equation (3) we find

.V
r B ./
(Sf) = TILL exp[-(fi- rJHnluh 2:1 5'} expl—rH0) drSfJ/x/ Tr[exp(-/3H.,)]. (6)

The tracing operation ensures that only the i = 1' term in the sum contributes. F urther-
more, we n0te that all but two terms of Ho (which are -.IS;'_1.S‘,z -JS,’S,‘- 1) commute
with 5}“ so that equation (6) can be rewritten as

:x_. VV‘

55-15311 3v

_ , , TSf—sz-i Tgf- 57" T52-l5,:-2 . . . ngsg/‘Tr[exp(-l3Ho)] (7)

l

where 25,: stands for the sum over the 23 i- 1 possible values of S 5. We have defined the
standard transfer matrix (Suzuki er al 1967)

Tsfsf’i = exp(flJSES.-‘-l) (8).

and have defined a new transfer matrix

13
T3: 5: I = (uh 2L exp[.l(fi- r)(5f-15f + Sfo-1)]

X Sf exp[Jr(S,-‘-IS;' + Sfo-1)] drSf. (9)

Interpreting the sums in equation (7) as a trace and cyclically permuting, we obtain

(5;) = Tr{T”'2T‘]/Tr[T”]. (10)

We now define the maximum eigenvalue of T and the corresponding eigenvector to
be Am and lm) respectively, so that in the limit of N-> so,

<.5f)=(mlT‘lm)/A§n. 111)
Thus the zero field perpendicular susceptibility is

u(m| T’lm)
z. = .V‘—-———,—.

. (12)
h A?“

The (25 *- 1)2 elements of T‘ may be determined analytically from equation (9) by
examining the 25 ~.- 1)3 states involved. However. since T§;_15;-, is a t'uncrion only of
55f-1 4-3,?.. 1 I, it can be shown that there are just 25 + 1 independent elements of T"r and
hence only (25 + 1)2 states need be examined. A short computer algorithm was devel-
oped to evaluate equation (9) analytically for larger spins. For 5 = l, we find

it p ,8] sh
T 7"“? ‘2‘ (1.,
x = .3
h
_/:h_ 75inh/371

11

Letter to the Editor L239

while when S = 1
h sinh 2p] 2/1 sinh x3]

.,
J J ‘fih
2h ' 2 '
T‘ ___ smh l3] 2 1312 h Slnh l3! ( 14)
I J
7 ' - 1
2&1 ..h 51:11 6] h Sinlli .31

Combining well known transfer matrix results with (12) and (13). F isher’s ( 1963) result
for S = éis recovered:

LL: 1' 3fl+____tanh(fiJ/4)‘l ‘
N112 i[,fisech 4 . j J (13)

while for S = 1 we obtain
it- _ (cosh l3] - l 4- iq) sinh 213] - 85inh (3] + ZBJq

.V‘u‘: !q(cosh p1 + l + sq )2 “6)

 

where

q = [(2 cosh 51 - 1); 4- 3]“. (17)
With Tand T" both known analytically. equation ( 11) may be evaluated numerically for
larger spins. The results for S = i through 5 = 3 are displayed in figure 1. By evaluating
equation (6) considering the two low lying states, all spins pointing up and all spins
pointing up with the exception that Sf =5 - 1, it can be shown that?

1X a (T = 0)

.v.u= =1 (18)

independent of S. as is shown in figure 1.

 

0.7

0.6

 

 

 

 

 

Figure l. The perpendicular susceptibility (IL, .VuJ) pictted against reduced temperature
(kgT/J) forS = i, 1. i, Zandi.

1‘ This result can be seen most easily by using ordinary second-order perturbation theory at T = i) to give
_\E = 1V(uhl/2)2[(VIZS)Z/ZIS] = fix; hl.

12

L240 Letter to the Editor

 

 

 

 

0.1.0 1 L 1 l l l

 

O 0.1 ).2 13.3 DJ. 3.3 1.6 3.7

Figure 2.. The perpendicular susceptibility (Jr -. .Vir') plotted against the reduced tempera-
ture (kgT/[lSlS - l) I). This 15 the same as figure 1 except for the rescaling or the temoer~
attire axes. N0te that the vertical scale does nor begin at zero. Values or 5 are shown for each
curve.

The high-temperature limit

,‘(_ 5(5 1- l)

_, = ___—.. (19
W 13 .3 )
suggests rescaling the temperature. The plot in figure 2 of x_ against kBTx’SlS ~ 1)
indicates that in addition the location of the peak scales roughly with 5(5 -- 1) for
S 2 1.

The perpendicular susceptibility for an infinite spin é chain in the presence of a
non-zero parallel field has already been calculated exactly by Grest and Rajagopal
(1974). By including parallel field terms in the exponents in equations (8) and (9 ). similar
results may be obtained numerically for larger spins. Finally. we note that some modi-
fications after equation (7) will yield the perpendicular susceptibility for finite open Ising
chains and hence for dilute Ising chains (Thomsen and Thorpe 1983).

We are indebted to Dr S D Mahanti for several useful discussions.

References

Fisher M E 1963!. Math. Phys. 4124

Grest G S and Raiagopal 1974 J. Math. Phys. 15 589

Katsura S 1962 Phys. Rev. 127 1508

Suzuki M. Tsuliyama B and Katsura S 1967 J. Math. Phys. 8 124
Thomsen M and Thorpe M F 1983 to be published

13

J. Phys. C; Solid State Phys., 16 (1983) 4191-4198. Printed in Great Britain

Dilute Ising chains with spin Si‘

M Thomsen: and M F Thorpe

Department of Physics and Astronomy. Michigan State University. East Lansing. MI
18824, USA

Received 17 February 1983

Abstract. Randomly diluted one-dimensional Ising chains of spin 5 are studied by breaking
the systems down into non-interacting finite segments. Exact closed-form expressions are
obtained for the specific heat and perpendicularsusceptibility of 5 = Echains and numerically
exact results are given for spins up to it.

1. Introduction

The problem of a dilute one-dimensional spin system with nearest-neighbour interac-
tions is tractable becauseiit can be reduced to a collection of non-interacting finite
segments. Matsubara and Katsura (1973) and Thorpe and Miyazima (1981) have
exploited this fact to calculate thermodynamic properties for dilute XY quantum chains.
Matsubara er al (1973) have investigated the specific heat and parallel susceptibility of
dilute spin-i Ising chains whilst Matsubara and Yoshimura (1973) have investigated the
same quantities for S > i.

In this paper we use this technique to calculate the transverse susceptibility for
spin-S dilute chains for the firsr time. This is an extension of previous work on pure
spin-S chains by Thorpe and Thomsen (1983). The main problem to overcome is that
the transverse field term in the Hamiltonian does not commute with the exchange part
of the Hamiltonian. This can be handled as shown previously (Thorpe and Thomsen
1983) by introducing a new transfer matrix T". In addition we have to introduce another
transfer matrix T“ in the dilute case to handle the situation when the field is on a spin at
the end of a segment. Details of this calculation and results for 5 = i, 3 and it for
concentrations c— = 0. i. %. fland 1 are given in § 3.

In § 2 we calculate the specific heat using the transfer matrix method. This serves to
introduce much of the notation although results have previously been obtained for the
specific heat and parallel susceptibility for S > i by Matsubara and Yoshimura (1973).
In neither of these calculations does one have to worry about non-commuting variables.

For 5 = i it is well known that the specific heat per bond is independent of c. This is
because each bond acts independently. We show that the specific heat per bond is weakly
dependent on concentration for 5 2 1.

* Work supported in part by the National Science Foundation.
2: Exxon Fellow.

© 1983 The Institute of Physics 4191

14
4192 M Thomsen and M F Thorpe
2. Specific heat

The transfer matrix technique was introduced by Kramers and Wannier (1941) to solve
infinite and finite closed chains in the Ising model. One starts with an Ising Hamiltonian
with a field in the z direction:

=-leSi‘55-1-,uh 215.; (1)
where if S..1 a S the chain is closed. The partition function is
Q: = 1315le = Tr (D; exp{fl[‘,5izsiz*l " iuhle *' 55-1)”) (2
If one defines a (25 + 1) X (25 + 1) transfer matrix Tby
73:5!“ = eXPifilJSi’Si’q '1' iii/“Si: ‘1’ 55.1)” 13)
where the 5,? refer to eigenvalues. nor operators. then ('2) becomes
Q: = 2 E . . . Z (SilT’5§f><S§iTTS§>- . . cars.» .4)
5g 5g 5%,
= Tr[T']. (5‘)
Taking it,- to be the jth eigenvalue of T.
1.5.1
Q: = 2. A;- <6)
,.

To modify this for an open chain. we note that we need only set the bond energy of
the 5fo pair to zero. In analogy to the derivation of (5) we find for the open chain

Q? = Tr[T" l7"] (7)
with
T§fsf = T55;(J = 1)) '3 CXPHBUJHSE "’ $13)]. {8‘

The partition function ('7) is generally evaluated by rorating T and T' into the represen-
tation which diagonalises T.

In this section we wish to find the zero field specific heat of a finite open segment.
NOting that in zero field, T533 =1 for all 55 and S 1’ , and defining M as the matrix that
diagonalises T. it is easy to show that

ZSrI 234-1 3
o... ' ‘-li ‘
Qr - S A: '2 -Will- ‘9)

1']. i [‘1
For 5 = i, this yields
Q? = 2(2 cosh ilil)" ‘. ('10)

while for S = 1 we find

3-'" .- it‘-3 .-
Q9=ZC _A._lA)'*7'f:-T:(A )' (11)
A. A. A.

 

15

Dilute Ising chains with spin 5 4193
where
= = coshfiJ i-l:[(cosh13’J—l)2 + 2]“? (12)
Using
1 _ a 2 a )
hie-”(r H In Q. (13)

it follows that in the S = 1 case
Cr/ke = (r - 1) (Bl/4)z sechzffll/4) (14)

(see also Stanley (1971)), while application of (13) to (11) yields an expression too
complex to be of much use analytically. For general spin. the matrix calculations involved
in (9) are performed by computer at three points close to each other so that the
differentiation in (13) can be done numerically.

4 1 41 4,4 714
//’°/°/// ° °//// °
W —-I M W
r=2 r=l r=3 r=4

Figure l. A dilute spin-S chain with segments of length r = 1. Z, 3 and 4.

In a dilute chain of length N (N -> an) and concentration c, the number of seg-
ments‘containing exactly r consecutive spins (as shown in figure 1) is

P, = N(1 - C)ZC’. (15)
Since
C=2ec do

we find that the specific heat per bond in the spin-l case is
C/Nkac‘ = (191/4) seen: (131/4). (17)

which agrees with Matsubara er al (1973).

We note that the specific heat per bond is independent of the concentration c because
the spin-i Ising chain can be mapped onto an equivalent system of independent bonds.
This is done by defining a new variable r,- = 5,5..1. Because of the absence of loops the
t.- can take on values 1'.- = t i and so in the absence of a field h, the Hamiltonian ( 1) splits
up into non-interacting Hamiltonians each associated with a bond with a spin at each
end.

For 5 > i, the specific heat per bond is notindependent of c, as can be seen in figure 2.
However, the specific heat per band depends only weakly on the concentration and this
dependence does not become greater as 5 increases.

A useful check is provided by noting that

(”gnaw-5(a) <18)
6

where 5( 00) is the entrOpy of the system at infinite temperature and 5(6) is the entropy

16

 

 

 

 

 

 

 

 

 

 

 

   

 

 

 

 

4194 M Thomsen and M F Thorpe
10+ 10 .
-. r,
S a .55." ‘K S a 1
i If: a"; “1%
05 ) 05 3,),- “K
B n I
..“ o 10 20 o 10 20 so
i
T 7-01 - 1 2a - -
1.
15> 5:: . 1_g- 5::
-. z : u. 2
~15
:5: 1:15.
10» If! l
as
1‘ ix
\
‘ A A —__.Jl
1.0 7.0 310 o 1.0 7.0 so
kJ/JS’

Figure 2. The specific heat per bond (C/Nc‘k.) as a function of reduced temperature
(kg 1713’) for various values of spin 5 as indicated and concentration c - 0 ( ----- ), l
(- - - - -),§(-—--),i(-—--), 1( ). ForS - lthespecificheatperbondisindependent
ofcasshown. Thee=011mitcorrespondstotwospinsjoinedbyonebondandotherwise
isolated.

 

at a small positive temperature 6. Since the total number of states is
(2.9 + 1) m“*‘" = (23 + 1)”‘
S( on) = Nkac ln(ZS + 1). (19)

On the other hand, the degeneracy of the ground state (in zero field) is 2”“(25 +1)”1
where N1 is the number of isolated spins and N2 is the number of segments of length
greater than one. This gives

3(6) = Nk3(1 - c)c2 1n 2 + Nk3(1 - c)=c 1n(2.s + 1) _ (20)

and hence
[:ng = Nkaczflz - c) ln(2.S' + l) - (1 - c) In 2]. (21)

The low-temperature behaviour of the specific heat is dominated by the lowest
excited states. For S = 1, this corresponds to flipping all the spins simultaneously in any
part of that segment, provided that part includes an end spin, as shown in figure 3. For
S 2 1, the low-lying excitations come only from reducing by one unit the 2 component
of an end spin which is also shown in figure 31. It can thus be shown that at low

1‘ For spin 1 , additional degeneracy arises for the lowest excited states in segments of length 2. However, it is
still true that all the lowest-lying excitations are governed by end segment spins for S - 1.

17

Dilute [sing chains with spin S 4195
temperatures.

C/ngcz == ([315): e“"’5 S = .1 (22a)

C/ngc2 z [3 - lc](1 - c)(,BJ.S)2 e‘w S‘= 1 (22b)

C/ngcz 2' 2(1 - c)(l3JS)2 e‘ws S 2 '2‘ 22C)

llloolllllllloo
ll°°\llllllll°o

Figure .3. The lowesr excitations for 5 = 5 and S 2 1 dilute chains as described in the text.

When c-> 1, the right-hand sides of equations (22b) and 225) vanish since end-of-
segment spin excitations are no longer possible. which can be seen in figure 2 where the
c = 1 specific heats go to zero much more rapidly as T-> 0 for S 2 1. This is because
there is an extra factor 2 in the exponent of the next terms in (22b) and (2c) reflecting
the fact that every spin has two neighbours rather than one.

We conclude this section by noting that when one evaluates the high-temperature
expansion of Q = Trfe ”3”] in powers of )6. it is found that

C/Nk3c3 2 [15(5 + 1))3112 (23)

which is independent of c for all S.

3. Perpendicular susceptibility

Calculation of the perpendicular susceptibility is complicated by the fact that parts of
the Hamiltonian do not commute with each Other. In the absence of any parallel field

r-ol r

H=-J,§.‘lstss-.-,Lm‘ZStaHo—H, 124)
ta :-1

for an open segment of length r.

Using linear response theory. Grest and Rajagopal ( 1974) calculated the perpen-
dicular susceptibility for a closed spin-Ii chain. With transfer matrices, this was gener-
alised to higher spins (Thorpe and Thomsen 1983). It was shown that the contribution
of each spin to the magnetisation is

r I -5H 1‘s
X 7.51-351-) T span. 7.55.1554 . - - Tsist/Tl’ie ‘ ”l 14-5)

18

4196 M Thomsen and M F Thorpe

where T53),l is as defined in (3) and
3
first. = .qu 2f [0 expmr - r)(57-15,~‘ + s;s;..>is;

x expws;., s; + s; 51.1)] dTSf. (26)

UsingSf 1(5 * + S,‘ ), the matrix T’ can be evaluated analytically for any 5. In par-
ticular, when S= i

((W/J) sinh llil With"

“3“”! (WI/051““ “31 .
We now open the chain by setting 1 = O for the SSS? term. If S,- is nor at the end of a

segment, then (25) becomes

(SX)=SE§§.. 2 2 ”@7333 T593...

3-1 Sin

(27)

X 7.57.2.3“ Tx5f_l3f.l Tng 55.2 . . . T3145; T'gsg/TI'R -dHo] (28)
where T335? is defined in equation (8). In the case of an end spin,

(51) = (55) = 35.; 2:. - - - $2 Tstngsgsg- - - 755-33-. Trams/[THC "W1 .29)
1 2 ’1
with the new transfer matrix T" defined as
p
11's.“.1 = 52 we I exp[(fl - r)Js*-ls:]5: exp[fiS‘_( 5:] dist. (30)
0

We n0te that the subscript Sf appears for formal reasons only—7‘55” 5g is. independent
of Si. For 5 = l, we obtain

= (OW/J) Sinh W (2147!) sinh is!)
(ZW/Dsinhffl] (zmx/Dsinhffi] '

Using the above equations. we find the zero-field perpendicular susceptibility of a finite
open Ising chain to be

(31)

x”: ;[(2'1'1'[ richer-1' 1T1) +2TrT"’T’]/TIT' ‘7‘]. (32)

This expression is valid for r 2 2. For r = 1, one uses the free-spin Curie susceptibility.
With all four matrices in equation (32) known analytically, it can be shown that for
S = 1

X1! = [(r - 2)/2](/3/4) Sit—ch2 $131 + [(r + Zl/Z/l tanh H31- (33)

Performing a sum similar to equation (16). but singling out the special case of r = 1, we
find that for a dilute spin-l Ising chain

Jib/Neill = (1 - 021137 + 1181 c2 sech‘llil + 1(4c - 3c?) tanhltil. (34)
Setting c = 1 reproduces Fisher’s (1963) result for a spin-l pure Ising chain

JxJ/Nis2 = 1,3! sech2 llil + l tanh 1131. (35)

19

Dilute [sing chains with spin S 4197

As c approaches 0 only the first term in equation (34) survives so that
JxJNuzc e W

as expected for No free spins.

All four matrices in equation (32) are known analytically for larger 5. but the final
evaluation of that expression must in general be performed numerically. The results are
displayed in figure 4 where we have plotted the inverse susceptibility per spin as a

 

 

 

 

”Xi/Hip“.|

 

 

‘2

 

 

 

 

C.)
2
,)

1 £0 .10 i 2.0 1
k3 I'USJ

Figure 4. The inverse perpendicular susceptibility per spin [ILA/ctr)“ as a funcn’on of

reduced temperature {kg 1715"] for various values of spin 5 as indicated and concentranon
C’0( """" Li‘l---'°).i('~'-).i(---).ll l.Thec=0limitcorrespondsroa

single isolated spin.

 

function of temperature for four different spins: S = is, 1. 5%, 3. We note that 1:1
approaches zero linearly for c < l, reflecting the effect of .Vcil - c)z isolated spins
[X1 = iNuzd 1 - c)25(S «'- 1)/k3 T]. For c = 1, we recover the result for the pure chain
that x; = lN/f/J. At high temperatures we obtain the Curie susceptibility of NC isolated
spins [1 g = lNuch(S 4- l)/k3 T] so that [J x J’Nc‘tfl‘l becomes independent of c.

4. Conclusion

The perpendicular susceptibility per spin and the specific heat per bond do nor change
dramatically upon dilution over much of the temperature range studied. The leading

20

4198 M Thomsen and M F Thorpe

term in the high-temperature expansion of each of these quantities is independent of
concentration. For the special case of spin l, the specific heat per bond is independent
of concentration for all temperatures. The most n0ticeable effect of diluting an Ising
chain is in the way these thermodynamic quantities behave in the zero-temperature
limit. The specific heat per bond differs in the nature of its exponential approach to zero
upon dilution (for S 2 1). Most striking is the divergence of X- at T = 0 when the chain
is diluted, owing to the finite probability of finding isolated free spins.

We note that the results for both the specific heat and the perpendicular susceptibility
depend only on Ill and so are identical for an antiferromagnetic Ising chain. This is most
easily seen by performing a rotation of 180° about the x axis.

One can expand on this study by turning on a weak interaction between segments,
particularly those separated by just one non-magnetic site. However. the
independent-segment assumption and hence equation ( 16) are not valid in this model
and the problem becomes significantly more complex.

Acknowledgment

The authors wish to thank S D Mahanti and T A Kaplan for several valuable discussions.

References

Fisher M E 1963 J. Math. Phys. 4 124

Grest G S and Rajagopal A K 1974 J. Math. Phys. 15 589

Kramers-H A and Wannier G H 1941 Phys. Rev. 60 252

Matsubara F and Katsura S 1973 Prog. Theor. Phys. 49 367

Matsubara F and Yoshimura K 1973 Prog. Theor. Phys. 50 1824

Matsubara F. Yoshimura K and Katsura S 1973 Can. J. Phys. 51 1053

Stanley H E 1971 lnrroducrion to Phase Transitions and Critical Phenomena (Oxford: Oxford University
Press)

Thorpe M F and Miyazima S 1981 Phys. Rev. B 24 6686

Thorpe M F and Thomsen M 1983 J. Phys. C. Sou'd State Phys. 16 2237

21

Much of the preceding analysis can be carried out in the
presence of an external parallel field, hz. For instance, using
equation (4) the standard transfer matrix for spin S is determined
in a parallel field. From this one can obtain the partition function,
although for S>%, the transfer matrix is not easily diagonalized
analytically.

Similarly, one can write down the elements of the TX transfer
matrices, introduced in the X1 calculation, in the presence of

a parallel field:

 

 

 

 

 

 

r3
X X 3-. .- _ X “X
T32 Sz - h g2 JO expliB r)T1] Si exp(rfl) 51 dr (17)
i—l 1+1 1
TX, - n‘ z (a x [(B-t)F 1 3X ex (If ) sX dT (18)
2 S2 C Sz 0 e p 2 r .p 2 r
r—1 1 r
where
Z Z Z Z Z Z Z
F1 — JSi(Si_l + Si+l) + h [(Si-l + Si+l)/2 + Si] (19)
Z Z Z Z Z Z
r2 = JSr_1Sr + h [(sr_1 + 51)/2 + sr] . (20)
For instance, for S=%,
exp(BhZ/2)sinh[%g(hz+d)] sinh(BhZ/2)
J+hz h2
.l; TX = (21)
h sinh(BhZ/2) exp(-th/2)sinh[8(J-hz)/2]
hz J-hz
and
2exp(3hz/2)sinh[8(%J+%hz)] ZSinhLB(%J+%hz)1
J+2hz J+2hz
l__Tx'= (22)
hz 23inh[8(%J+%hZ)] 2exp(-th/2)sinh[8(%J-%hz)]

 

 

J+2hz J-th

22

To calculate the transverse susceptibility for a pure chain,
one needs to calculate the maximum eigenvalue, Am, and corresponding
. x . .
eigenvector lm> of T . Then by equation (12) of the preceding
letter, "Transverse susceptibility in spin S Ising chains," one

can obtain

Xl exp(th/2)sinh[B(J+hz)/2] + exp(-th/2)sinh[B(J-hz)/2l

 

 

 

 

= Z Z
Nui 2(J+h ) 2(J-h )
+ eprBJ/Ajsinh(BhZ/2) eprBhZ/2)sinh[B(J+hz)/2]
2q1 J+hz
_ expg-BhZ/2)sinh[B(J-hz)/21 + 2expg-BJ/2) ]
J-hz hz J
X [exp(BJ/4)cosn(BhZ/2) + qll-Z (23)
where
q1 = [exp(BJ/2)cosh2(BhZ/2) - 23inh(BJ/2)] % . (24)

This result was obtained earlier by Crest and Rajagopal (1974).

The modified transfer matrix, TX, can be obtained straight-
forwardly for higher spins in a parallel field, but analytic cal-
culation of Xl is no longer tractable. Nevertheless, numerical
results may be obtained easily.

Inspection of equation (23) reveals that for hz=iJ, there is
the possibility of a divergence due to vanishing denominators. In
fact, an antiferromagnetic chain in this field will have a divergent
zero temperature transverse susceptibility, as noted by Crest and
Rajagopal (1974). This can be seen for general spin by considering

the Hamiltonian associated with a single bond:

23

22 Z Z 2
H12 = Josls2 - h (31+32)/2 (25)
with
JO 2 —J > 0 . (26)
Taking hz = 2JOS,
Z Z Z
H12 - (Jos2 - J05)s1 - Joss2 . (27)

H12 can be minimized with respect to S? by setting S§=S, yielding
a ground state energy
E12 = -JOS“ , (28)
independent of 8:.
As long as one spin in every bond is parallel to the field,
the system is in a ground state. This state leaves a large number
of spins free to point in any direction, which can be seen by noting

that E12 is independent of S These free spins are then able to

z
2.
line up in the direction of an infinitessimal transverse field,
hence creating the divergence.

In particular, from equation (23), one can show that at H=iJ

and as T + O,

xi/Nui = s(1+/§)/(10+6/§) . (29)
The dilute chain transverse susceptibility in the presence
of a parallel field may also be obtained using the formalism developed
in the preceding paper. It should be noted that now the contribu—
tion of each spin in an open segment is a function of its location
in the segment. This is a direct consequence of the fact that
T BHd T, are not simultaneously diagonaliéable in the presence of

a field. Thus for dilute chains, analytic work stops at the writing

24

? x X,

down of the four matrices T, T , T , T , and the matrix multiplica-
tion and diagonalization is performed numerically.

It was seen earlier that the isolated spins contributed a
divergence at T=O to Xl' Application of a finite parallel field
destroys this divergence. However, for the same reasons outlined
for the pure chains, an antiferromagnetic chain will have another
T=O divergence in Xl when hz=32JS. A third divergence occurs at
hz=.JS due to the spins at the end of open segments. Figure 1 il-
lustrates one of these divergences, at hz=JS, with a plot of xll
versus temperature. Since the divergence is proportional to 3,
x11 approaches zero linearly, as shown.

One limitation of this formalism that has been developed for
the calculation of XI is that it is based on linear response theory
and hence is valid only for infinitessimal hx. However, in exchange
for this restriction, one has gained the ability to calculate Xl

for a variety of models that are not tractable by other methods.

3.5

2.5

iJIXL-l 2

Ncu2

 

0.5

Figure l.

25

 

 

 

 

DIOT
.’ .4
d ."’ '.".
a, -.'
0’ 0..
0’ "
a’ 0'
'I '0‘
d I./ x”
O" "
a’.’ ..o"
J “-—'-"'—. .0 """
'l"- oooo '-
-4 ,1"- “‘----------"".."
‘ ‘—‘H/Ulfl15
_/’ ----. r1/ 1,: 3:0 .4
- --- HI (.1 [20.6
I I I I I
0 0.5 l 1.5 2 2.5
kBT/IJISZ

The inverse of the transverse susceptibility
per spin of a dilute (c=O.25) Ising chain in
the presence of three different parallel fields.

3. P(h): The Local Field Probability Distribution Function

The local field probability distribution function, P(h), has
been studied at zero temperature for a variety of spin glass models
(see, e.g., Palmer and Pond 1979 or Walker and Walstedt 1980).

One of the distinctive features of P(h) these studies bring out
is a dip at h=O. The nature of P(h), measuring the distribution
of effective fields acting on a single spin, suggests that it be
relegated to mean field theory or to low temperature studies. It
is probably for this reason that the temperature dependence of
the distribution has not been studied much before this. It turns
out, however, that for Ising systems one may obtain exact thermo-
dynamic relations using P(h).

This chapter is devoted to a formal discussion of the properties
of P(h). While Appendix A contains details for a very broad class
of Ising systems, for simplicity a more restrictive Hamiltonian
is used here:

H = —% Z J..o.o. (1)
i,j 13 l J

where the 01's (=il) are spin % Ising variables. The matrix de-

fined by the Jij's is symmetric with vanishing diagonal elements
but is otherwise arbitrary.

One begins by defining a local field operator, hi:
h- = Z J..o,
1 j 13 J . (2)

The missing factor of % in equation (2) ensures that hioi contains

26

27

all the Oi terms in the original Hamiltonian. The local field
probability distribution function at site i is defined by

Pi(h) = <5(h-hi)> . (3)

The notation <...> refers to a thermal average. The site averaged

distribution for a system of N spins is

P(h) = tPi(h) . (4)

l

L

Alr—l

This function is of course normalized:
f P(h) dh i 1 . (5)

Associated with P(h) are the following three exact results:

M = N f tanh(Bh) P(h) dh ' (6)
E = ?§ I h tanh(Bh) P(h) db (7)
S(k,w) = %[P(w/2) + P(-w/2)]/[1 + exp(-Bw)] . (8)

These last three equations have simple physical interpretations.
The magnetization of an isolated spin in a field h is tanh(Bh).
Equation (6) performs a weighted average with respect to h. Since
tanh(Bh) is odd, M couples to the antisymmetric part of P(h):
Pa(h). The second equation has a similar interpretation with an
additional factor of % present to correct for double counting of
bonds. E couples to Ps(h), the symmetric part of P(h). Finally,
a beam of neutrons polarized in a direction perpendicular to the
z axis interacts with the system by flipping individual spins at
an energy cost proportional to the local field. Since these neutrons
cannot distinguish between +2 and -z, they couple to Ps(h).

The statement these equations make is significant. Without

any explicit knowledge of the Jij's, one.can obtain the thermodynamics

28

of a system directly from P(h). Furthermore, one can measure Ps(h)
directly, and that is sufficient to determine E.

Detailed proofs of equations (6) - (8) appear in Appendix A.
However, to give the reader a flavor for the proofs, the energy
equation, (7), is derived here. One begins with the energy of an
isolated spin in a field h:

€(h) = —h tanh(Bh) . (9)

In its more primitive form, this is
€(h) = TrOIeXp(-BHO(h)) H0(h)I/TrOIeXP(-BHO(h))I (10)

where

H0(h) = -hoO

and TrOmeans trace over 00 = :1. Denoting the energy of all the
bonds terminating at site i by 2Ei (E1 is then the energy associated

with site i), it needs to be shown that

f €(h) P(h) db = 21?.1 . (11)

The following definitions need to be introduced:

Trl = trace over Oi = t1
Tr' = trace over all spins except Oi
H. = -h.O.
1 1 1
H' = H — Hi’ which is independent of 01
Z = Tr[exp(-BH)]

Then,

I P(h) €(h) dh

 

I h Tr[exp(_BH)5(h-hi)]Tr0[exp(-BHO(h)) H0(h)] (12)
d

Z Tr0[exp(-BHO(h))I

29

exp(—BH) Tr0[exp(-BHO(hi)) HO(hi)]

 

 

1 (13)

= -Tr
Z Troiexp(-BHO(hi>)1

=-l-Tr' exp(-BH') Tri[exp(-BHi)] Tr0[exP(-BHO(hi))HO(hi)] (14)
z

Tr0[exp(-BHO(hi))I

In the last step, it is crucial that [H',Hi] = 0 so that the expo—

nential can be split up. Since
Tri[exp<-eHi>1 = Trotexp<-BHO<hi>>1 , (15)

one can cancel a term from the denominator and the numerator. Now

00 can be relabelled Oi and the traces can be regrouped, giving

I P(h) €(h) dh %-Tr[exp(-BH) Hi] (16)

2Ei , (17)
which was to be shown.

P(h) is fundamentally a lggal quantity coupling to local pro-
~perties such as Bi and Hi. The specific heat and other quantities
which involve long range correlations require differentiation of
P(h). The free energy, meanwhile, is obtained by integration of
E in equation (7). Thus P(h) contains not only all the information
found in the free energy, but it also contains additional information
such as S(k,m).

This formalism does not of course provide short cuts to solving
for the thermodynamics of a system. The price one pays for the
additional information in P(h) is that it in general is harder to

calculate than the free energy. The remainder of this chapter will

be devoted to methods for calculating P(h).

30

This part of the discussion will be limited to ferromagnetic
systems with nearest neighbor interactions (although some of the
techniques will be generalizable). Setting the interaction parameter

J = l, the Hamiltonian becomes

H = - X 0.0. (18)
<ij> 1 3

where X indicates a sum over nearest neighbor pairs and the
<ij>

spins are assumed to lie on a lattice. One of the advantages of

this simplification is that P1(h) = P(h).

The first approach to calculating the distribution is a formal

one starting from the definition,

Pi(h) = <5(h—hi)> . (19)

An integral representation of the 6 function is used, and a little

algebra yields the form

2
Pi(h) = E- wS 6(h—s) . (20)

s z
z is the lattice coordination and the sum is over the allowed

values of the local field, -z, -z+2, ..., 2—2, 2. Details of

this approach are in Appendix A. It is shown there that the weights
ws involve numerical factors, which depend only on coordination
number (not on the lattice itself), and all possible correlations
among the nearest neighbors of site 1. Hence the real work involved
with this calculation becomes the determination of those correlation
functions. This method will be demonstrated in Chapter 5.

For a system of coordination of two or three the correlations

on the nearest neighbor shell can be eliminated in favor of the

31

nearest neighbor correlation, <OOOl> and the magnetization, <OO>.
If, for example, z=3, then the possible field values are :3 and :1
and so four equations are needed to determine the four weights, ws.
These equations are obtained by using the normalization requirement,

(5), the energy and magnetization relations, (7) and (6), and the

fact that <h> = 3<O >. One then obtains

0
(W3 + w-3) + (w1 + w_1) = 1 (21)
(w3 + w_3)3tanh(3B) + (w1 + w_1)tanh(B) = 3<oOol> (22)
(w3 - w_3)tanh(38) + (w1 - w_1)tanh(8) = <oO> (23)
3(w3 — w_3) + 1(w1 - w_l) = 3<oO> . (24)

The first two equations yield Ps(h) and the last two, Pa(h).

These equations can be solved explicitly to yield

 

 

 

 

1 3(tanh(3B) - (0001)) 3<OO>(tanh(3B) - 1)
w = -- i (25)
:1 2 3tanh(3B) - tanh(B) tanh(38) - 3tanh(B)
1 3<0001> - tanh(B) <OO>(l-3tanh(B))
w = -' i" (26)
:3 2 _3tanh(38) - tanh(B) tanh(3B) - 3tanh(B) ’

 

showing explicitly in this case that the w's depend only on (0001)
and <OO>. No reference is made here to the particular network
other than it has coordination three and is compatible with a
Hamiltonian of the form given in (18).

Finally, some systems lend themselves to calculation of the

moments <hn> by standard thermodynamic methods. If one knows all

the moments up to <hz>, then the 2+1 equations,

<hn> = Z w s (27)

32

can be inverted to find the ws directly.

These methods may be generalized to sgmg other systems, such
as antiferromagnets, but their generalization to disordered systems
is nontrivial. Different and specialized techniques must be de-
veloped for calculation of P(h) on these systems and so discussion
of such techniques will be postponed until a specific disordered

system is introduced (Chapter 7).

4. P(h) for Bethe Lattices

A simple system for which P(h) can be determined consists
of S=% Ising spins on a Bethe lattice with constant nearest neigh-
bor ferromagnetic interactions. The Bethe "lattice" is really
a pseudolattice (see Figure 2). It has the unfortunate property
that, even in the thermodynamic limit, a finite fraction of its
sites lie on the surface. On the other hand, it has the advan-
tage that its coordination number, 2, is easily varied. The z=2
case, a linear chain, will not be considered here since many of
its properties are distinctly different from the z>2 Bethe lattices.
Because of the large surface area, calculations must be ap-
proached with caution. For instance, Eggarter (1973) pointed out
that by assuming free boundaries, the Bethe lattice does not
exhibit a phase transition. However, the more traditional Bethe-
Peierls approach (Bethe 1935, Peierls 1936) which assumes boundary
conditions determined by self-consistent fields, does find a phase
transition. This latter approach will be used in this discussion.
The Bethe-Peierls method, leading to the self-consistent
field equation (8) below, is described in Domb (1960) and is
shown to be an exact solution for the system. An effective Hamil-
tonian is assumed for a cluster consisting of site 0 and its
2 nearest neighbors:

Z Z
H = —J 2 o o. — H o - H 2 o. . (1)
z j=1 0 3 O O lj=1 J

33

34

> <2
2 é ’

Figure 2. Section of a Bethe lattice for z=3(upper)
and z=4(lower).

 

 

35

H0 is an external field acting on the system and will be set to
zero later in the calculation. H1 is an effective, temperature
dependent field acting on the Oj's accounting for both HO and
the z-l neighbors (other than 00) of Oj. It can be shown that

the partition function for H is

z
-%,—-1—-% 222 22-. 2122
Q2 = A0 ( 1 x + Alx ) + AO (XIX + Al x ) (2)
where
10 = exp(-2H0/kBT) (3)
11 = exp(—2H1/kBT) (4)
x = exp(-2J/kBT) . (5)

To obtain H1, one requires

<0 > = <0.)
0 J ’

or equivalently,

A
3 1 8
A0 575 £an = z— Efi £qu - (6)

Using equation (2), some algebra yields

11/10 = [(11 + x)/(1 + 111012“1 . (7)

The external field HO will now be set to zero, having served its
purpose. Some rearrangement gives

tanh[H1/(z-1)kBT]
tanh(H1/kBT)

 

= tanh(J/kBT) , (8)

and this specifies H1 as a function of temperature. An expansion

for small H1 gives a transition temperature defined by

tanh(J/kBTC) = 1/(2-1) . (9)

36

Above Tc’ H = O, and the effective partition function becomes

1

Q2 = 2[cosh(J/kBT)]z , (10)

corresponding to 2 independent bonds. This result is typical for
spin % Ising systems with no loops, in the paramagnetic phase
(see, e.g., Thorpe 1982). It should be noted that Qz determines
only local quantities such as the energy associated with 00 or

the local magnetization <OO>. It is 22; true that the free energy

F is given by -BF = Zan.

Below TC, one needs to determine H1. Introducing two new

variables

Y = All/(2‘1) (11)

rt
ll

tanh(J/kBT) , (12)

it is now possible to rewrite the self-consistent equation (8) as

-1
l-Y 1+1z

—'—:-=t <13)

1+Y 1_Yz 1

01'

—1 2'2 i

(l-t)Yz - 2t 2 Y + (1-c) = o . (14)
i=1

By symmetry (H1 ++ -H1), one expects that if YO is a solution to
this equation, Yal is also a solution. This fact can be exploited
to reduce the degree of the polynomial by a factor of two. Since
a fourth degree polynomial is solvable (Burington 1946), one can

in principle obtain analytic expressions for H1 up to z = 10,

In practice, the complexity of these solutions makes some of the

37

larger 2 results not worth pursuing analytically.
For completeness sake, analytic expressions giving X1
(=exp(-2H1/kBT)) as a function of x (=exp(-2J/kBT)) are shown

below for z=3 up to 2:6:

1
z=3 A1 = %[x - 2x_1 — 1 + (xFZ-x_1)(l-2x-3x2)2] (15)

- .102] <16)

2: A = %[x

1

1
2:5 Define C = [1 - x + (5x2+2x+l)2]/2x

Then 11 = {[c + (c2-4fi1/2}4 (17)
2:6 Define D = [x_1 + (x-2+4)%]/2
Then A1 = {[D + (DZ-4)%]/2}5 . (18)

A system with infinite range interactions produces a mean
field result (Stanley 1971). On the Bethe lattice this corresponds
to the z +'m limit, which is readily accessible. In order to
keep the energy finite, this limit must be taken holding Jz con-

stant. Since Hl should sCale with Jz,

l L
xii xi = 1 2 , (19)

NIH
,__. 1+

reducing the partition function (2) to
2
Q2 = 2cosh(HO/kBT)[2cosh(H1/kBT)] . (20)

The self-consistent field equation, (8), becomes

Hl/Jz = tanh(H1/kBT) . (21)

Which has a nontrivial solution if and only if kBT < Jz. Hence

kBTC = Jz . (22)

With the fundamental equations in place, it is now possible

to calculate P(h) for Bethe lattices. This is most easily done

38

by calculating the moments of the distribution. Inspection of
the Hamiltonian (1) reveals that these moments may be obtained

by differentiation of Q2. Taking
h/J = Z 0. , ' (23)
- J
J
(assuming zero external field), one finds

a n
<(h/J)n> = —' [-ZA T] Q . (24)

D—4
Q)

As pointed out in the previous chapter, given the moments
up to <(h/J)z>, one can solve for the weights in the 2+1 6 functions
associated with P(h). This process, while straightforward, yields
very complex analytic expressions which by themselves do not provide
much insight. For this reason, they are omitted here in favor
of the numerical results below.

The mean field result is easily derived. iStarting from the
partition function, (20), and applying a modified form of equation

(24), one obtains for H = O,

 

 

0
<(h/J)n> = (2cosh8H1)‘z [SBHI n(2coshBH1)z (25)
= z(z-1)...(z-n-1)tanh“BH1 + 0(zn‘1) (26)
Thus
<(h/Jz)“> = tanhnBHl + 0(1/2) (27)

This implies the 6 function distribution

P(h) = 5(h - JztanhBHl) . (28)

Using equation (21), this can be rewritten

P(h) = 6(h — H1) , (29)

39

which is another way of saying that in the large 2 limit, mean
field theory becomes exact.

In Figure 3, P(h) is plotted against h/z for z = 3, 4, 6, and
infinity. It is symmetric for T 2 Tc’ reflecting the fact that
the magnetization is zero. At infinite temperature, the weights
are determined by counting the number of ways to make a particular
value of h. Since there are more ways to make h = O, the distri-
bution is peaked there. When T = TC , the distribution is still
peaked at small h, with more weight being at h = O as the mean
field limit is approached. The symmetry is broken below TC when
spontaneous magnetization sets in. Again the approach to mean
field theory is seen as the maximum in P(h) moves from h = Jz (z = 3)
to h = sz for larger 2.

It is apparent that there is strong dependence of P(h) on
coordination number for Bethe lattices. While the number of 5
functions is clearly determined by z, the profile of the distri-
bution, particularly at Tc’ is also sensitive to 2. It is interesting
to note that the nearest neighbor correlation at TC is also strongly
2 dependent:

<oooj>T=TC = 1/(2—1) . (30)

It will be argued later that this correlation plays an important

role in the profile of P(h).

4O

 

 

 

.8 U53

’3

 

 

 

 

 

 

 

6
8n H\B

N.

o
An E\H

.0

.q .m :ofiumcfivuooo mo mmoflsumfi onumm com ALVA .m ouzwflm

N\:

 

 

 

 

 

 

 

 

 

 

 

 

 

o
e.gn %\E

o.ouue\e

 

o
0.0" E\H

 

1 «.3

A 3.3

«.9

...,.

. 3.:

sum

N

\0
ll

Aces

5. P(h) for Two Dimensional Ferromagnets

There are several two dimensional Ising systems for which
P(h) can be calculated exactly. These analytic results can provide
useful insight into systems with phase transitions not dominated
by surface effects (in contrast to the Bethe lattice).

The first two dimensional Ising system to be solved was the
square net with uniform nearest neighbor interactions (Onsager
1944). Taking the ferromagnetic case and setting J = 1 yields
a Hamiltonian of the form

H=~Z 0.0. . (1)
<ij> l 3

Since the square net has coordination four, the possible values

of h are :4, t2, and O. The weights of the corresponding 5 func-
tions in P(h) are easily expressed in terms of correlation functions.
using equations (38) and (43) and Table 1 of Appendix A. The
nearest neighbors of an arbitrary site, 0, have been numbered

cyclically. One finds:

”:4 = “1'6 1 I“? + 2150102> + ém‘lcb>
i-%<010203> + IE<01°20304> (2)
“:2 = T * 2‘01) 1 %<"1"203> ' 2‘0102’304) ° (3)
W0 = g ‘ 23102> " F010? + %<O102‘53°4> ' (4)

Apparently. <010203> is not in the literature, nor is it easily

calculable. This problem is dealt with by using the magnetization

41

42

equation ( (6) in Chapter 3) to eliminate <010203> in favor of
<00),

The even spin correlations can be reduced to elliptic integrals
using the Pfaffian method described by Montroll, Potts and Ward
(1963). While all of the two spin correlations needed appear ex-
plicitly in that paper, the four spin correlation does not. Details
concerning the calculation of the latter in terms of elliptic
integrals are given in Appendix B. The final step, evaluation
of the integrals, must be performed numerically.

In Figure 4, excerpted from the preprint in Appendix A, P(h)
is shown at five different temperatures. The infinite temperature
distribution, as for the Bethe lattices, has roughly a Gaussian
profile reflecting the fact that there are more h a 0 states than
h = 14 states. However, in contrast to the Bethe lattice, a pro-
nounced dip has developed at Tc' This dip will be discussed in
more detail later in this chapter. Below Tc’ the asymmetry in
P(h) develops very quickly, reflecting the rapid rise of the
magnetization.

Results for the honeycomb lattice are more easily obtained.
Since the lattice sites are threefold coordinated, equations (25)
and (26) of Chapter 3 apply and one needs just the magnetization,
<OO>, and the nearest neighbor correlation, <0001>. The former
has been calculated by Naya (1954), while the latter may be obtained
by differentiation of the partition function in Green and Hurst

(1964). The expression for <0001> thus obtained involves a double

integral, making it somewhat awkward to evaluate numerically.

43

 

 

 

 

 

 

 

13/1}:=(313
0.5-
c) ,_
-F/-EC= C)f3
0.5
O 1 A
P(h) T/TC=I.O
0.5-
V
C) E3 E3 E3 El a;
T/Tc= 1.2
0.5‘
0 12
T/Tc=00
0.5-.
‘4 '2 O 2 4
h
Figure 4. P(h) for the ferromagnetic square net.

 

44

Analytic work reduced this form to a single integral, details of
which appear in Appendix B. A plot of P(h) for the honeycomb lattice
in Figure 5, excerpted from Appendix A, reveals that it shares many
of the features found in the square net, including having a dip

at Tc'

One can also investigate P(h) for a triangular net. A problem
arises however since for a sixfold coordinated lattice, many more
correlations are needed. While the even spin correlations can be
obtained (Choy 1984), not all of the odd spin correlations can be.
Hence the results in Figure 6, also excerpted from Appendix A, are
incomplete. Nevertheless, one can clearly see similarities between
these results and the earlier two. In particular, there is once
again a pronounced dip at Tc'

The existence of the h = O dip means that at TC most spins
find themselves in the presence of a nonzero local field. This
suggests that perhaps small domains form in the system before
long range order sets in. The existence of strong short range order
is consistent with the value at TC of <Oool>, which for these lat-
tices is 0.71:.06 (Domb 1974). One might expect that for systems
which require short range order stronger than this before long
range order sets in, the h = O dip should be even larger at Tc'

A lattice with a larger value of <0001> at TC can be formed
following a procedure described by Syozi (1972). One starts by
decorating the bonds of a honeycomb lattice with two sites and then
one performs a star-triangle transformation, eliminating the original

honeycomb sites (see Figure 7). This process can be carried out

45

 

 

 

   

 

 

 

 

 

 

T/TC=0.0
0.5 *-
(D ,
'fi
T/TC=0.9 g
0.5-- ?
é
O ..l /_
P(h) T/T -
—|.0
0.5 - C 4
9
O a z 21 Z
T/TC=I.2
0.5 ‘
7 r
O 2 a a
13/1}:=‘33
0.5- «
O E] g g 711
'3 'l I 3
h

Figure 5. P(h) for the ferromagnetic honeycomb lattice.

46

 

   

 

 

 

 

 

 

 

E
I
T/TC=0.O g
0.5 ~ g ..
é
/
O ,
0.5 ~ -
0
P(h) T/TC=I.0
0.5 - -
0.5 L T/Tc=l.2 .
0 n m a m m n m
T/TC=CD
0.5 - .
on... m a E a .1 ..__1
’6 -4 -2 0 2 4 6
h

Figure 6. P(h) for the ferromagnetic triangular lattice.

47

-Decoration+

 

Star-Triangle
+

Figure 7. Generating the 3-12 lattice from the honeycomb
lattice.

48

mathematically, generating exact results from analagous honeycomb
results. The new lattice, called the 3-12 lattice, is more open
than the previous one and hence requires more short range order,

<0001>T = 0.8751, before long range order sets in. Note that
c

for the honeycomb, <0001>T = 0.7698. A comparison of P(h) at TC
c

for these two lattices is made in Figure 8. As anticipated, the
dip is bigger for the 3-12 lattice. '

From the two dimensional lattices, one sees that the profile
of P(h) above TC reflects short range order in the system; there
appears to be little effect due to coordination number. Below
Tc’ the rapid onset of the magnetization dominates. Finally, it
should be noted that many of these calculations can easily be
modified for antiferromagnets. For instance, a square net or
honeycomb structure can be broken up into two sublattices so that
at T = 0 the system can order by designating one sublattice as a
spin up sublattice and the other as a spin down sublattice. The
spin up sublattice will have a P(h) precisely as in Figure 4
(for a square net), while the down sublattice will have the same
profile but with peaks on the -h side below Tc' The resultant
P(h) for the entire system is symmetric, but below Tc’

Pi(h) z P(h) . (5)

The antiferromagnetic triangular net does not decompose into
two sublattices and hence cannot order. Calculations by Choy (1984)

indicate that P(h) for this system changes very little with tempera-

ture, which is consistent with its inability to order.

49

 

0.3 .
P(h)

0.2 ‘

 

 

 

 

 

 

 

'3 -l 1 3

 

 

 

 

 

 

 

 

’1
..4

-3 -1

H E

Figure 8. Showing P(h) at T for the honeycomb(upper) and 3-12
(lower) lattices.

6. P(h) for Ferromagnetic Systems of Higher Dimension

There are no exact solutions for Ising models on standard
three or four dimensional lattices. However, it would still be
interesting to obtain some information, if only approximate, on
systems of dimension greater that two in order to gauge the effect
lattice dimension has on P(h).

Frank et. al. (1982) developed an approximation method which
effectively allowed them to estimate the moments of P(h) for a variety
of three dimensional lattices at TC. In principle, one can take
their moments and recover the original distribution. In practice,
the error bars associated with those moments were too large, so that
negative weights appeared in the reconstructed distribution.

Since it appeared that little progress would be made analytically,
Monte Carlo simulations were developed.

The temperature scale for these simulations was set by the
values for TC based on high temperature series and reported else—
where: kBTC/J = 2.7044i.OOOl (diamond), 4.5108i.0002 (simple cubic),
6.3533 1.0010 (body centered cubic), and 9.79521.0005 (face centered
cubic) from Domb (1974), and kBTC/J = 6.6817i.0015 (four dimensional
simple cubic) from Gaunt et. al. (1979). No attempt was made to
locate TC precisely with these simulations since detailed information
(such as critical exponents) was not sought.

Simulations generally have the most trouble converging when

T = TC. However, since it is known that P(h) is symmetric at Tc’

50

51

only PS(h) need be extracted at the critical temperature. This
function couples to the energy, or equivalently to the nearest
neighbor correlation <0001>, which is less sensitive to critical
fluctuations than the magnetization and Pa(h). The criterion for
convergence at TC was thus taken to be that <0001> should equili-
brate. This correlation can then be compared to estimates obtained
from high temperature series.

The numerical simulations were performed on lattices of varying
sizes and with periodic boundary conditions, employing the standard
Metropolis algorithm (Metropolis et. al. 1953). Technical details
associated with these studies appear in Appendix C. A given system
would initially be brought into equilibrium and then be allowed
to remain there while data was accumulated for thermal averages.
Statistical errors quoted refer to the fluctuations in the thermal
average during the second half of the averaging part of the simula-
tion. The simulation was performed on three different lattice sizes
(four at Tc)’ with the largest system containing roughly 30,000
spins. Then the data was plotted against the inverse of the num—
ber of spins in the system to allow extrapolation to infinite
size. Figure 9 shows an example of an extrapolation of <0001>
for a simple cubic at Tc' The size of the data points is of the
same order as the statistical error.

The results of these extrapolations can now be compared to
high temperature series results (see Domb 1974 and Sykes 1979)
to gauge the accuracy of the simulation. This comparison, reported

in Table 1, shows reasonable agreement between the two results.

52

.404

<OOOI>

.34-l

-32‘

 

 

0.0 .2 I .0 .8 1.0

Figure 9. Extrapolation of <0 0 > to the infinite lattice
limit for a simple cubic at Tc'

53

Table l. <OOOI>T as obtained from Monte Carlo simulations and
c

compared to high temperature series results.

 

Lattice High T Series 1 Simulation
Diamond .437 t .003 .44 i .01
SC .3307 t .0001 .34 t .01
BCC .2732 i .0002 .28 t .01
FCC .2474 t .0001 .26 t .01
Hypercubic .188 t .003 .20 i .01

 

 

54

Based on this, one would expect the data for P(h) to be accurate
to within a few (3 or 4) percent, which is sufficient to determine
the general nature of P(h).

Figure 10 displays P(h) for four different three dimensional
lattices. All of these lattices show an extremely flat profile
near Tc' Whether there is a slight dip or bump at h = 0 cannot
be conclusively determined from this data, but it is clear that
these profiles differ from those of the two dimensional lattices.
It is also evident that in three dimensions as in two, coordination
number plays very little role in the profile.

To see the dimensional dependence explicitly, it is construc-
tive to look at the simple cubics in dimensions one through four
(i.e., the linear chain, the square net, the simple cubic, and
the hypercubic). These lattices all have coordination number equal
to twice their dimension. The infinite dimension limit yields
the mean field result discussed in Chapter 4, and will be included
with this group too. These lattices have critical temperatures
scaling roughly with coordination 2, providing a consistent way
to normalize all the temperatures at which P(h) was measured (either

analytically or by simulation). More exactly, one finds

O S TC/z S l (1)

where the lower limit corresponds to the linear chain and the up—
per limit to mean field theory.

The graphs in Figure 11 show P(h) for each of the five systems
measured at each of the five (renormalized) transition temperatures.

The profiles for any group of lattices above their transition

0.5 1
0.0 J
0.5 1
0.0 —
0.5 1
0.0 4

0.5 ‘

 

0.0 n

0.5 y

0.0 J

0.5 1
0.01
0.5 1
0.0 -
0.5 1
0.0«
0.5 1
0.0«

0.5 4

 

0.0 -1

Figure 10. P(h) for four different three dimensional ferromagnets.

55

DIAMOND

T/‘l’c = 0.0

 

I * r

T/Tc=0.9 a

T/Tc = 1.0

a—U—J—Z

T/Tc=1.2
(.1 a 3 9L2:

171'c = on

-4-2024

 

 

0.5

'0.0d

0.5 ‘

0.0

0.5

0.0

0.5

0.0

0.5

0.0

0.5 ‘

0.0

0.5 ‘

0.0

0.5

0.0

0.5

0.0

0.5 1

0.0

4

J

 

 

 

l

 

 

T/Tc = ..

42-8404 312

h

56

.uHEHH tame“ came mnu 00m tam 430w :mzoucu mco mcowmcmsflc

OJ

'1

 

 

 

N\__

 

 

 

 

 

—_ o
—1— '-
_—_ o
——9-
__1_
fi— 0-
_J_
__L .,
—

 

 

 

 

 

 

 

 

l _ fl _ l T111 4 .11--.._I__ _ l _ _ fillrli
I _ m _ 1. 11.11 _‘-._-..-4_.--_-,-_1 _ _. _ L . .
1.34011 111-31-. 4.1.111 _ . _ 1 _ L
p _ L 11L 1!! x. p - L III. P L _ _ I.

 

 

 

 

. 1.30.. u is ovum." as...» u is 22.1 22°» 3. 22... 1 Each 1 S. o a 3? u w:

IO

(p
IO

{p
10

Fr

1O

f

no

 

mo moansu oaaefim now Acvm .HH oszwflm

E

am

an

3

”E.

57

temperatures are very similar provided they are measured at the same
renormalized temperature. The most dramatic example of this is

seen by looking at the one, two, and three dimensional data at

a temperature corresponding to TC in three dimensions. All three
profiles are extremely flat.

Below Tc’ the symmetry is broken and the profiles lose much
of their similarity. In the case of the four dimensional hyper-
cubic, the approach to mean field theory can be seen by noting the
development of a peak away from h = Jz.

Clearly the coordination of a lattice affects P(h) in that it
determines the number of 0 functions. It additionally helps set
the temperature scale. Beyond that, the profile of P(h) is sensi—
tive more to the dimension than to the coordination. At Tc’ this

dimension dependence comes chiefly through <0001>T , reflecting
c

the amount of local order required before long range order can

be established. Thus in four dimensions, P(h) at Tc is more nearly
Gaussian than in one or two dimensions because less short range
order is required in four dimensions for long range order to set

in.

7. The Sherrington-Kirkpatrick Spin Glass Model

Having looked at P(h) for a wide variety of ordered systems,
it is now useful to turn to a disordered system, in particular, to
a spin glass model. The model used is based on the Edwards and
Anderson (1975) concept of a spin glass arising from some form of
random interactions among spins. In the spin glass phase, the spins
appear to be randomly oriented, but they are in fact frozen in. If
one defines a quantity q by

q = lim <Si(O)Si(t)> , (l)

t+m
then it becomes nonzero in the spin glass phase. Edwards and Ander-
son argued that at least within mean field theory, this model is
capable of reproducing some of the characteristics of experimentally
observed spin glasses such as a cusp in the susceptibility which

is rounded by the application of a field. Whether or not the Edwards
and Anderson approach to spin glasses is an accurate description

of experimentally observed spin glasses is not the focus of this
investigation. Rather, as before, known results of a model will

be exploited to calculate P(h). This form of P(h) will look substan-
tially different from those found in the ferromagnets. Having ob-
tained P(h), one can then perform a neutron scattering experiment

on a disordered Ising system (preferably a spin glass) to link up
with experiment.

The spin glass model to be investigated is that introduced by

58

59

Sherrington and Kirkpatrick (1975). They assumed a Hamiltonian of
the form

H = -% Z Jijoio. (2)
i,j 3

with the matrix of Jij's symmetric and having vanishing diagonal
elements. The Jij's are otherwise independently distributed random
variables with mean Jo and variance J2. Since this system has infi—
nite range interactions, one might hope some sort of mean field theory
would work. The infinite range interaction also means some variables
must be rescaled with the number of spins, N:

32 JZN (3)

O JON . (4)

3

The energy scale is then set by 3 or 30.
When this model was introduced, Sherrington and Kirkpatrick

proposed a solution based on the replica method. Starting from the

identity

in z = %ig(Zn—l)/n (5)
where

Z = TrIeXP(-BH)] (6)

the calculation involves looking at Zn for integral n. Some algebra
reduces this to n replicas of the original system with an effective
interaction between a spin in one system and the same spin in a rep-
lica.. The new system, exactly solvable, depends on the integer n.
One then analytically continues Zn to n=O, recovering the free energy

of the original system.

While there appear to be no problems with this approach in the

6O

paramagnetic phase, it yields incorrect results in the spin glass
phase (e.g., negative entropy). The analytic continuation of Zn
with respect to n and replica symmetry breaking are thought to be

at the heart of the problem. The phase diagram obtained by the rep-
lica method is shown in Figure 12. Due to the problems with the
calculation outside the paramagnetic phase, the spin glass-ferromag-
netic phase boundary may not be correct. However, since most subse-
quent studies have been limited to JO=O, this appears to be the

most accurate phase diagram available.

Although Parisi (1980a) has devised a method for obtaining exact
convergent series for various thermodynamic quantities in the spin
glass phase, the less complete solution due to Thouless, Anderson
and Palmer (1977), hereafter referred to as TAP, is more relevant
to this study of P(h). Looking at the SK model for 30:0, TAP compared
it to a Bethe lattice with the coordination number, 2, tending to
infinity. The interactions on the Bethe lattice were taken to be
nearest neighbor only and with a random distribution identical to
that for the standard SK model. In the paramagnetic phase, TAP showed
that the difference between the free energy per spin associated with
such a Bethe lattice and the free energy per spin of the SK model
with z spins, vanishes in the infinite 2 limit. They were also able
to show this is true just below the transition temperature and for
temperatures near zero.

The TAP equations reproduce the replica result in the paramagnetic
phase,

E/N = —%BJZ . (7)

1-4

1.2

0.3
kBT/J
0.6

Figure 12.

61

 

 

 

 

PARA
FERRO
‘ SPIN CLASS
I l I I I l I
0.2 0.4 0.6 0.8 1 1.2 1.4

30/3

The phase diagram of the Sherrington-Kirkpatrick model
(reproduced from Sherrington and Kirkpatrick 1975).

 

62

Unfortunately, their equations do not lead to exact, closed form
solutions below TSG' They do however yield some information about
the zero temperature form of P(h) for small h.
First, a digression is required to a discussion which appeared
in a paper by Palmer and Pond (1979). They assumed that
P(h,T=O) a ht ' (8)

and tried to determine t. Starting with the system in a local mini-

 

mum, the N sites are relabelled in order of increasing [hi . The
spins at sites 1 through n are flipped, where
1<<n<<N . (9)

Before the spins were flipped, the energy associated with all the
bonds comming into a given site was just —ihii. By flipping the
spins, the energy has changed by an amount
n n n
AB = 2 Z lhil — 2 2 Z Ji.oio. . (10)
i=1 1:]. j=1 J J
If the system was originally in a local minimum, then
AE 2 O . (11)

0
Consider I Z J..O.l: one expects that the signs of J.. and O.
j=1 1J J 1J J

are not correlated so that the sum should yield

0 1 1 1
| z Ji.0.| z n2J = 36/112
j=1 J J
n
However, the signs of Z Jijoj and Oi are correlated in the ground
i=1

state so that

63

II n

X Z J..0.0. ~ 303/2/N
i=1 j=l 11 l J

1/2 (12)

Returning to the first term in (10), assume that for small h,
P(Ihl) = Aht . (13)

Then,

11 h t
- I m Ah dh

N 0
A t+1
- t+1 hm ’ (14)

where hm is the largest field associated with the n flipped spins.

 

Thus
h _ .Eil.fl 1/(t+l) <15)
m - A N ’
so that
n hm t+1
z lhil - N J Ah dh
i=1 0
A t+2
N t+2 hm
' t+2 A N
= const X N %'(n/N)l/(t+l)
= const X n (n/N)1/(t+1) (16)
If the constants are dropped, requiring that AB 2 0 yields
n(n/N)1/(t+1) 2 nun/N)“2
or,
L £n(n/N) > i 2n(n/N) (17)
t+1 - 2 °

64
Since £n(n/N)<O, this requires
1221.21:
or,
t 2 1 . (18)
This is as far as Palmer and Pond took their analysis. It shows
that it is plausible that
P(h,T=O) a h
In fact their simulations of the SK model at zero temperature as
well as simulations by Kirkpatrick and Sherrington suggested that
this was the case (Kirkpatrick and Sherrington 1978).
This evidence makes reasonable the assumption by TAP that
for small h,
P(h,T=O) = 11/112 . (19)
They were able to obtain some information about the slope, H-Z.

It should be noted that their definition of P(h),

1
0mm) = If 6(h-l<hi>l) . (20)

is in general not equivalent to the definition used in equations

(3) and (4) of Chapter 3,

P(h) = Z<0(h—hi)>

l
N .

1
However, at T=0 for a system that orders, that is, when Pi(h) is
of the form

Pi(h) = 5(h-h2) . (21)

these two definitions yield essentially the same result.

PTAP(h) has the advantage that at low temperatures, it appears

65

to be less temperature sensitive than P(h) so that TAP can neglect
temperature dependence in their calculation. However, PTAP(h) is
a mean field quantity and will not necessarily yield exact results.

Returning to their calculation, they first define
mi = <Oi> . (22)
Using the fact that the Jij's scale with J = 3//; and are all hence
small, and using standard Bethe lattice techniques, they arrive at
2 2 -1
ZJ..m. - m.BZJ.. (l-m.) = k T tanh m. . (23)
J 1j j 1 j 13 j B 1

This equation is simplified by approximating

2 2 ~ .2. 2
E Jij (l-mj) - (1-mj) E Jij
J J
2 <1-q) 32 (24)

where q is the Edwards-Anderson order parameter. It will turn out

to be consistent to take the low temperature behavior of q to have

the form
q = 1 - a(kBT/3)2 . (25)
To further simplify notation, define
h: s <hi> = z Jijmj . (26)

J
Then equation (23) can be rewritten

* -1
hi = akBTmi + T tanh mi . (27)

While an exact expression for q is

'72 1 2
q = mi = N?- [fpi(h) tanh Bh dh] 9

it is plausible within the mean field limit to write

q = f m2<h*> PTAP<h*> dh* , (28)

66

where m(h) is obtained by inverting (27). No attempt is made either
here or in TAP to rigorously justify (28). Combining equations

(28), (25), and (19) one arrives at

H
1 2 h* dh*
= I [l-m ] -§' aE—'dm . (29)
O H

*
Equation (27) will provide h (m) and its derivative so that

 

2 2
k T k T 1
0[—§—} = L%%J ‘. [l-mz] [0m + tanh-1m ] [ a + ] dm . (30)

IO 1-m2
Integration and rearrangement yields

22n2+1 + 2n2
3 a

 

<wm2=%+ on

While this equation relates H2 to a, it does not determine
either of them. It is consistent with the results of Monte Carlo
simulations reported in Palmer and Pond (1979) to take the minimum

value of H2 compatible with (31). This occurs when

0 = 2V2n2 (32)
giving
(H/J)2 = 3&9%il-+ «knz . (33)

Taking care of a factor of two discrepancy arising from the dif-
ferent definitions of P(h) used, one would then expect at T=0
P(h) = 0.307h (small h) , (34)

where P(h) is as defined in Chapter 3.

67

Clearly, the approach used in TAP to calculate this slope has
been mean field. Nevertheless, as reported in the next chapter,
it meets with much success when compared to Monte Carlo simulations.
Finally, because of TAP's definition of P(h), there appears to be
no straightforward way to use this argument to investigate the low

temperature properties of the P(h) defined in Chapter 3.

8. P(h) for the Sherrington—Kirkpatrick Model

As pointed out in the previous chapter, most of the attention
paid to P(h) has been concentrated on the zero temperature problem.
However, using the TAP trick of replacing the SK model with a
Bethe lattice, it is possible to obtain P(h) in the entire para-
magnetic regime of this model. To do this, one makes use of the
fact that for a Bethe lattice in the paramagnetic phase, that is

when boundary conditions are of no concern,

<0.0. > = <0.0. ><0. 0. >
1 1+2 1 1+1 1+1 1+2

(1)

tanhBJi+1,1+2

tanhBJi’1+1

(2)
This is valid when site 1+1 is the only site crossed by the path
from i to 1+2 (see Chapter 4 or Thorpe 1982). In this model, the

temperature scale is set by

kBT ~ 3 = J/N (3)
so that
BJiJ. ~ 1M? (4)
and one can approximate
tanhBJij 3 BJij . (5)

Starting from the integral representation of the 0 function,
Pi(h) = Ida exp(-ih0) < exp(ihi9) > (6)
or

Pi(6) = <exp(ihie)> , (7)

68

69

and using

h1 = E Jijoj . (8)

(where the sum is over the z neighbors of site i), it follows that

Pi(8) = <g exp(iJijOj6)>
= <g [COS(J1je) + iOJsin(Jij9)]>
= <§ cos(Jije)[l + intan(Jij5)]>
= <§ cos(Jij6)[1 + inJij6]> . (9)
P(B) is obtained by averaging over all configurations of {Jij}'

That average will be denoted by an overbar. Thus
P(B) = 51(8)

<H cos(J..8)[l + iO.J..6]>
j 13 J 1J

 

 

----- z .
cos(Jij8) [ 1 + 12(Jij6<0j>)

 

+ 12.212211 J..J. 62<0.0 >
2! 1j 1k 3 k

 

 

13 z(z—l)(z—2) J..J J. 03

 

 

+ 3! 13 ik 1l <OijOl>
.4 z(z:l)(z-2)(z-3) 4
+ 1 4! JijJikJilJime <0j0k010m>
+ ...] + 0(l/z) . (10)

The order l/z corrections, arising from averaging over all the
cos(Jij9) terms independently of the Jij<...Oj...> terms, are
neglected.

In the paramagnetic phase, the thermal average of an odd number

of spins is zero. Furthermore, using equations (2) and (5),

70

~ 2
<0j0k> - B JijJik (11)
Thus, in the large 2 limit,
2 ___ 2 4 ___ 4
----- z z 2 2 2 z 4 2 4
P(B) — cos(Jij8) [1 - 578 [Jij 6 + 278 [Jij] 6 - ...] . (12)

 

This last step has used

 

fi——— 2
2 2 [J2 ]

JijJik 3 ij (13)
which holds since Jij and Jik are independent random variables.
The infinite sum in (12) is the Taylor expansion of cosine,
so that
----- z 2
P(G) - cos(Jij6) cos(BzJij9) . (14)
Using notation and scaling introduced earlier,
J2. = J2 + J2
ij 0
= 32/2 + 33/22 . (15)
For large 2,
'72 2
2 ij - 3 (16)
Furthermore,
cos(J 6) z = (1 - ABZI— + )2
ij 2 ii
2 2
6 3 2 z
[1 - 22 + 0(1/z )]
= exp(-0232/2) (17)

and hence

P(e) = exP(-6232/2) cos(8326) . (18)

71

Taking the inverse fourier transform, one obtains

1

exp(—0232/2 - 02/232) cosh(Bh) . (19)
m3

P(h) =

 

Equation (19) is valid throughout the entire paramagnetic phase
and interestingly contains no explicit Jo dependence.
Figure 13 shows a plot of P(h) at three temperatures. It

was assumed that 30 < 3 so that (19) would be valid at k T = 3,

B

which is then the spin glass temperature. The infinite temperature

limit is, as expected, a Gaussian:

P(h,T=w) = exp(—h2/232) . (20)

 

V203
When the temperature is decreased, the Gaussian broadens and flat—

tens, a phenomenom made clearer by looking at the small h expansion,

 

2
2 2 h 2 2
P(h) = exp(-B J /2) [1 --———(1 - B 3 )
F211.) 232
AA 2 2 4 4
+ 32(1/8 - B 3 /4 + B 3 /24) + ...] . (21)

The fact that the h2 term vanishes precisely at TSG indicates
how extremely flat this function becomes for small h.
This calculation of P(h) cannot be repeated in an ordered

phase for several reasons, among which are

<0.) 1 0
J

<0.0 > 1 <0.0.><0.0 > .
j k j 1 1 k

Without being able to express higher order correlation functions in
terms of the nearest neighbor correlation functions, there appears

to be little hope for finding an analytic expression for P(h) in

72

 

0.5
— kB'r/CI = 1
0 4 ---- kBT/3'= 2
0.3«
P (h)
0.2 «
0.1 1
0

 

 

 

 

h/J

Figure 13. P(h) from the analytic result for the Sherrin ton-
Kirkpatrick model in the paramagnetic phase ( o < 3).

73

the spin glass phase of the model. For this reason, a Monte Carlo
simulation was employed to study an ordered phase, specifically,
the spin glass regime.

The probability distribution for the Jij's was taken to be
of the form

p(Jij) = a0(Jij - J) + (l - a)0(Jij + J) , (22)

where

a = %(1 + Jo/J) . (23)

This form ensures

J.. = J
1j o
and
?.=.I2
13

Since it was explicitly seen in the paramagnetic phase that once
Jo and 32 have been fixed, the details of the Jij distribution are
irrelevant, it is clear that in this phase one is free to choose
the particular form of p(Jij) for convenience. The same statement
cannot be proved in the spin glass phase but is believed to be
true there owing to the infinite range of the interaction.

The advantage of the distribution defined in (22) is that
now Jij can be represented as a single binary digit. This fact
was exploited to increase the effective memory available for the
computer simulation and to increase the execution speed of the
simulation. The details of this technique are provided in Appendix
C. Unless otherwise noted, the data displayed in this chapter

results from an average over four different 1020 spin systems.

74

In Figure 14 the results of a simulation for Jo = O at T = 2T
are compared to the analytic calculation, and a similar comparison
at T = TSG is made in Figure 15. While there is some noise in these
simulations, it is clear that they are in essential agreement with
the analytic result.

1 In Figure 16, P(h) as obtained from the simulation is plotted
at four temperatures. In keeping with the trend established above
TSG’ a dip at h = O has now developed below TSG'

The data at T = O was accumulated by taking forty different
{Jij} configurations and searching for one state in each that was
stable against single spin flips. Kirkpatrick and Sherrington
(1978) and Palmer and Pond (1979) have already studied the zero
temperature problem in detail, so no attempt was made to be more
careful in searching for local minima. Those authors, in fact, report
lower values for P(h=O,T=O), and it is believed that for an infinite
system, P(0,0) = 0. However, even using states stable against
only single spin flips, the small h slope obtained from this simu-
lation, 0.311/32, is in excellent agreement with the TAP (1977)
prediction of O.3O7/J2 discussed in the previous chapter. Hence
the slope appears to be less sensitive than the intercept to the
particular local minimum one is in.

As the temperature is reduced from TSG/Z to 0, most of the
change in P(h) occurs for small h. The large h regions are relatively
stable. This suggests that the formation of the h = O dip is an
important characteristic of the spin glass phase.

To study P(O) in more detail, first note that above TSG’

75

0.4

 

—— SIMULATION
---- EXACT

0.3 -1

P(h) 0.2 -1 \

 

 

Figure 14. A comparison between the Monte Carlo simulation
of the Sherrington-Kirkpatrick model and the

analytic result: T = ZTSG'

.4

 

76

 

 

 

 

0.4
— SIMULATION
‘ ---- EXACT
0.31
P(h)
0.2 -
0.1 —
0 .
-6 6

 

Figure 15. A comparison between the Monte Carlo simulation
of the Sherrington—Kirkpatrick model and the analytic

result: T = TSG'

77

 

   
      

 

 

 

0.4
- — kBT/J = 2.0
""’" kBT/J =1.0
0.3 _‘ """ k317i: 0.5
P (h) " ,
0.2 “ 1' . "
J
0.1 -
C) I 1’50 1 I I ‘1
-6 -4 -2 0 2 4

Figure 16. Results of the Monte Carlo simulations of the
Sherrington-Kirkpatrick model: P(h) at four dif-

ferent temperatures. The T=O data has been
symmetrized.

78

 

P(h=O,T) = 1 exp(+82J2/2) (24)
/203
so that at TSG = J/kB,
P(h=O T ) = -l——-e-% (25)
30 55:1
It is also easily shown that
_ _1
BPgh-OzT) = l e 2 (26)
18C V_fi

Since it is believed that the specific heat and energy are
continuous at TSG for this model, these values hold for just under
TSG also. If one adds to this the assumption that

P(0,0) = O , (27)

there are three constriants on P(h=O,T) for O 5 T g T G' One solution,

S
but certainly not the only solution, is to try the linear form
k
_1
-l——-e 2 B T . (28)

P(h=O,T<T
7213

so) ‘

Figure 17 compares this function to the data from Monte Carlo simu—

lations. The deviation at T = O is probably, as discussed before,

a finite size effect. Attempts to study the size dependence in

the rest of the regime were not conclusive. It is clear, however,

that even if (28) is not exact, it is at least a good approximation.
It was noted in Chapter 5 that for the antiferromagnetic square

net below Tc’ P(h) is symmetric while Pi(h) is not. That is, when

the phase transition occurs, a symmetry is broken so that there are

two distinct sublattices, one with a positively biased Pi(h) and

the other with a negatively biased Pi(h)' One might expect a similar

symmetry breaking for this model. While an extensive study was not

P (0)

 

79

 

I
0.3-
O/
02 9”
n “'1 /9/
/Q/
9/
o//
0.1 - 9/
0”
//
q o I)/
/
/’
0.0 I j
0.0 0.5 1.0
kBT/J

Figure 17. P(O) as a function of temperature in the spin

glass phase of the Sherrington-Kirkpatrick model.
Open circles indicate simulation results. The
dashed line represents the linear hypothesis.

80

made, a few simulations of a 540 spin system (JO = 0) yield the
tentative conclusion that

Pi(h) = P(h) T > T (29)

SG

Pi(h) 1 P(h) T < TSG . (30)

The detailed nature of this broken symmetry is not as of yet clear,
although one can say that it is substantially more complicated than
the two sublattice state of the antiferromagnet.

Following the temperature development of P(h) has allowed the
tracking of its development from a Gaussian function to a function
with a hole in the center. It has been seen that the dip develops
precisely at TSG and the system very quickly establishes the large
h states. Below about TSG/2, most of the action in P(h) occurs
in the small h region.

There has been some recent work (Mézard et. al. 1984) sugges-
ting that the relevant order parameter in this system is not q
(the Edwards—Anderson order parameter), nor is it q(x), a generalized
probability function introduced by Parisi (1980a). Rather it has
been suggested that the system is best described by a probability
distribution governing the distribution q(x). That is, by performing
a configurational average to calculate q(x), one is throwing away
valuable information in the spin glass state. One might be con—
cerned then that the Monte Carlo simulations performed configurational
averages and that what should have been measured was the probability

distribution of the distribution, P(h). However, it has been noted

81

earlier that one cannot calculate q directly from P(h); one needs

to know Pi(h) at all sites. Thus q is not expected to couple strongly
to P(h). It is not unreasonable that even though q and q(x) may

not be well defined thermodynamic variables, P(h) may nevertheless

be well defined for this model.

9. Future Work

While the nature of disorder in Ising systems has been examined
here, the major development has been that of a new tool, P(h), to
study Ising systems at all temperatures. This function leads to
exact thermodynamic results as well as having a measurable symmetric
component.

An understanding of this function for a wide variety of Ising
ferromagnets is in hand. However, only one disordered system, the
Sherrington—Kirkpatrick spin glass model, has been studied. While
some useful conclusions were drawn, it was not possible to determine
which general characteristics of P(h) were common to Ising spin
glasses in general and which were peculiar to infinite range spin
glasses. Therefore, it would be useful to look at an Ising spin
glass model with short range interactions. It is reasonable to
expect that a system that orders will have a dip at h = 0, but
beyond that it is difficult to generalize on the basis of the
Sherrington-Kirkpatrick model.

One can extend the definition of P(h) for Ising systems to
a P(H) for classical KY and Heisenberg systems. While the relation—
ship between this function and the inelastic neutron scattering law
is not clear, inspection of the proof in Chapter 2 reveals that a
similar relation between the energy and P(h) should hold. It was
crucial in that proof that all parts of the Hamiltonian commute with

each other, and this is the case for these classical systems.

82

83

Extension to these systems would allow for the simulation of
a Heisenberg spin glass. One possibility would be the EuxSr1_XS
system. Eu2+ ions have spin 7/2 and lie on an FCC lattice. They
have ferromagnetic nearest neighbor interactions and antiferromag-
netic next nearest neighbor interactions. EuS acts as a ferromagnet,
but dilution with the nonmagnetic Sr2+ ions enhances the competition
between the nearest and next nearest neighbor interactions, creating
a spin glass at appropriate concentration and temperature. This
system has been well studied both experimentally and by simulation
(see, e.g., Krey 1981). To the extent that spin 7/2 can be approxi-
mated as classical, EuxSrl_XS then provides a prime candidate for
studying the temperature dependence of P(h) in a well characterized
system. A connection between S(k,w) and P(h) could be extremely
useful in comparing simulation results to experiment.

Another system worth investigating consists of classical XY
spins on a triangular lattice with nearest neighbor antiferromag-
netic interations. A study by Lee et. al. (1984) indicates a very
rich phase diagram in the H-T plane. With this single system it
would be possible to sample several different phases.

It may also be possible to extend the definition of P(h) to
the spherical model (Berlin and Kac 1952). This would provide
access to the well studied spherical infinite range spin glass
model (Kosterlitz et. al. 1976).

A more difficult extension of this formalism is to quantum
Heisenberg or XY models. Problems arise early on since the local

field operator,

K. = 22.1.15. (1)
1 . 1] j
J
is not diagonalizable. Hence,
P10?) = <60? - Ki» - (2)

is not well defined. Preliminary efforts to relate the Ising P(h)
to P(IHI) or to P(hz) for quantum systems have not yielded a use-
ful formalism. It could be that there exists a relation between

the energy and some generalization of the Ising P(h) that is more
complicated than a simple integral. Any such hypothesized relation—
ship can be checked against the known results for the infinite quan-
tum KY chain (see, e.g., Lieb et. al. 1961) or for finite length
quantum Heisenbero systems (Bonner and Fisher 1964). An extension
of the P(h) formalism to quantum systems, while difficult, has the

potential for providing new insight into these models.

Appendix A.
Local Magnetic Field Distributions I.

Two-Dimensional Ising Models

Local Magnetic Field Distributions I.

Two-Dimensional Ising Models

M. Thomsen+
Department of Physics and Astronomy
Michigan State University

East Lansing, MI 48824 U.S.A.

M. F. Thorpe++
Cavendish Laboratory, Madingley Road

Cambridge CBB OHE, U.K.

T. C. Choy and D. Sherrington
Physics Department, Imperial College

London SW7 ZBZ, U.K.

Copyright 1984 The American Physical Society
+Exxon Fellow

++Permanent address: Department of Physics and Astronomy, Michigan State

University, East Lansing, MI 48824, U.S.A.
85

86

Abstract

We show that the Statistical mechanics of Ising models :an be conveniently
reformulated in terms of the local magnetic field probability distribution
function P(h). It is shown that for arbitrary exchange interactions J44,

~J
which may or may not be random, boch thermodynamic quantities like magneti-
zation, specific heat etc.; and the neutron scattering law S(k,0) can be
obtained from P(h). Indeed S(k,m) provides a direct measurement of the
symmetric part of P(h) which also determines the energy, specific heat etc.
while the magnetization can be obtained from the antisvmmetric part of P(h).
As an example, specific results for P(h) are presented for the honeycomb,
square and triangular lattices with constant nearest neighbor interactions.

All three lattices exhibit a pronounced dip in the center of P(h) at the

transition temperature.

87

I. INTRODUCTION

At various stages in the development of the theory of magnetism, the
local field probability distribution function P(h) has been used}.5 Despite
the fact that P(h) has a simple physical interpretation, it has 323 been
widely adopted. This is probably because its usage has been associated with
mean field type theories and as such it is regarded as an approximation.

In this paper we show that an exact, useful and complete description
of the thermodynamics of a rather general class of Ising models can be
obtained in terms of P(h)- These models have localized spins Si associated
with sites R1. The spins can have different magnitudes and the sites R1 do

not have to form a crystalline lattice. The Hamiltonian involves only the

S2 components and.contains arbitrary exchange interactions Ji and external

1 1
fields Hi' The range of the Jij is also quite general. All these results
are scattered around in the literature of the early 1960'31-4 and are

A

summarized by Southern.3 We collect them together in this paper.

We show that a measurement of the inelastic neutron scattering law
S(k,m) at a particular temperature provides a direct measurement of the
symmetric part of P(h) at that temperature. The total magnetization of
the system can be obtained directly from the antisymmetric part of P(h),
although the converse is not true. Finally the total energy can be obtained

from the symmetric part of P(h), if all the Hi = O. Other thermodynamic

 

quantities of interest, like the free energy, specific heat etc. can be

‘H

88

obtained from the energy. The above connections are very powerful because
they require a knowledge of only P(h) and the temperature and no other infor-
mation about the system. In the case where the Hi # 0, the energy, and

heats the other thermal thermodynamic quantities, can still be found from
P(h) in simple cases. Such an example would be a uniform field where

all the Hi a H, and a knowledge of P(h), the temperature and H would lead

to the energy.

The layout of this paper is as follows. In the next section we develop
the general formalism for S s 3. This is for pedagogical reasons and the
rather obvious generalizations to arbitrary spin are given in the Appendix.
In Sec. III we illustrate these results for various spin 3 nearest neighbor
Ising models; the one-dimensional chain and the honeycomb, square and tri-
angular net in two dimensions. These results are exact and could have
been obtained at any time in the past 20 years. Surprisingly they do
not appear to be in the literature. The results in two dimensions are
rather interesting as P(h) changes from a roughly Gaussian shape at high
temperatures to develop a pronounced dip in the center at the transition
temperature Tc. It is-unclear how general this phenomenon is. We define
a quantity that measures the fluctuations in the local field h and show,

not surprisingly, that it has a maximum at Tc.

89

The attraction of P(h) is that it provides an alternative description
of a thermodynamic system with perhaps a more intuitive interpretation than
the free energy. It contains more information than just the free energy and,
for example,the neutron scattering law S(k,m) can also be derived from it.
We note that S(k,m) cannot be derived from the free energy.

Our original motivation for this work was that computer experiments on
spin-glass models at zero temperature often focus on P(h). Indeed it is
a "natural" quantity to compute. Such experiments have universally demon-
scrated a zero-field minimum in spin-glass models with competing interactions
in their ground or lowblying metastable scate5.: A similar minimum has
. . . . . - .- 7 3 .
oeen found in Simulations of amorphous antiferromagnets. ’ Because or the
interescing new results found for simple Ising models, we have deferred a

discussion of more complex spin systems with competing interactions until a

subsequent paper.
II. GENERAL FORMALISM

We define a spin 8 Ising system with the Hamiltonian

H. 0. (l)

H.
U
H
U
P:
LJ
HIvI
H
H

where the Ji" Hi are arbitrary and the factor 3 is to prevent double counting.
This Hamiltonian describes a completely general spin 8 Ising system with

0i - : 1, an arbitrary range of exchange parameters Jij and an external field
Hi that can vary from site to site. We define the Jii * O. The sites can be
inequivalent and we make no assumptions about the existence of a crystalline

lattice or translational invariance. Although we develop these results for

spin 8, all the results of this section easily generalize to arbitrary spin

. ..r"

90

as shown in the Appendix.
There are a number of ways of decomposing the Hamiltonian (l). The

most obvious is

where

L4. 1

This is not useful in discussing local fields whereit is necessary to decompose

(l) as
H--hioi+H' (4)
where

The first term in Eq. (4) contains all terms involving 3i and other terms
are lumped together in 3'. Notice that hi as defined in (3) and hi as defined

. l .

1n (5) differ in the factor-E in front of the exchange term.. This means that
some care must be taken in calculating the energy of the system later in this
section.

We consider the thermal average ( 00, > where O is any operator not

involving site i;

'1

91

< 0": > - < Tjr [0°11 exp (311131)] / TirEexpmhicifl >

- < 0 tanh Bhi > . (6)

7
Putting 0 equal to the unit operator we find that

9

< 31 > - < tanh Bhi > a < tanh 3(é Jij Jj + Hi) > . (7)

J

We note in passing that replacing the thermal average of the tanh in (7) by
’3
the tanh or the thermal average leads to mean field theory. Putting 0 equal

to hi in (6) leads to

Jijgi‘j 11

La. ('1. 1

a < hi tanh 8h, > (8)

It is not possible to obtain the energy from (8) alone for H1 # 0 because

of the missing factor-% in front of the exchange as required for the

decomposition in Eqs. (2) and (3). However the energy can be obtained as
(9)

l ,
- 2 < (hi + H1) tanh 3hi > . x

The magnetization M and the energy E are given by

I".

92

K
n
Hd\4
A
Q
P.
V

(10)

We introduce the probability distribution Pi(h) for the local magnetic

field at site i by

Pi(h) - < 0(h - hi) > (ll)

where both the "internal" and external fields are included in hi via Eq. (5).

It is of course possible to define a similar quantity Pi'(h) with only the

internal field contribution counted;

' I — I —
P1 (h) <50. §Jijoj)> Pi(h 111) . (12)
For the rest of this paper we shall use (11). From the definition of
P,(h), we see that
jPi(h)dh- 1 , (13)
<'TJ J+H)n>=‘r'nP(h)dh (1')
‘1; ij 1 1 ,1“ 1 1 “
and,from (7) and (9),
< oi > a I tanh 3h Pi(h) dh (13)
J
(16)

A.

h tanh 3h Pi-I’h) dh - 1% Hi < I; >

m
I
I
NIH

’4-

93

The distribution function P(h) for the whole system with N sites is

obtained from (11) via

P(h) - §- 2 Pi(h) . (17)
1

It can be seen from (15) that a knowledge of P(h) and the temperature is

sufficient to determine the magnetization ;

u . N ; tanh 3h P(h) 1h . (18)

It is not necessary to know anything else about the syscem; in particular
it is not necessary to know the Jij or the Hi. It is this sense in which
P(h) provides a description of the system in the same way as the free energy

does. Because tanh 8h is an odd function of 8h, only the antisymmetric part

Pa(h) of P(h) contributes to M:
l .
Pam) =- -2- [P(h) - 14-10] . (19)
(20)

u a N j tanh 3h Pa(h)dh

In general the energy cannot be obtained directly from P(h). However ii,

all the Hi are known to be zero, then from (10) and (16)
N I
E a --5 J h tanh 3h P(h) dh (21)

and E can be calculated directly from P(h) without any knowledge of the Jij'

 

94

Because (h tanh 8h) is an even function of h, only the symmetric part Ps(h) of

P(h) contributes to E ;

930:) %IP(h)+P(-h)1 . (22)

le

I
I h tanh 3h 9,01) dh . (23)

J

In this case the free energy F is given by

L)

If the Hi are sufficiently simple, the free energy can still be found.

For example if the magnetic field is constant at every site Hi - H,the
energy given in (9) becomes
\1 1’
E - -13 J (h + H) tanh 3h P(h) dh (25)

and the free energy F(H,T) is given by {24) with the integration done at
constant H.

It can be seen that under appropriate conditions the local magnetic field
probability distribution P(h) leads to the magnetization M and the energy B.
All other thermodynamic quantities like the free energy, specific heat etc.
have to be found by integrating or differentiating M or E. This is because
M and E are determined from the expectation values of local operators
whereas other thermodynamic functions are not. Thus M and E play a special
role when the statistical mechanics is done via P(h). This point appears not to
have been appreciated by Klein and Brout4 who attempted to write the free energy

as an integral over P(h) [see their equation (2.5)].

‘u

The function P(h) contains more than just thermodynamic information; it

. 10
is directly related to the inelastic neutron scattering cross section,

+09
.. 1 I" _' q .. .5 -> * .
Sac...» a ,—‘- - e ‘M at; exp[i 14- (R. - R.)]< J.‘ LYCc) > (26)
in j ij i 3 i J
-Q

where x is any direction perpendicular to the z axis. In the previous part

. z
of this section we suppressed the z superscript on the 3i operators. This

cross section is k independent because Ising systems have no dynamics, so

T

- . . . .1
that transrorming to raiSing and lowering operators ,

< o. 370;) + 37.7fm > . <27)
1. l 1.

cl.
(‘1

1
a "' .L
( , ) a"

.0

The time dependence of these operators is given by

01 (t) - exp(th) oi exp(-th) - oi exp(2ithioi) (28)

where the h, are given by (S). Combining (27) and {28) we find that

f“
s(i€,w) - '21? e‘i‘“ dt 2 < m(Zithioi) >
1
+09
3 ‘L' g efil‘”t dt 7 < cos 2th. + i 0, sin 2th. > .
2H J i i i 1

Using equation (6), this can be rewritten as

 

96

M
3(E,w) =«J= { e‘l‘”c dt 3 < cos 2th. + i tanh 3h. sin 2th.>
2w } 1 i 1 i
+ao
a 3%.j [ e-1wt dt E(cos 2th + i tanh 8h sin 2th ) Pi(h)dh
J 1
. é} [P(w/Z) + P(-«n/2)] / [1 + exp(-3w)] . (29)

which shows that the neutron scattering law provides a direct measurement of

the symmetric part P‘Ch) of P(h);
3
3(E,d) . u PSCm/2) / £1 + sxp(-dw}3 (30)

The thermal factor [1 + exp(--13m)].l ensures detailed balance.
From the properties of P(h), it is easy to show that
1 r w
EE' 1 S(k,m) [1 + exp(-Bm)] dw . 1 (31)

and if all the Hi = O

r
'% } S(k,») [l - exp(-Bw)] ndw 8 -E - (32)
Both (31) and (32) hold for arbitrary Jij' The expression (31) could be useful

in normalizing neutron data taken at different temperatures, if this data were

to be used to find E.

97

III. TWO-DIMENSIONAL ISING MODELS

In this section we present results for P(h) for some pure two dimensional
Ising systems.12 This is straightforward but does serve to illustrate some of
the points made in the previous section and provides some interesting results.

It is convenient to redefine the Hamiltonian (1) with only a nearest

neighbor Jij which is set equal to l ;

where the sum goes over nearest neighbor pairs only and the factor % is to

prevent double counting. In what follows we use a slightly modified form of

1

the notation of Choy and Sherrington.L3 As all sites are equivalent we can
ignore the distinction between Pi(h) and P(h). Therefore

P(h) I < 6(h - G.) > (34)

l J

I «\«1 H

La.

where the sum over j goes over the z nearest neighbors of atom i. This can be

conveniently rewritten as

m
p

z
P(h) 8 fL' ! exp(-ihd) dB < exp(i 9 S J.) >
” } j;l

-Q

2
3.1; . exp(-ih5) d9 < H (l + i 0, tan a) > cos2 9 (35)

-@
where we have used the an integral representation of the Dirac delta function

with the usual convergence factors implied at :w.

‘1'}.

98

Multiplying out the brackets in the product, this becomes

a ct 6(h - s) (36)

where the r sum is in steps of l and the 3 sum is in steps of l. The final

result has the form

2
P(h)- E w~ 5(h-s) (37)
33-2 °
with
2
"s“; 3 3:5 C (33)
2 r-O r

z
and the ars are given by the generating function

a
2% J exp(-ih6) d6 cosz-r a sinr a 1r

1 3 z 1
=-- ' a o’h - s . 39
2 g rs \ ) ( )

2 s--z

The factor l/2z is included so that the a:s are integers. They are easily
obtained for a given 2 by doing the simple integrals involved in (39) and are

. . ' . - . l3 -. .
given in Table l. The Cr are correlation functions defined by

q

- 3 < 3. J. . . . o > (40)

c »
..‘ . 1. 2.
(1]...L) J

r

7‘1

where the
among the

neighbors

0

r4

0

‘99

round brackets denote that only the distinct products of r operators
2 nearest neighbors are to be taken. If we label the z nearest

of an atom cyclically from 1 to 2, then for the linear chain

1

r4

' E < a; > a 2 < J, > (41)
a Z < 3 3. > = < I, 3, > .

For the honeycomb lattice 12 = 3)

(J

 

a : < 7 > a 3 < :1 >
(i)
> (32)
' ;. <<:J 3 > a 3 <31 :2 >
(13)
a E < 31 3 3‘ > a < 3 37 :3 > .
<in) J ‘ ' J

v

100

For the square net (2 a 4)

 

C a l W
o
-l = e < Ji > a 4 < 51 >
(1)
c a R" >al ’ 2 F
r. < I“ IJ 4 < 3* 32 > + < 31 73 >
(13)
c a 7 < 3 I 3, > a a < 1 3q 3 >
3 (13k) l J a l - 3
C; a f < 3L Jj 3‘ it > a < Jl 32 J3 Ia > J .
kljki)

and for the triangular net (2 a 6)

e -Z (oi>-6<ol>

l (i)

C =- y > a 1 1

2 (f‘ < 3i Jj 6 < 01 02 > + 6 < 31 J3 > + 3 < 31 04 >

J) .
c a .; < giJij > a 6 < 315233 > + 12 < 013234 > + 2 < 31 03
(13k)

C = 7 < o. o. o, o, > = 6 < o o o o > + 6 < o J o J-

A (ijkl) i 3 a Z l 2 3 4 l 2 3 3
+ 3 < Jl 02 3 Us

‘2 a E < o. o. o, J, o > a 6 < 3 o, J J; 3- >

3 (ijkim) i J a L m l - 3 4 3

c6 = ) <J_ Jj kazgm3n>=<3132 33 '34 35 36>

(ijklmn)

o- >
3

(43)

 

5 (44)

 

101

z . - . . .
From the ars and the cr’ the amplitudes or the various fields ws occurring
in (37) can be calculated using (38). Many of these correlation functions

Cr can be obtained from the literature for two-dimensional lattices. No

exact calculations of the Cr exist in higher dimensions. We include the linear

'

chain only because it is very simple. In what follows h 3 SJ, where J has

previously been set equal to one.

14
For the linear chain

,
s, 8 tanh” K

so that from (38) we have

1 2
w = - l - tan . K
o 2< h )
(4+6)
W 'W -$<l+canh21<)
2 -2 4
There is no phase transition in the linear chain so C. = O at all tempera-
tures and P(h) is symmetric at all temperatures. At infinite temperature
(K a 0) we have we . g and w2 . w ) ' %,whilst at zero temperature
(K 8 a) we have W' i O and V7 8 w 7 ‘ L.
o - -H
For the hongyggmb Lattige all the Cr can be expressed in terms of the
reduced energy E and the reduced magnetization m
S a < o J >
o l
(47)

m = < 3 >
o

102

where 00 is the spin at the site of interest. This is achieved by rewriting

(6) as 1.5-1.8.

< 0 00 > 8 A < 0(0l + s + 03) > + B < O o J 3 > (48)

2 l 2 3

where

[tanh 3K + tanh K]

J-‘It-J

(49)

- . ,, _ ..,1
B 8-: Ltann sh - 3 tanh 34

Putting 0 = 1 yields

< 31 J, 33 > = m(l - 3A) / 3 (50)
and putting 0 . 01 yields
< 31 32 > = (a - A) / (2A + 3) . (51)

Combining all these results we find that19

 

 

 

 

h . l_ 3(tanh 3K - s) * 3m(tanh 3K - 1)
i1 2 3 tanh 3K - tanh K ‘ tanh 3K - 3 tanh K
(52)
h 3.1 35 — tanh K - m(l - 3 tanh K)
=3 2 3 tanh 3K - tanh K ‘ tanh 3K - 3 tanh K
12,13

The reduced energy 5 can be written as an elliptic integral and the

90
reduced magnetization m is a known function or K. These results are

illustrated in Fig. l.

103

For any lattice at high temperatures P(h) obeys a binomial distribu-
tion. The field is h if %(z + h) nearest neighbor spins are up and %(z - h)
are down. This occurs with probability 2C1“Z + h) / 2Z so that at infinite

temperature

P(h) . E J? zcl 5(h - s) .
2 Ԥ(z + s)

s = -z
As 2 gets large this approaches closer to a Gaussian distribution.
When the temperature is lowered P(h) at first flattens and then develops
a pronounced dip just above TC for the honeycomb lattice as shown in Fig. 1.
At TC, the values of VS can be obtained by putting 5c - 4vr79, mc a O and
exp(2Kc) - 2 + /3 in (52). The distribution function rapidly develops an

3 with

asymmetry below Tc as the magnetization increases from zero as (Tc - T)
B - 1/8. At zero temperature there is a single peak at h a 3 as we would expect
when the alignment is complete.

Note that P(h) contains the critical exponents a (through a) and 3
(through m) in Eqs. (52). All the critical exponents of interest can be
obtained from these two.21

As the coordination increases, so does the number of correlation functions
that must be known. For the sguare net, the three spin correlation function
is still proportional to the reduced magnetization via a relationship

l5-18 . . .
However.it lS necessary to compute two other indepen-

similar to (50).
dent even spin correlation functions as they cannot all be reduced to just 5
via relations like (48) and (49). These involve evaluating elliptic integrals.

For the triangular net the situation is much more complex but a similar

 

reduction to elliptic integrals for even spin correlations can be made starting
- - . 22 . .
trom the Pfaffian form given by Stephenson . The details are complicated and

‘33 2
are given elsewhere‘ . While the magnetization is given by Potts , the remain-

 

104

ing odd spin correlations can be found using the methods of Barry, Minera and
Tanakazs. However, we found their results to contain errors and have not
pursued this further. Consequently, Fig. 3 for the triangular net is not
complete for T < Tc'

The results for the square net and the triangular net are shown in
Figs. 2 and 3. They are very similar to those for the honeycomb lattice
given in Fig. 1, showing a rather universal behavior for P(h) for all two-
dimensional lattices.

Results have been published previously13 by two of us for the antiferroe
magnetic triangular net which show that there is very little change in P(h)
between T 8 n and T 8 O apart from some flattening. This is related to the fact
that this model does not have a phase transition and has a finite entropy at
zero temperature. The antiferromagnetic honeycomb and square lattices map on
to their ferromagnetic counterparts as they are both bichromatic. Using the
definition (ll) for Pi(h)’ these functions will be identical above Tc for both
sublattices and equal to the ferromagnetic counterparts. However below TC, the
asymmetry will develop on opposite sides of h 8 O for the two sublattices and
so P(h) for the whole syscem is juSt the symmetric part PSCh) of P(h) defined
in (22) for the ferromagnetic counterpart. This leads to zero net magnetization
via (18) as expected. In order to get the staggered magnetization it would be
necessary to know Pi(h) at the up and down sites separately and not just the
average.

The one-dimensional Ising model can also be regarded as having a dip in
P(h) at Tc if TC is identified with zero temperature. At zero temperature
wO = 0 and w2 = w_2 8 8 from (46).

If the distribution P(h) is to be characterized by a single parameter,

the most useful is

V8[h2-(h)2]/z (s3)

105

where the bar denotes an average over P(h). From the definition of the ct,

this can be rewritten as

V8l+(2c-c2/z <54)

2 l )

and is shown in Fig. 4 for ferromagnetic interactions in the honeycomb,

 

square net and triangular net. At high temperatures v - l for all lattices

as the correlation functions c1, c2 become zero. Not surprisingly this quantity,
which measures the local fluctuations, has a maximum at T5 and then decreases
rapidly to zero as the temperature goes to zero. The increase in v as the
temperature is lowered from infinity towards Tc is due to the flattening and
then the dip in P(h). The rapid drop in vbelowTé is due to the asymmetry
caused by c1 and hence the magnetization.

For antiferromagnetic interactions, the quantity v keeps going up as the
temperature is lowered. For two sublattice antiferromagnets, like the honey-
comb and square net, the local magnetic field distribution function is given
by the symmetric part (see Eq. 22) of the distribution for the corresponding
ferromagnet, (i.e. all exchange interactions changing sign). Thus v increases
monotonically from 1 at high temperatures to z at zero temperature. For the
triangular net antiferromagnet v increases from 1 at high temperatures to
1.097809 at zero temperature.13

Finally in this section. we consider the neutron scattering law (26) and
(30). The scattering only takes place at discrete energies. In Table 2 we
have included the thermal factor in (30) to give S(k,w) for the square net.

At very low temperatures, only the "spin wave peak" is seen at w 8 8.10 This

corresponds to flipping a spin where all its nearest neighbors are parallel.

As the temperature is raised, other peaks have nonzero weights.

'1!‘

106

Even at Tc’ we note that 63% of the weight is still in the spin wave peak.
However there is some weight at w 8 4, which corresponds to a spin flip
where 3 nearest neighbors are up and l is down. The peak at w 8 0 will
combine with the elastic peak and the peaks at negative frequencies are
related to those at positive frequencies by detailed balance factors. At
infinite temperature the thermal factor is 8 for all frequencies and the
neutron scattering law becomes symmetric. There is a useful sum rule for
spin 8 (additional to (31) and (32)) which shows that the total integrated

intensity is independent of temperature, as is clear from Table 2;

J 302,11) dw =- 1 . ' (55)

This is easily proved from the definition (26).
IV. CONCLUSIONS

We have shown that the statiscical mechanics of Ising models can be
described through the local magnetic field probability distribution function
P(h). This function determines both the neutron scattering law S(k,n) and
the thermodynamic quantities of interest for a large class of Ising systems.
We have calculated P(h) for the ferromagnetic honeycomb, square net and tri—
angular net and have shown that in all cases a pronounced dip develops at TC.

In a subsequent paper, we plan to extend these results to random spin
systems and spin-glass models and also to look at systems described by

classical rather than Ising spins.

 

107

ACKNOWLEDGMENTS

We should like to thank J. H. Barry, R. Haydock, T. A. Kaplan and

J. Stephenson for useful discussions and the N.S.F. and S.E.R.C. for

financial support. One of us (M.F.T.) would like to thank the Cavendish

 

Laboratory for its hospitality.

 

108

APPENDIX

The results of Sec. II can be generalized to case of arbitrary spins
Si whose magnitude can vary from site to site. We will adopt a numbering
scheme for equations such that equations in the text and the Appendix corre-

spond. The Hamiltonian is

aa-éf J..S?S?-)H.S‘T‘ . (Al)
- .. ij i j 1 i i
ij 1

Note that this J;. differs from the jij in the main text by a factor 4 in the

limit when the spin becomes 3. This is because we found it more convenient
to use Pauli operators in the main text. In a similar way there is a factor
2 different in the definition of Hi' Rather than redefine these quantities
we prefer to start from (Al). The reader will easily discover places where
factors of 2 appear between equations in the main text and the corresponding

equations in the Appendix. In order to proceed it is necessary to insist

that the Jii 8 O. The most obvious decomposition of (Al) is

h S? ~ (A2)

Pl‘l
H.
P

where

J 57‘ + 3. (A3)

(_J.‘ ‘v'l

but the most useful for our purposes is

H = -hi 3",: + H' (A4)
where

J.. S§+H. . (A5)

P
LJ. p4
"l
L:

H

109

. . . . . z .
Note that it is important that H be linear in all the 31 so that terms like

2 2 z z z 2 . . . . .
i) S , (SiSj) etc. are not allowed in the damiltonian if the

J

present formalism is to be used.

2 2
(Si) . (s

z . . . .
The thermal average < O 51.- > where O is any operator not involv1ng Slte

i is given by

z — z 2 g 3 z
< o si > =- < Tjr Lo 31 exp(8hiSi)] / i: [exp(ph5i)] >

 

8 < O 3- (,Sh.)> , (A6)
a. 1
1
- . 26 . - .
where the modified Brillouin :unction is defined by
3()=(S+l)rh’f3+1=)7 3 n’/"‘
S x 2 tot L, 2 x j 2 cot (x _, ,
and hence
B ( l h /2 \
3 x) 2 tan x .
Putting 0 equal to the unit operator in (A6) we find that
Z a
. B sh.)> ; A7
<si> < s,( 1 ( )
1
putting 0 equal to h; leads to
<h.Sz.>=- <EJ..s‘?sf+a.sf>
i 1 J ij i 1 i
a ' B c > ,
<ni S (phi) (A8)

110

There is the same problem with the factor 8 in obtaining the energy as

before and we have

Ei ' - < hi 83. (Bhi) >
1.
l

The magnetization M and the energy E are given by (10) as before.
It is convenient to define a P(h) for the whole system as

- fl

P(h) aL 33 (3111 91m) 1 r . 3 (3h) . (All)

1 1 1 1
and, further, an average Brillouin function 3(Bh) by

4 l

Mam-E 233 (an) .

i i

This takes account of a varying spin magnitude. It is only necessary to know
the number of sites with each spin magnitude in order to know B(Bh).

From (7) and (9) we have

7- -9,
< 51 > } 3 (Sb) Pi(h) dh (A15)

1

S

and

F: a _ i f, H’ 3 - 1'. z
“‘1 2 } h psiflsh) Pi(h) dh 2 Hi < si > . (A16)

The magnetization M is given by

M = N f 3 (8h) P(h) dh (Al8)

111

and,as before,it is not necessary to know the Ji or the Hi in order to find

1
M from P(h).

Again as before, the energy cannot be obtained directly from P(h). How-

ever if all the H are known to be zero, then from (10) and (A16)

1

E 8 -‘% t h 3 (8h) P(h) dh , (A21)

the free energy F is given by

1“

SF . - : m(lsi + 1) - ;' :de , (A14)

l J
O

and equation (ll) generalizes in the obvious way with tanh 3h - 3(8h).

The inelastic neutron scattering cross section is given by

2

+m
502,11) . 7 [

-co

where x is any direction perpendicular to 2. Using arguments similar to

those in Sec. II, we find that

* l ' -iwt ~ + - - +
S(k,w) - _— I dt » < s. s (c) + s. s. (c) >
Zn i i i i i
=- 2 j [ 9,.(11) + Pi(-'u) J 33 (8m) / [ 1 - exp(-8w) ]
1 " 1

= 2N E P(w) + P(-w) ] 3(Bw) / [ 1 - exp(-.6w) 1 (A29)

112
or

S(k,m) 8 4NPs(w) 3(Bw) / E l - exp(-3w) ]

These lead to the sum rules

 

(A30)
_l_ v . l - exp(-8m) .
2N )1 S(k,.u) 2 30311)) d 1 (A31)
and if all the Hi 8 O
1 — _ .
-§ S(k,n) L 1 - exp(-dw) j a dw . 5 (A32)
These sum rules are very general and can be proved without introducing
P(h). The second one is particularly easy to establish as the magnitude of the
spin does not enter explicitly.
Using the identity‘

< A(O) B(t) > 8 < B(t - 18) A(0) >

for any operators A,B; we see that

3(k,o) : 1 - exp(-ém) I

113

Only the diagonal terms contribute for Ising models so that

 

 

+<n
‘5, ._l_ {f 81.001: + -
S(k,m)[l exp(8w)] 2nje dc§<[si,si(c)]
- + 1
+ t ’ * .
Therefore
SCI”); i: 1 - exp(-dm) 3 add)
J
- +-
.. + 33‘. .78.
a - .‘ " a J. + " - b 1
l “T < L Di 9 5C 1 le 9 3C .1 >
i
,+ - - +
E < 1 >1 . his1 J E 5i , hiSi ] >
=- 4 ,‘7 < 11.5.2 >
i i i

\J
O

8.

10.

16.

17.

114

REFERENCES

W. Marshall, Phys. Rev. 118, 1520 (1960).

H- B. Callen, Phys. Lett. 4, 161 (1963).

B- Mfihlschlegel and H. Zittartz, Zeit. fflr Physik 122, 553 (1963).

M. W. Klein and R. Brout, Phys. Rev. 122, 2412 (1963).

B. W. Southern, J. Phys. C‘g, 4011 (1976).

D. Sherrington, A.I.P. Conf. Proc. 22, 224 (1975) reporting work by

S. Kirkpatrick, L. R. Walker and R. H. Walstedt, Phys. Rev. Lett. 28,

314 (1977); Phys. Rev. 3 22, 3816 (1980); K. Binder, Z. ffir Physik 3 26,
339 (1977); R. . Palmer and C. .. Pond, J. Phys. P‘g, 1451 (1979):

P. T. Bantilan and R. G. Palmer, J. Phys. P 11, 261 (1981).

S. N. Khanna and D. Sherrington, Solid State Comm. 36, 653 (1980).

J. R. McLenaghan and D. Sherrington, to be published.

T. C. Choy, D. Sherrington, M. Thomsen and M. F. Thorpe, to be published.
W. Marshall.and S. W. Lovesey, "Theory of Thermal Neutron Scattering"
(Oxford University Press, London, 1971).

See for example, L. Schiff in "Quantum Mechanics" (McGraw-Hill, NY, 1968).

L. Onsager, Phys. Rev. 62, 117 (1944).

T. C. Choy and D. Sherrington, J. Phys. A 16, L265 (1983» J. Phys.111§, 3691 (1983)-

See for example, M. F. Thorpe in "Excitations in Disordered Systems,"
N.A.T.0. Advanced Study Institute Series 878, edited by M. F. Thorpe
(Plenum, NY, 1982), p. 91.

These kinds of identities were first considered by M. E. Fisher, J. Math.
Phys. 4, 124 (1963).

M. F. Thorpe and A. R. McGurn, Phys. Rev. B fig, 2142 (1979).

P. H. Stillinger, Phys. Rev. 131, 2027 (1963).

"t

18.

19.

20.

26.

27.

115

J. Stephenson, J. Math. Phys. 1, 1123 (1966).
These results can also be derived by using the four equations:

[ P(h) dh . 1 , f h tanh 8h P(h) dh . 3s ,

f h P(h) dh - 3m and f tanh 8h P(h) dh - m .
See for example, L. Syozi in "Phase Transitions and Critical Phenomena"
edited by C. Domb and M. S. Green (Academic Press, NY, 1972) Vol. 1, p° 269.
See for example, 3. Stanley in "Introduction to Phase Transitions and
Critical Phenomena” (Clarendon Press, Oxford 1971).

J. Stephenson, J. Math. Phys. 2, 1009 (1964).

T. C. Choy (to be published).

R. B. Potts, Phys. Rev. §§_392 (1952) .

These involve using the star-triangle transformation on the correlation

functions, see J. H. Barry, C. H. Mfihera and T. Tanaka, Physica 1111’
367 (1982).

See for example, C. Kittel ”Introduction to Solid State Physics"

(J. Wiley, NY, 1971), p. 506. This modified Brillouin function BS(x)

differs from that usually used, which is

s i2§_:_£l 3§_:;11. _ .1. .JL
Bs(x) ZS coth [ 28 x ] ZS coth ( ZS ) ,

by including a prefactor of S, and a factor S in both arguments so that
85(x) 8 S-Bs(Sx). This saves having to write a lot of factors Si in
subsequent formulae.

D. N. Zubarev, Soviet Phys. Usp. 3, 320 (1960)..

116

Table Captions

 

Table l
The coefficients a:s defined in Eq. (39) are given for z 8 2, 3,

4, and 6.

Table 2
o o n o a n ’+ a
The inelastic neutron scattering intenSity per site S(k,m)/N at various
special temperatures for the square net. The frequency is in units of J and
the numbers in the table give the weights 2Wg/[l + exp(-3m)] of the five
3

delta functions.

117

_ _2 _
02 .t m2
2 m2 _2
22 o s
2 m T
c2 .7 -
_ _ .
o m s
dl
m
MI
1
m

T. —

o N

n d:

c #1

ma .1

c N

_ —

-2n121 212w! -
_ _2 _
_2 _2 m
_s _ m
_ _ _

 

 

_ c..-
e o2
m. N- a
SN 3 .VI
2 N c
e .V on
_ o d
-s - :72..- .2...
o u u
.2
.2
_
a
1.7..- ,-
. n n

M233.

 

c

 

 

118

 

d\T T 8 0 T 8 Tc T 8 m
8 1 0.6324 0.0625
4 0 0.2254 0.25
0 0 0.0849 0.375

-4 0 0.0387 0.25

-8. 0 0.0186 0.0625

 

Fig.

Fig.

H

rig.

Fig.

119

Eigure Captions
For the honeycomb lattice, the bar graph at the left and the line
graph in the upper right show the behavior of P(h) as a function
of temperature. The ordinate gives the weight wh in the delta
functions. The table in the lower right gives the value of wh
(defined in Eq. 37) at special temperatures.
Same as Fig. 1 except for the square net. The parameter x 8 l/w - l/rz.
Same as Fig. 1 except for the triangular net. The results for T < Tc
are not complete as we have not been able to compute the required
odd spin correlation functions.
The fluctuation v in the local field as a function of reduced

temperature T/Tc for ferromagnetic interactions in the honeycomb,

square net and triangular net.

120

a auawqm

 

 

 

 

 

 

 

 

;
m _ .. m-
8.8.2. £315 0 m-.. s m m s o
men-one» 38.2. o .-., . t r .00
memoirs. 38.8. o .3 81¢:
8.8.9. memoirs _ n; - . s -- o
8; up; one m m A m
. .00
minute
w S S w o
a a
o _ .
is . no
. a
. «0
V. «N l O
a
m no
\ -
m mousse.
A.
W O
m
A .8
m const»
A
\.

 

 

 

 

 

 

 

 

121

 

a yuan—m

 

 

 

 

 

 

 

 

 

 

..
s N o N- s-
88.8.3. 853.0183. o ...z s m a! o
8.9!. $828.3-.. 0 «-3 x m
men-onus 33881-6 o o; . o .00
3.0.3. $82.01..-. 0 N3 8"»;
ammo-0.3. mmvnmmdumkn . .i
Isl: . . o
8i are c... s s
. .80
o
-00
2...
o
.00
. o
...,-o

 

 

 

 

 

 

 

 

 

 

 

8090.032. $005.0 0 P;
2800.an 0003.0 0 .7;
338.032.... 03.2.00 0 «-3
3.00.05 08000.0 0 0;
233.0303. 0330.0 0 «3
05000.33 0003.0 0 3..
8090.049. «0003.0 . 03
8A or; 0.»
m a»:
1 x. N. 0.0 0.0
. 0
.. N0
2 - «.0
I mwuflfl_b lull.ll i mw.flv
weanlil
Nun: ..... .
I Our-Ill]. 1 mo
_ _ _

 

 

 

« aupaqm

 

 

 

Busts
“812“. _m an m. _m _u
. 1
m. use:
crush -
a L
0.0 use P

 

 

 

mAu

mAV

mgu

szn

0A0

nAu

123

q 323 _.

uh:
ca 0..

 

d d

bwz
mduaozfimh Ill

...mz wm<DOm Ilill
mioo>wzor ...........

 

 

 

 

 

 

>

Appendix B. Some 2D Correlations

In this appendix, some details are given for the calculation
of correlation functions needed in Chapter 5. First the calculation
of a four spin correlation on the square net will be discussed
and then the second part of this appendix will discuss the one honey—
comb correlation needed.

As mentioned in Chapter 3, most of the square net correlations
needed are found in Montroll, Potts, and Ward (1963), the exception
being the four spin correlation. Rather than rederive or rewrite
what has already appeared in MPW, their notation will be adopted here

and then enough information will be presented so that ig_conjunction

 

with MPW, the interested reader can calculate this particular
correlation.

MPW number the spins using rectangular coordinates so that
in their notation, what is needed is

ca 5 <012010021001> ' (1)

It is useful to rewrite this as

C4 = <012011011010021011011001> ' (2)

In analogy to MPW equation (49), one obtains the form

c = 24 P(Y‘1+Q) P(Y) (3)

4
where
z = tanh(BJ) , (4)

P(A) stands for the Pfaffian of matrix A, and Y and Q are matrices

124

125

which are defined below. It can easily be shown, as in MPW equation
(50), that

P(Y) = x4 (5)
where

x=z -z . (6)

The matrix (Y-1+Q) can be written down almost directly. The
columns and rows of the matrix are labelled in MPW's notation to
aid the reader in seeing how this particular form came about:
<Y‘1+Q) = <7)
(0,1)U (1,1)D (1,1)U (2,1)D (1,0)R (1,1)L (1,1)R (1,2)L

(0,1)U O A O B C D -D -C
(1,1)D -A O F O -D E E -D
(1,1)U O -F O A -D -E E D
(2,1)D -B O -A O C -D -D C
(1,0)R -C D D -C O A O B
(1,1)L -D -E E D -A O F O
(1,1)R D -E -E D O -F O A
(1,2)L C D -D -C -B O -A O

The six different entries are defined by:

A = (1-322)F10 - 2(1+22)FOO - x‘1 (8)
B = (1-22)F20 - z(1+zz)FlO - 222Fll (9)
c = ZZFIO + ZZFOO — F11 (10)
D = (1-22)F10 — zFll - zFOO (11)
E = -F + 22F + 22F (12)

00 10 ll

126

2 2 2
F = z(1+z )F10 + 22 Fll - (l-z )FOO (13)

The Fij's depend on elliptic integrals and are defined in the
appendix of MPW. In practice they are calculated numerically before
the Pfaffian is taken. Since the Pfaffian is the square root of

the determinant, the desired correlation is

c4 = 24x4[DET<Y‘1+Q>J

_1_
2

(14)

Turning to the honeycomb lattice, one starts with the partition
function (Green and Hurst 1964)

JN°Z' (15)

Z = C°[cosh3(BJ)
or .

Zn Z = In c + 3N2n cosh(BJ) + Zn 2' (16)

where c is a constant and

2N 2n 4
J d6 I d¢ Rn{1 + 3y

in Z' =-;%z
O 0

8H-
2 2
- 2y (l-y )[cose + cosm + cos(6-¢)]} ; (17)
y E tanh(BJ) . (18)
Noting that, unlike the convention used in Chapter 5, Green and

Hurst double count their bonds,

 

 

1 Bin 2'
<Oool> — tanh(BJ) +-§ 3(BJ)
'
= tanh(BJ) +-l sech2(BJ) §&2_§_ . (19)
3 By
One can show
Bin 2' _ -l-J2fld8 [2"d¢ A + Bcose (20)
By 81T2 O O C + Dcose + EsinB

with

127

 

A = 12y3 (21)
B = -3<4y - 8y3) (22)
C = (1+3y4) - 2y2(1-y2)cos¢ (23)
D = -2y2(1—y2)(1+cos¢) (24)
E = -2y2(1-y2)sin¢ - . (25)
Stephenson (1964) pointed out that
TT
_L_ exp(in8) _ n 2 _ 2 _ 2 -% .
2n [_nde c + Dcose + EsinB “ a (C D E ) ' (26)
_1_
a 3 [(c2 - D2 - E2)2 - C1/(D - 1E) . (27)

This relation is used to perform one of the integrals in equation

(20), giving

, ZTT _1_
an 2 L] de [A + BCD ] (CZ—DZ—E2)-2

 

By 4n 0 D2+E2
2n
+A%;-[ d6 ED 2 (28)
O D +E

Now <OOOl> is expressed (through equations (28) and (19)) in terms
of a pair of single integrals, a more suitable form for obtaining

numerical results.

Appendix C. Details of Monte Carlo Simulations

The basis for the Monte Carlo simulations of the spin glass
and ferromagnetic systems was the algorithm introduced by Metropolis
et. al. (1953). Given a system at a temperature T, the algorithm
generates an ensemble of states {Si}, which can be used to deter-
mine the thermal average of properties of the system. From the
ith member of an ensemble, Si’ the next member is generated by
creating a new state St. In Ising systems this is typically done

by flipping one spin in Si' Given

Et = energy of St (1)
E1 = energy of Si (2)
AB = Et - Ei , (3)
Si+l is defined in the following way:
= st

1) AE < O + Si+1

2) AE > 0

Choose a random number p such that O < p < 1.

t
a) p < exp(—AE/kBT) + Si+l = S

b) p > eXp(-AE/kBT) + 81+ S

1 i .
The initial state, 81’ may be chosen at random or taken to

be the last state generated by a simulation at a nearby temperature.

In either case, one generally excludes some of the initial members:

of the ensemble from the averaging process, which is equivalent

to allowing the system to come into equilibrium before making

128

129

a measurement.
In the simulations described below, a system change consists
of flipping a single spin. A Monte Carlo step is defined by N
attempted spin flips, where N is the number of spins in the system.
Most of the rest of this appendix will describe the spin
glass simulation, with some details of the ferromagnetic simulations
appearing at the end.
The Hamiltonian for the Ising Sherrington-Kirkpatrick model
is

H = -

NIH

Z J.,o.o. (4)
izj lJ l J

where as usual the superscript, z, is suppressed for convenience.
p . th
Let Oj— denote the value of OJ in the p state of the ensemble.

Suppose a trial state is arrived at by flipping the ith spin so

that its new value is O: = -05. Then

AB = -Z J..o?(oF - oi)

. l 1
J J J
= 2( E J..oP) o?
. lJ J l
- P P
‘ Zhi Oi (5)

At this point, the difficulties involved in doing a simulation
of the SK model are apparent. First one has to have stored in
an accessible way, the N(N-1)/2 independent values of Ji" Second
one must be able to calculate hi, which has N-l terms, efficiently.
Modifying an array packing scheme suggested by Jacobs and Rebbi
(1981), these difficulties can be overcome to a degree.

Neglecting the fact that the Jii's are zero, the information

130

describing this system is binary in nature:

J.. = t1
13
09::1

1

The Jii = 0 problem is easily treated separately in the computer
program. Rahter than taking an entire computer word, which is sixty
bits on a Cyber 750 system, to store one value of, say, 0:, one

can store sixty spin values in one word. These individual bits

can be accessed and modified through such FORTRAN commands as SHIFT,
OR, XOR, and AND. The last three correspond to the logical functions

whose operations are defined below.

AND(A,B) OR(A,B) XOR(A,B)
301 301 BAOl
000 001 001 (6)
101 111 110

The corresponding FORTRAN functions, operating on two words,
perform these logical operations on a bit by bit basis. For example,
AND(1,7) = 1 (7)
because 7 = 1112 and 1 = 0012. Only in the least significant bit
(LSB) do both words have a 1 so that 1 appears only in the LSB of
the answer. In fact, a projection operator, projecting out the
LSB of a word M, can be defined by
P1(M) = AND(1,M) . (8)
The SHIFT(M,Nb) function performs a cyclic permutation of the bits

Of word M. Nb bits to the left. Hence one can define a projection

operator for the Nth LSB of a word M by

131
PN(M) = AND(SHIFT(1,N-1),M) . (9)

Similar use of the OR function will allow the definition of a par-
ticular bit in a word. Hence each bit of the sixty bit word is
both definable and accessible.

The values of N spins can be stored in an array with N/60 words.
This is why N was taken to be a multiple of 60 for the spin glass
simulations. Similarly, one could store the N(N-1)/2 values of
Jij in an array of approximately N2/12O words. However, for ease
in calculation of the h? , the full matrix of Jij values was stored
in an array of N2/60 words. A simple rescaling of the energy allows
one to let a binary 0 represent —1 in the physical system.

The calculation of the hi is achieved by noting that when Jij
and 0? have the same value, 0? contributes positively to h: , while
negative contributions arise when Jij and 0? have opposite values.
Hence the product Jijog can be represented by the XOR function.
Furthermore, the FORTRAN XOR function will perform these multipli-
cations sixty at a time. What remains is to sum up over all j.
Apparently there is no easy way of doing this, and so each bit of
the product word is projected out separately and the results are
added as integers.

The memory limit of the Cyber 750, 127,000 words, is reached
for an SK system of about 2400 spins. The practical computation
limit, based on cpu time, is 2040 spins, while systems of 1020
spins were sufficiently large for most calculations. A typical
measurement at one temperature involved averaging over four 1020

spin systems; that is, four different Jij matrices were generated.

132

Each system was brought into equilibrium by running for sixty Monte
Carlo steps and then thermal averages were performed over two hundred
steps. Such a run requires about 37,000 words of storage and five
minutes of cpu time on the Cyber 750.

Not all-of the advantages to using binary methods for simulating
the infinite range spin glass model hold for simulating short range
ferromagnetic systems. First, for uniform nearest neighbor systems,
it is more efficient to devise a way of locating nearest neighbors
than to store the entire Jij array. Secondly, the number of multi-
plications involved in calculating 111 no longer scales with the system
size. However, with long range order being important in ferromag-
netic systems, it is useful to be able to work with larger systems,
particularly near Tc' Hence the amount of central memory associated
with storing the spin configurations becomes important.

There is a trade-off between efficiency of use of storage space
and ease with which spin values can be accessed. One compromise
between these two optimizations is to utilize the sixty bit length
of each word as an extra dimension. As an example, the spin con-
figuration of a 3OX3OX30 simple cubic lattice was stored in an array
of 3OX30 words, with each word in the array containing information
about 30 spins.

The chief difficulty in programming a Monte Carlo simulation
of these ferromagnetic systems is finding an efficient way to locate
the nearest neighbors of a spin. For a D dimensional simple cubic

system, the method is intuitive. A spin located at

(11.12.....1D> .

133

where iD refers to the bit number in word
(i1,i2,...iD_1) ,

has nearest neighbors located at
(ili1,iz,...,iD),
(il,i211,...,iD) , ...,.

A body centered cubic lattice, in three dimensions, can be
represented by two interpenetrating simple cubic lattices. Locations
on the BCC lattice are thus described by four integers:

(L,i1,i2,i3)
where L=1,2 is the sublattice number and

(i1,12,i3)
locates the position on that sublattice. With this decomposition
it is straightforward to write out a table which describes the
nearest neighbors of a given spin. Similarly, a face centered cubic
structure can be decomposed into four simple cubic sublattices and‘
a diamond structure into eight. The nearest neighbor locations
which are obtained in this decomposition are shown in Table 2.
The periodic boundary conditions used in these simulations are imposed
with the aid of modular arithmetic.

Typical runs on ferromagnetic systems began with 1000 Monte
Carlo steps to bring the system into equilibrium. This was followed
by 1000 steps (2000 steps at TC) during which thermal averages were

computed.

134

Table 2. Nearest neighbor locations for three dimensional lattices.

Site

(1,i,j,k)

(2.i.j.k)

Site

(l.i.j.k)

(2,i,j,k)

(3,i,j,k)

(4,i,j,k)

Site

(1:1!jak)
(2,i,j,k)
(3,i,j,k)
(4,i,j,k)
(5.1.j.k)
(6,i,j,k)
(7,i,j.k)
(8,i,j,k)

BCC
Nearest Neighbors

(2 .k)(211.Jk)(2111k)(2 j.k-1)
(2 ,j-—1,k) (2, i- -1 ,j, k- 1) (2, i, j- *1 fl) (2,i—1,j-1,k-1)

(1,1, j,k) (l, i+1, j, k) (1, i, j+1, k) (1, i, j,k+1)
(1, i+1, j+1, k) (1, i+l ,j, k+1) (L i, j+1, k+1) (1, 1+1 ,j+1, k+1)

FCC

Nearest Neighbors

(2,i,j,k) 2.i-l.j.k) (2.1.3.k-l) (2.1-1 j,k- l)
(Byisjyk) (39i-19jak) (3,i’j“lvk) (3’i-11j-12k)
(Avivjvk) (Arivj’l’k) (4,i,j9k-l) (4,i9j- “lik 1)

,i,j,k+1) (1,1
,i,j,k+l) (3,1
aivj-10k) (Aai

',j,k) (l,i+1,j,k) +1 ,j, k+1)
' ' '1,j-1,k+1)
+1 ,j— —1 ,k)

(l,1,j,k) (1,i+1,j,k) (1,i,j+l,k) (l,i+1,j+l,k)

(2,i,j,k) (2,i,j+19k) (Ztlijik-l) (2,i,j+19k-1)

(4,i,j,k) (401+19j0k) (4919J9k'1) (491+19j9k‘1)

(1,i,j,k) (1,i,j+l,k) (1,i,j,k+1) (1,1, j+L k+1)

(2,i,j,k) (2,i- 1, j,k) (2,i,j+l,k) (2.1-1, j+1, k)

(3,i,j,k) (3,i-1,j,k) (3,i,j,k+l) (3,i-1,j,k+1)
DIAMOND

Nearest Neighbors

(S,i,j-1,k) (6,i-1,j-1,k-1) (7,i-1,j-1,k) (8,i,j—l,k-1)
(5,i,j-1,k) (6,i,j-1,k) (7,i,j-1,k) (8,i,j-1,k)
(S,i,j,k) (6,i,j,k-1) (7,i,j-1,k) (8,i,j-1,k-1)
(S,i,j,k) (6,i-1,j,k) (7,i-1,j-1,k) (8,i,j-1,k)
(1,i,j+l,k) (2,i,j+l,k) (3,i,j,k) (4,i,j,k)
(l,i+1,j+1,k+1) (2,i,j+l,k) (3,i,j,k+1) (4,i+l,j,k)
(1,i+1,j+1,k) (2,i,j+l,k) (3,i,j+l,k) (4,i+1,j+1,k)
(1,i,j+l,k+l) (2,i,j+l,k) (3,i,j+1,k+1) (4,i,j+1,k)

LIST OF REFERENCES

LIST OF REFERENCES

Bantilan F T and Palmer R G 1981 J. Phys. F.ll 261

Barry J H, Mfinera C H, and Tanaka T 1982 Physica 113A 367
Berlin T H and Kac M 1952 Phys. Rev. §§_821

Bethe H A 1935 Proc. R. Soc. London Ser. A 15_ 552

Binder K 1977 Z. Physik B 26 339

Bonner J and Fisher M 1964 Phys. Rev. $355 640

Burington R S 1946 Handbook 9f_Mathematical Tables and Formulas
(Sandusky: Handbook Publishers, Inc.)

Callen H B 1963 Phys. Lett. 4_161

Choy T C 1984 submitted to J. Math. Phys.

Choy T C and Sherrington D 1983a J. Phys. A‘lg L265
Choy T C and Sherrington D 1983b J. Phys. A.16 3691

Choy T C, Sherrington D, Thomsen M, and Thorpe M F 1984
to be published

Domb C 1960 Adv. Phys. 2_251

Domb C 1974 "Ising Model" in Phase Transitions and Critical Phenomena
vol. 3, C Domb and M S Green eds. (New York: Academic Press)

Edwards S F and Anderson P W 1975 J. Phys. F.§ 965

Eggarter T P 1974 Phys. Rev. Big 2989

Fisher M E 1963 J. Math. Phys. 4 124

Frank B, Cheung C Y, and Mouritsen O G 1982 J. Phys. C.14 1233
Gaunt D S, Sykes M F , and McKenzie S 1979 J. Phys. A l2_871

Green H S and Hurst C A 1964 Order-Disorder Phenomena (New York:
John Wiley and Sons)

 

135

136

Grest G S and Rajagopal A K 1974 J. Math. Phys. l§_589
Hurst C A 1963 J. Chem. Phys. 38 2558

Ising E 1925 Z. Physik §1_253

Jacobs L and Rebbi C 1981 J. Comp. Phys. 4; 203

Katsura S 1962 Phys. Rev. 121 1508

Khanna S N and Sherrington D 1980 Solid State Comm. 36 653
Kirkpatrick S and Sherrington D 1978 Phys. Rev. B 11 4384

Kirkpatrick S, Walker L R, and Walstedt R E 1977 Phys. Rev. Lett.
38 514 ‘

Kittel C 1971 Introduction £2_Solid State Physics (New York:
John Wiley and Sons)

 

 

Klein M W and Brout R 1963 Phys. Rev. 132 2412

Kosterlitz J M, Thouless D J, and Jones R C 1976 Phys. Rev. Lett.
36 1217

Kramers H A and Wannier G H 1941 Phys. Rev. §Q_252
Krey U 1981 Z. Physik 3.3g 231

Lee D H, Joannopoulos J D, and Negele J W 1984 Phys. Rev. Lett.
52 433

Lieb E, Schultz T, and Mattis D 1961 Annals of Physics l§_466
Marshall W 1960 Phys. Rev. 118 1520

Marshall W and Lovesey S W 1971 Theory 2: Thermal Neutron Scattering
(Oxford: Oxford University Press)

 

Matsubara F and Katsura S 1973 Prog. Theor. Phys. 49 367

Matsubara F and Yoshimura K 1973 Prog. Theor. Phys. §Q_1824
Matsubara F, Yoshimura K, and Katsura S 1973 J. Can. Phys. §l_1053
McLenaghan J R and Sherrington D 1984 to be published

Mess K W, Lagendijk E, Curtis D A, and Huiskamp W J 1967
Physica.§4_126

137

Metropolis N, Rosenbluth A W, Rosenbluth M N, Teller A H, and
Teller E 1953 J. Chem. Phys. 21 1087

Mézard M, Parisi G, Sourlas N, Toulouse G, and Virasoro M 1984
Phys. Rev. Lett. §2_1156

Montroll E W, Potts R B, and Ward J C 1963 J. Math. Phys. 4 308
MUhlschlegel and Zittartz H 1963 Z. Physik 11; 553
Naya S 1954 Prog. Theor. Phys. _1 53

Onsager L 1944 Phys. Rev. §§_117

Palmer R G and Pond C M 1979 J. Phys. F 2 1451
Parisi G 1980a J. Phys. A 1§_L115

Parisi G 1980b J. Phys. A.13 1101

Parisi G 1980c J. Phys. A.£3 1887

Parisi G 1983 Phys. Rev. Lett. 29 1946

Peierls R E 1936 Proc. Cam. Phil. Soc. §2_477
Potts R B 1952 Phys. Rev. §§_392

Schiff L 1968 Quantum Mechanics (New York: McGraw—Hill)

 

Sherrington D 1975 A.I.P. Conf. Proc. 22 224
Sherrington D and Kirkpatrick S 1975 Phys. Rev. Lett. 35 1792
Southern B W 1976 J. Phys. C 2_4011

Stanley H E 1971 Introduction to Phase Transitions and Critical
Phenomena (Oxford: Oxford University Press)

Stephenson J 1964 J. Math. Phys. §_1009

Stephenson J 1966 J. Math. Phys..Z 1123

Stillinger F H 1963 Phys Rev 131 2027

Suzuki M, Tsujiyama B, and Katsura S 1967 J. Math. Phys. §_124

Sykes M F 1979 J. Phys. A 12_879

138
Syozi I 1972 in Phase Transitions and Critical Phenomena vol. 1
C Domb and M S Green,eds. (New York: Academic Press)
Takeda K, Matsukawa S, and Haseda T 1971 J. Phys. Soc. Jap. 39 1330
Thomsen M and Thorpe M F 1983 J. Phys. C.1§ 4191

Thomsen M, Thorpe M F, Choy T C, and Sherrington D 1984 to be
published (in Phys. Rev. B)

Thorpe M F 1982 in Excitations in_Disordered Systems NATO Advanced
Study Institute Series B78, M F Thorpe,ed. (New York: Plenum)

 

 

Thorpe M F and McGurn A R 1979 Phys. Rev. 3.29 2142

Thorpe M F and Miyazima S 1981 Phys. Rev. B 24 6686

Thorpe M F and Thomsen M 1983 J. Phys. C 1§_L237

Thouless D J , Anderson P W, and Palmer R G 1977 Phil. Mag. 35 593
Walker L R and Walstedt R E 1980 Phys Rev B 22 3816

White R M 1983 Quantum Theory of Magnetism (New York: Springer-
Verlag)

Wright J C, Moos H W. Colwell J H, Mangum B W, and Thornton D D
1971 Phys Rev B.3 843

Zubarev D N 1960 Soviet Phys. Usp. 320