A STUDY OF THE EFFECT or CRYSTAL-FIELD ANISOTROPY 0N MAGNETICALLY ORDERED SYSTEMS ’ Thesis for the Degree of Ph. D. MICHIGAN STATE UNIVERSITY JOHN F. DEVLLN 1970’ "-w‘ > »- w.....—.nu~q. \HES‘S‘ 3 1293 10108 8031 . 5 L I B R 11 R Y Nlizij“ ' if? '3. State Unwasity This is to certify that the thesis entitled A STUDY OF THE EFFECT OF CRYSTAL-FIELD ANISOTROPY ON MAGNETICALLY ORDERED SYSTEMS p eeeee ted by John F. Devlin has been accepted towards fulfillment of the requirements for _Eh_._D_.__. degree in Miss Bglflfi‘l ‘D. Sf‘fl" F Major professor Date 10/15/70 ABSTRACT A STUDY OF THE EFFECT OF CRYSTAL-FIELD ANISOTROPY ON MAGNETICALLY ORDERED SYSTEMS BY John F. Devlin A study has been made of the effect of crystal- field anisotropy on magnetically ordered systems by using the Green function techniques. The anisotrOpy was repre- sented by a term D(Sz)2 added to the usual Heisenberg Hamiltonian. For a spin S system we have found a convenient set of Operators, ZS in number, which are appropriate for completely solving, within the Green function form- alism, a system of non-interacting spins in a crystal- line field environment.'The inclusion of exchange coupling between the spins introduced higher order Green func- tions. These higher order Green functions were found to be only higher order exchange Green functions, and not higher order anisotropy Green functions. Therefore there was no need to approximate any of the anisotrOpy functions. We have found an approximation scheme which reduces the higher order Green functions to lower order John F. Devlin ones, thereby decoupling the set of equations. This scheme is a generalization of the random phase approxi- mation (RPA). For D=0, our set of equations reduces to the set obtained by Tahir-Kheli and ter Haar when they used the random phase approximation. Our tech- nique is valid for any size spin system, and also for antiferromagnetic as well as ferromagnetic systems. We have numerically solved the necessary set of equations for the spin S=l, ferromagnetic and antiferro- magnetic systems with simple cubic and body centered lat- tices, and for the spin S=3/2 ferromagnetic system with a simple cubic lattice. The decoupling approximation we have presented here is superior to those that were used in the earlier work on this subject. The transi- tion temperature Tc(D) as a function of the crystal- field parameter D is finite for all values of D. Our particular decoupling approximation shows that TC(D) depends more strongly on D, for low anisotropy, than does the molecular field theory (MFT). For anisotrOpies of the order of the exchange, J, or smaller, our results are equivalent to those of Lines. At the transition temperature the ensemble averages of the different powers of S2 as predicted by our RPA calculation and as predicted by the MFT are the same. The temperature dependence of these ensemble averages, below Tc, shows a greater variation with changes in D/J than does the MFT prediction. For very large anisotropies our results John F. Devlin approach those of the MFT, and at D/J = w both theories agree exactly. This is to be contrasted with the earlier single excitation schemes which predict TC = m at D/J = w. A STUDY OF THE EFFECT OF CRYSTAL-FIELD ANISOTROPY ON MAGNETICALLY ORDERED SYSTEMS BY John F. Devlin A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics 1970 DEDICATION To my wife, Patricia, who spent four long years in a town she cared little about, listening to a subject she knew nothing about, waiting for me to get a decent job, I dedicate this thesis. ii ACKNOWLEDGMENTS My sincere thanks go to Professor Robert D. Spence, who first suggested the topic for this thesis, and who effectively guided me to its conclusion. I would also like to thank Dr. Jack H. Hetherington for providing me with a c0py of his computer routine ANQB, which per- formed all of the necessary numerical integrations. Finally, I would like to express my appreciation to the National Science Foundation for supporting me with the three-year Traineeship under which this work was performed. iii TABLE OF CONTENTS LIST OF FIGURES Chapter 1. INTRODUCTION 2. THE MOLECULAR FIELD THEORY 3. THE GREEN FUNCTION TECHNIQUE 4. THE SINGLE EXCITATION THEORIES 5. THE MANY-EXCITATION SOLUTION FOR A FERROMAGNETIC SYSTEM A. The Non-Interacting Spin System B. The Exchange Term Decoupling C. The Complete Solution D. The Transition Temperature E. A Decoupling Approximation for the Molecular Field Theory 6. THE MANY-EXCITATION THEORY FOR AN ANTI- FERROMAGNETIC SYSTEM 7. THE RESULTS FOR SOME SPECIAL CASES A. B. C. The Spin S=l Ferromagnet The Spin S=l Antiferromagnet The Spin S=3/2 Ferromagnet SUMMARY AND CONCLUSIONS iv Page vi 13 l7 17 22 3O 33 34 36 47 47 57 7O 76 Chapter Page REFERENCES 79 APPENDIX A 81 APPENDIX B 84 Figure 1. LIST OF FIGURES Page The ratio T (D)/T (0) versus D/J for a simple cugic spin S=l ferromagnet. (a) Narath's results. (b) Anderson and Callen's results. (c) Lines's results. (d) The results of this thesis. (e) The molecular field theory results. 51 The ratio TC(D)/TC(0) versus D/J for a body centered cubic spin S=l ferromag- net. (a) Narath's results. (b) Anderson and Callen's results. (c) Lines's results. (d) The results of this thesis. (e) The molecular field theory results. 53 The sublattice magnetization versus T/TC for a simple cubic spin S=l ferromagnet. Curves (a), (b), and (c) are the Green function calculations of this thesis. Curve (d) is thermolecular field theory for the same three values of D/J, but in this case the three lines are indis- tinguishable. 56 The ensemble average <(Sz)2> versus T/TC for a simple cubic spin S=l ferro- magnet. (a) The Green function calcula- tion of this thesis for D/J=l.0. (b) The molecular field theory prediction for D/J=l.0. Note that at T=Tc the curves cross each other. 59 The sublattice magnetization versus T/TC for a simple cubic spin S=l antiferromag- net. (a) D/J=0.0l. (b) D/J=O.l. (c) D/J=l.0. 67 The ensemble average <(SZ)2> versus T/TC for a simple cubic spin S=l antiferromag- net. (a) D/J=0.0l. (b) D/J=O.l. (c) D/J=l.0. 69 vi Figure Page The ratio TC(D)/T (0) versus D/J for a simple cubic Spin S=3/2 ferromagnet. (a) Narath's results. (b) Anderson and Callen's results. (c) Lines's results. (d) The results of this thesis. (e) The molecular field theory results. 75 vii l . INTRODUCTION A magnetically ordered system may be defined as a system in which the ensemble averages of the magnetic moment operators, Ei' or equivalently the spin operators, Si, associated with the lattice sites labeled by the index i, are non-zero, even when there is no external magnetic field. Since the ensemble average is obtained by taking the trace of the product of Si and the density Operator, exp (-H/KT), one first needs to know the Hamiltonian, H, of this system before attempting to calculate any of these statistical mechanical quantities. However, for such a system the true Hamiltonian is much too complicated to ever be of any direct use. One generally resorts to using a Hamiltonian which is appr0priate only for the spins of the system, and not for the whole crystal, and all its constituents. The spin Hamiltonian, as it is called, is rarely, if ever, derived from first principles; rather it is a construct of reasonable and convenient forms of quantum mechanical spin operators. The assumption most frequently made is that the spin Hamiltonian depends only linearly or bilinearly on the spin operators. In fact, for systems in which there is no external magnetic field one invokes time reversal invariance of the Hamiltonian to say that products of an odd number of spin operators cannot appear at all, since the Si change sign under time reversal, whereas H cannot. The allowed forms of the spin operators are generally of two distinct types: single ion terms and exchange terms. The single ion spin Hamiltonian is of the form, = -guBZ {Rafi a + Z 28Aa8%5 5“ 58 (1-1) a HSI where the index i is for the lattice sites, and a and 8 represent the components of the spin operator vector. The coefficients g and “B are the gyromagnetic ratio and the Bohr magneton respectively. The first term is the Zeeman term, and it is just the interaction of the spins with the external magnetic field, HO. The second term represents the splitting of the low lying energy levels of the system due to the presence of a crystalline elec- tric field around each ion(l). This term has the informa- tion of the symmetry about each ion included in the para- a8 meters Ai The interaction between the spins may be written as, _ _Z Z Z ZJSBS sis SB (102) i j a B 3 j GB The coupling coefficient J31] is commonly refered to as the exchange energy or exchange integral. A special case of the Hamiltonian is the much used Heisenberg-Dirac Hamiltonian, where we now assume Jij=J. and Jii=0. This is a completely ji isotropic Hamiltonian, i.e. the ath component of the total spin vector commutes with H for all a, whereas it does HD not necessarily do so for the Hamiltonian H As such, I. HHD presents a problem in trying to describe a magnetically ordered system. By the symmetry of H = Trace HD' a exp(-H/KT)) should be the same for all a. But this can (Si only be if =0 for all a. Hence in attempting to de- scribe any magnetically ordered system with the Hamiltonian HHD one must include terms which are anisotropic, i.e. terms which produce a preferred direction in space along which the spins may align. Two simple alternatives are possible at this point: we may include terms of the form, Jgjsgsg in which Jgj is different for different a, or else we may include terms like D:(S§)2 where the D? are different for different a. The first possibility is merely a more com- plicated form of exchange and has a weaker theoretical basis than even HHD itself. The second possibility is more reasonable since, as we have already mentioned, it may be present when the effects of crystalline field are taken into account. In this thesis we will be concerned only with a spin Hamiltonian that has cylindrical symmetry. The preferred direction will be the z axis. Then the apprOpriate Hamil- tonian for our problem, within the limits we have already set, is as follows: _ _ _ _ z 2 H — E ZJij~i sj ZD(Si) , (1.4) l j i where we have set the external magnetic field to zero since it will not be needed for any of our calculations. We have also assumed that the crystal field parameter D is the same at each lattice site. The physical effect of the crystal field anisotrOpy term is clear: spin states with large absolute values of the quantum number, ms, corresponding to the eigenvalue of the 2 component of the spin operator, will have the lower energies (assuming D Z O) and hence they will be more highly pOpulated than spin states with small Imsl values. This is completely equivalent to saying that the z direction is preferred over the x and y directions. Note that the Hamiltonian we have chosen makes no statements about pre- ferences between the x and y directions, hence ==O. Furthermore, the ensemble average of any function of S? must be the same if S? is replaced by SE. Once the spin Hamiltonian has been established for a given system, the next problem is to determine the transi- tion temperature, i.e. the temperature above which the system is no longer ordered when the external field is set to zero. In order to do this exactly for a many-body system, either one has to be able to write the Hamiltonian as a sum of one-particle Hamiltonians or else one has to diagonalize the Hamiltonian completely. Because of the presence of the exchange term in (1.4) it is not possible to write it as a sum of one-particle Hamiltonians. The second alternative is then to try to diagonalize H exactly. No one has ever done this except for some trivial one-dimensional lattices, and some N-particle lattices where N is a very small in- teger. Therefore one must rely on some sort of approxima- tion scheme for obtaining a reasonable estimate of transi- tion temperature. We shall discuss five previous attempts at obtaining this temperature from (1.4) before we present our own method. Each of these theories is derivable by appropriate "decoupling approximations" from the Green func- tion technique, which will be discussed in a later chapter. 2. THE MOLECULAR FIELD THEORY The molecular field theory (MFT) is by far the simp- lest of all the effective field theories to apply. It is clear from (1.4) that the Hamiltonian for a single spin Si is of the form, ~ H — z 2 i —{-2 gJij°§j}°§i'D(Si) . (2.1) The factor of 2 appears because we are no longer double counting over all lattice points as we did in (1.4). Hi would be a one-particle Hamiltonian if the quantity in braces were a real number rather than a quantum mechanical operator. This leads us to try and replace it by its ensemble average, )J..s. + {J... (2.2) j lJ~J j 13 3 This replacement is equivalent to saying that each spin interacts with an effective magnetic field, eff _ 2 H — 2§Jij/gu8, (2.3) :instead of with the individual spin operators Sj. Since (5:) is independent of the lattice site index i for a ferro- Inagnetic lattice, we denote S = . We also make the def- inition J(O) = EJ1j. Then the new Hamiltonian is, j _ _ - z _ z 2 Hi — ZS J(O)Si D(Si) . (2.4) 6 The eigenvalues of this Hamiltonian are, - 2 Em = -ZS J(0) m - Dm , (2.5) where the allowed values of the azimuthal quantum number m are, m = S, 8—1, 00., -S+l' “'So (206) Then the ensemble average of 5: becomes, = S = {2m exp[(2SJ(0)m + Dm2)/KT]}/ . m _ 2 (2.7) {Eexp[(25J(0)m + Dm )/KT]}. m For any arbitrary power of 8:, we have, <(S:)n> = {EmneXp[(2SJ(0)m + Dm2)/KT]}/ m (2.8) {Eexp[(2SJ(O)m + Dm2)/KT]}. m At the transition temperature S=0, since this is the temperature at which the magnetic ordering vanishes. By expanding the exponentials in the limit S+O, one may ob- tain the following expression for the transition temperature according to the MFT, KTC = 2J(0)[fmzexp(DmZ/KTC)1/[Zexp(um2/KTC)1. (2.9) m m VNhere we have defined TC to be the transition temperature, (pr, more commonly, the Curie temperature. This expression can also be written as, _ z 2 KTC — 2J(0)<(Si) >|T=TC (2.10) Two special cases of (2.9) are as follows, 28(S+1)J(0)/3, (2.11) for D=0, T 282J(0). (2.12) for D: co, T Graphs of TC(D)/TC(O) versus D/J are plotted in Figures 1, 2, and 7 in Chapter 7. Tc(D) is the Curie temperature for finite D and TC(O) is for D=0. The exchange integral in this case is between nearest neighbors only. The var- iable J(O) which appears in all the equations then reduces to J(O)=zJ, where z is the number of nearest neighbors. (z=6 for a simple cubic lattice, and z=8 for a body centered cubic lattice.) In chapter 5, we shall again present the MFT, but this time as a decoupling approximation within the Green func- tion formalism. 3. THE GREEN FUNCTION TECHNIQUE The method of double-time temperature-dependent Green functions was first applied to the problem of ferromagnet- ism by Bogolyubov and Tyablikov(2)'(3). Since that time a large number of problems in ferromagnetism, antiferromag- netism, and ferrimagnetism have been solved by this tech- nique(4)-(9). We will present here only a very brief (10) outline of the method. (See Zubarev for complete de- tails and example problems.) The retarded double-time temperature-dependent Green function is defined as, <>= -i9(t-t')<[A(t) ,B(t')]>, (3.1) where 6(t-t') is the unit step function, singly pointed brackets denoting the ensemble average as usual, and doubly pointed brackets denoting the Green function. The equation of motion of this Green function is, igE<>= 5(t-t')<[A(t),B(t)]?J~~* (3.2) +<<[A(t),H];B(tl)>>. FTNarier transforming over the time variable t into the energy Variable E (h=l) gives, lO E<> = E <[A,B]> +<< [A.H1;B>> (3.3) E! NH :3 where the E subscript on the double brackets denotes the Fourier transformed Green function. Once equation (3.3) has been solved for <>E, the double-time correlation function can be written down quickly with the aid of the Spectral Theorem (10): . +m << >> . -<< = ile A'B E=w+1€ A'B>>E=w-ie 5+0 -w exp(w/KT)-l (3.4) X exp(-iw(t-t')) do. There are two main difficulties to be overcome for every statistical problem to be solved by the Green function method: (1) The operators A and B have to be chosen very carefully, as not all choices are convenient and useful. (2) Once A and B are chosen, one needs to solve the equa- tion (3.3) for <> This last step is a non-trivial E' task since the Green function <<[A,H];B>> is also an unknown. If the equation for <<[A,H];B>> is written down one finds that it, in turn, depends on the Green function <<[[A,H],H];B>> which is also an unknown. Hence there is really a whole hierarchy of Green function equations implied by (3.3). For any given system with a many-body Hamiltonian it is impossible to solve the complete hierarchy of equations exactly. Approximate solutions to (3.3) may be obtained only by decoupling the chain of equations so as to reduce it to a much smaller number of equations. In general these decoupling approximations all assume the form of 11 writing a selected higher order Green function as a linear combination of lower order Green functions. Then all the Green functions which are of higher order than the one selected are ignored completely. The problem of decoup- ling is thus reduced to deciding at which point the hierarchy is to be truncated, and to decide just what the coefficients in the linear expansion should be. Since the chosen higher order Green function is writ- ten as a linear combination of the lower order ones, this finite set of equations will be a set of linear equations. The solutions for any one Green function in the set of N equations will then be of the form of a polynomial of de— gree N-l in the variable E, divided by a polynomial of degree N in E. A partial fraction expansion will then yield a Green function which has, at most, N simple poles in the complex B plane. For a real physical problem these poles will always lie on the real axis. The limit-taking pro- cess inside the integrand of (3.4) is easily performed using the identity, Lim l 1 _ _ . _ 8+0 {E“Eo+i€ - E-EO-ig} — 2N16(E E0). (305) The following corollary to the Spectral Theorem will be of use in later calculations: when the Green function <>E has the form, A1 E-E. l (3.6) 1 << 2—... A'B>>E 2n g 12 where the Ei are the poles, and the Ai their residues, the equal-time correlation function of B and A has the form, A. 1 expTEi/KT)-l' (3’7) = = Z i 4. THE SINGLE EXCITATION THEORIES All of the early attempts at solving the ferromagnetic system described by the Hamiltonian (1.4) were single ex- citation theories, i.e. their first order Green functions were approximated by a single simple pole on the real axis. This is the most elementary of all approximation schemes. The choices for the operators were: A=S; and B=Bh, where g and h denote lattice sites (not necessarily the same) and Bh is some function of s; and 8:. The com- plete equation of motion of the Green function is obtained by working out the commutator of S; and H, and then insert- ing this into (3.3). The result is, + z + o -' . ° > > E<> = F6 2n + 2 J. << 8 s. g.h/ E 19 ( g g (4.1) + D<<(s+s 828+);B >>, 9 9 9 2- g h where F=<[S;,Bg]>, and 69 h is the Kronecker delta function. I For convenience we have also dropped the subscript E on the Green function. Note that both of the Green functions on the right hand side of (4.1) are of higher order. Two very effective decouplings of the exchange term (3) Green function have been made, one by Tyablikov , the other by Callen(5). The Tyablikov decoupling is more 13 l4 commonly referred to as the random phase approximation (RPA). These decouplings are as follows: Tyablikov Decoupling (RPA): + z _ z +. . <>-<> for 1#g, (4.2) Callen Decoupling: <>=<>- 1 g g 1 h {3‘ (4.3) -l—> for ifig, 282 1 g - + S.><> = 2<>, (4.4) 9 9 SS9 h 9 9 h Anderson and Callen(7): <<(s+sz+szs+);B >> = 2<> 9 9 9 9 h 9 9 h (4.5) - 1 <5 z><(s S++ S+S )><>, 252 9 9 9 9 9' h Lines‘ls): + 52+ szS + << 8 ;B >> = ¢ <> (4.6) ( 9 59+ SS9 9) h 9 9 h ' where, <3(sz)z43(s+1)> _ Q = 9 for B = S , (4.7) g <82> h h 9 and z 3 2 z 2 z <2(s; ) 4+3(Sg ) -(ZS +28—l)(Sg) -s(s+l)sg> ¢ = g <(sg 3)3- (s ;)2-S(s+1)s; > z - (4.8) for Bh=shsh' Plots of TC(D)/TC(O) versus D/J for each of these decouplings 16 are contained in Figures 1, 2 and 7. In each case the exchange term Green function was decoupled by the RPA. In the limit D/J+®, each of these single excitation decouplings predicts an infinite transition temperature. This is clearly an unphysical result. For extremely large anisotropy the effective spin states at each site become |S> and |-S>, and only transitions between these two states can occur. The intermediate spin states, such as IS-l> and |-S+1>, are simply not pOpulated at all because they are so widely separated in energy from the states IS> and |-S>. A system with only two effective spin states still has a finite, though larger, transition temperature than a system with some populated intermediate states. Therefore, the prediction of TC:00 is completely erroneous. These single excitation decouplings are, at best, valid only for very small anisotrOpy. The Narath and Anderson and Callen decouplings also contain another fallacy. Both predict <(S;)2>=S(S+l)/3 at, and above Tc‘ But this can only be true for an iso- tropic system which has D=0. Lines's results do not suf- fer from this defect for D20; he has <(S:)2> 2 S(S+l)/3 at and above Tc' 5. THE MANY-EXCITATION SOLUTION FOR A FERROMAGNETIC SYSTEM A. The Non-Interacting Spin Problem Once it has been established that the single excita- tion approximations to the crystal field Green function will not work, one is tempted to try and use the higher order Green functions generated by the equation of motion of the lower order Green functions. To illustrate how this might be done we will consider the case of a lattice of non-interacting spins, which has the Hamiltonian, H = ZHi’ (5.1) l where _ _ z _ z 2 Hi _ ngHosi D(Si) . (5.2) As we saw with the MFT Hamiltonian (2.4), this problem can be solved very easily. Hi is already in diagonal form, hence exp(-Hi/KT) is quickly written down and one can then find all the relevant ensemble averages for this system. For example, for a spin S=l system the results are, exp(guBHO/KT) - exp(-guBHO/KT) - exp(-D/KT)+exp(guBHO/KT)+ex§7_guBHO/KT)a (5.3) 17 18 (5.4) ?)2 eXP(gUBHO/KT)+exp(«guBHo/KT) (S1 = exp(-D/KT)+exp(guBHO7kT):exp(-guBHO/KT)' It is also enlightening to solve this problem by the Green function technique. That this can be done exactly within the Green function fOrmalism is almost obvious. The set of all possible linearly independent spin Operators asso- ciated with a lattice site i for a given size spin S, is a closed set under commutation with the Hamiltonian Hi in (5.2). Hence as we attempt to generate the higher order Green functions <<[A,Hi];S;>>, we must eventually return to the lower order Green functions. Therefore we will only have a finite number of operators to work with, and cor- respondingly only a finite number of Green functions. As an illustration of this, consider, again, a Spin S=l system. For the first Green function choose A=S:, and 3:87. Then, 1 + — - + - ,° , > = 1 << .; .>> E< ZS/Zn + gnBHO S1 S1 _ (5.5) + D<<(s?sl + sls?);s.>>. l l l l 1 Now writing down the equation of motion of the higher order Green function that appears above we have, 2 + + z — z + + z - . .+ . . ; .>> = + << . .+ . . ; .>> E<<(slsl 8181) S1 Z/2n guBHO (slsl slsl) 31 (5.6) + _ +D<>, where Z = <[(S?Sf + STS?),ST]>. 1 1 1 1 1 19 This second equation has closed the set of equations upon itself. The solutions for the Green functions in these equations are, + - 1 { s+z §-z ) <> = —— \ _ _ _ , (5-7) 1 1 2n LE ngHO D E gUBHo+DJ z + + z - _ l . S+Z <<(sisi+sisi),si>> ‘ 2n )E-gu H -D k 8 o (5.8) Then using the corollary to the Spectral Theorem (3.6) we have, = (S+Z)/(exp((guBHO+D)/KT)-l) _ (5.9) + (S-Z)/(exp((guBHO-D)/KT)-l) = (S+Z)/(exp((guBHO-D)/KT)-l) _ (5.10) + (S-Z)/(exp((ngHO+D)/KT)-l) But since, = 2-§-<(s:)2>, (5.11) z = 4 - 2<(s:)2>, (5.12) = 2+5 - <(s?)2>, (5.13) 11111 J. equation (5.9) and (5.10) are really two coupled equations in the two unknowns S, and <(S:)2>. If we solve for these unknowns, we obtain the expressions (5.3) and (5.4) just as we had hoped to do. The point of this alternative derivation 20 was to learn how to deal with crystal field anisotropy in the Green function technique. It has essentially told us how to choose the operators, and hence the Green func- tions, for any size spin system. For a spin 8 system we will need 25 different opera- tors for the Operator A in (3.3). The most convenient set is as follows: 1 + A = S 5.14 g g, ( > i _ i-l + '3 where the braces L} + denote the anticommutator. Though the set {A;} contains only ZS operators it is still closed under commutation with the Hamiltonian in (5.2), since the following commutation relations hold, 1 i z . A S = -A for all 5.16 [9,! g] g 1! ( ) [A1,(sz)2] =—A1+1 for all 1525-1, (5.17) 9 9 9 and 25-1 . [A:S,(S;)2] =-—) ajAg, (5.18) 3 where the coefficients aj depend on the particular size of the spin of the system. By defining the Green functions to be Gl(E)=<>,the equations Of motion are, (E-guBHO)Gi(E)-DGi+1(E) = zi/zn, for 1525-1, (5.19) 23 23-1 j (E—guBHO)G (E)-D E ajG (E) = 225/29, (5.20) 1—1 (: .._. FY) 21 where Z. = <[Al,S_]>. 1 g 9 Or, if we write this in matrix form we have, E-guBHO -D Gl(E) Z E ‘1 2 -gu8HO -D 1 G (E) Z E-guBHO . ‘*w/ G3(E) Z -alD -a2D -a3D . . . E-guBHO G (E) Z (5.21) This can even be shortened to, (E g - gCF)°g(E) = g/z. (5.22) where I is the identity matrix, the column vector G(E) has components Gi(E), and Z has components 21' As we noted in the spin 1 case the ensemble averages and the Zi are just linear combinations of the mo- ments <(S:)n>. Hence equation (5.21) completely determines these moments for all temperatures. This process of solving for the moments, <(Sz)n>, via the Green function technique is certainly not any easier than directly diagonalizing the Hamiltonian H, k>ut it is very revealing. Suppose we add an exchange term ‘to the Hamiltonian which we will denote as H for the time El ZS 22 being. We can again write down the set of equations for the ZS Green functions as in (5.19) and (5.20), but now there will be extra terms of the form <<[A;,HE];S;>> associated with each of these equations. These exchange terms are not completely determined within this set of equations, i.e. they are higher order Green functions. The only way the exchange terms.can be known within this set of ZS equations is if we reduce them to lower order Green functions by an approximation scheme. In this way we will be able to determine the moments when there are exchange interactions present. Note that in this process we will only need to approximate the Green functions asso- ciated with the exchange. Hence, we may state that: the decoupling of the crystal field Green function in (4.1) is completely unnecessary. The real problem in trying to extract statistical information from an interacting system of Spins with crystal field anisotrOpy present is in de- coupling all of the exchange term Green functions, which arise from using a set Of Spin Operators for the A Opera- . . . . + tor 1n (3.3) rather than the Single Operator Sg' B. The Exchange Term Decoupling In this section we will present the decoupling of the exchange term for an arbitrary size spin system. For con- venience we will temporarily neglect the crystal field anisotrOpy. The Heisenberg-Dirac Hamiltonian (1.3) can be written as, 23 _ z z 1 + - - + ‘ ’Z XJij[sisj + 2(Sisj + Sisj)]' (5°23) HE Since the RPA is one of the most convenient and reasonably accurate decoupling approximations used in determining the transition temperature we would like to be able to decouple each of the Green functions <<[A;,HE];S;>> so as to obtain the RPA as originally stated in (4.2). It will first be necessary to define two new sets of Operators {8;} and {C3} , as follows, 131 =[Al,S+] for all i, (5.24) 9 9 9 i i - . C =[A ,S ] for all 1. (5.25) 9 9 9 In the usual matrix representation, the spin angular mo- mentum operators A; are non—diagonal, more specifically their only non-zero elements are on the diagonal just above the main diagonal of the matrix. ‘Therefore by the defini- tion (5.24) the {B3} matrices have their only non-zero elements on the second diagonal above the main diagonal. Hence the ensemble averages of the {Ag} and {8:} Operators must vanish at all temperatures: = = O for all i. (5.26) Next we note that the Operators {C3} are diagonal, and also that the quantities Zi defined in the previous section are of the form Zi={C;}. Since the {C3} are diagonal they can be? written as linear combinations of the powers of the HEitrix 8; and the identity matrix I. One can also show 24 that for i Odd the matrices {C3} can be written solely as linear combinations of the odd powers of 8;, and if i is even, the {C;} can be written solely as linear combina- tions of the even powers of 8;. These relations will be useful later when we consider antiferromagnetic systems. Using (5.16), (5.24), and (5.25), the commutator of A; with the Hamiltonian HE in (5.23) yields, [Ai.H 1 = —)Jf g(- ~2Agi g E f 8+2 BiS ‘ CiS where we have made use of the facts that J = J , J fg gf ff and Operators located at lattice sites 9 and f automat- =0, ically commute unless g=f. Hence we have separated the h>> into three parts: <>, <> and <>.Tak1ng our initiative from the RPA as stated in (4.2), we try decoup- Green function <<[A;,HE ];S- 1ings of the form: 1 z - _ z i. - <> — <>, (5.28) i— — .- < > — 0 (5.29) _ _ i +._ _ + .. <> — <> — Z. l<>' (5.30) The first and third decouplings are rather obvious: we have factored out the powers of s; and replaced them by their ensemble averages. In effect, then, we are saying that the powers of s: are uncorrelated with the operators A%, where g and f are different lattice sites. Since the operators A; are strictly non—diagonal they represent 25 fluctuations which are transverse to the quantization axis 2. The powers of 8:, of course, correspond to longi- tudinal fluctuations. Hence it is not unreasonable to expect the longitudinal fluctuations at site g to be almost completely uncorrelated with the transverse fluctuations at site f#g, as long as we are in, or near, the ordered state where the predominant correlation between the spins is a longitudinal-longitudinal correlation. The decoupling of (5.29) is less obvious, but nevertheless reasonable. f and replace it by its ensemble average, we get zero because If we attempt to factor out either the operator 3; or S = = O as we stated earlier. Both of the opera- f therefore the Green function in (5.29) represents transverse- i - . tors Bg and S correspond to transverse correlations and transverse correlations which are much smaller than even the transverse—longitudinal correlations at the transi- tion temperature. Clearly, then, this Green function is negligible. The equations of motions are then, Zi/2n + 22J <>= (S \QH + _ (5.31) - gziJfg<>. The sum over all the lattice sites f complicates somewhat the process of obtaining a solution of (5.31). This can be overcome by assuming that the magnetic ions all lie on a Bravais lattice. (If this is not the case one must introduce the concept of sublattices. This procedure will 26 be discussed in the section on antiferromagnetism.) h .depend only on the difference, Secondly one must assume that the Green functions <> and the exchange integrals Jgh Bg-Bh' in their spatial coordinates rather than on the individual coordinates, R and Rh' Then we may Fourier ~ transform as follows, Gl(k,E) = X <> exp(ik°(R -R )), (5.32) ” R -R g ‘ ”g ”h ~g ~11 J k = J ex -ik° R -R 5.33 (~) R ER gh p< ~ (~g ~h>). ( ) ~g ~h - i i . = 2 <5 A > eXp(-ik-(R -R )), (5.34) ~ Rh-R 9 ~ “ ~9 ~ ~g The inverses of these transforms are, i _ l i _. . _ <> _(§)g G (5,3) eXp( 15 ‘89 gh)), (5.35) _.£) - . _ Jgh —;NIE J(g) eXp(lE (89 3h)), (5.36) 11 - i . ={§)) exp(15o(Ah-§g)). (5.37) ~ where N is the number of magnetic ions on the Bravais lattice. The sums over k are over all points in the first Brillouin zone of the lattice. Then defining §=, we have, E 61(5,E) = Zi/Zn + 2§J(0)Gi(5,E)-ziJ(5)Gl(§,E). (5.38) 27 This may be written in matrix form as, E-2§(J(0)-J(k)) 61 25 - 2 zlJ(§) E ZSJ(0) G 22 Z2J(k) E-2§J(0) ' G3 23 ~ , . . _1 . 3 . . 7;; . L - 28 zZSJ(§) E-ZSJ(Q) G . 223 ’ (5.39) where Gl=Gl(k,E), and also as, (E I - a )-G(E,E) = g/zn, (5.40) where 9(E,E) and E are the column vectors as defined for (5.22) except that now 9 depends on E as well as E, and gE is the matrix containing just the exchange integrals. We will now show that this system of equations repro- (4) duces the results of Tahir-Kheli and tar Haar which was the first successful attempt at using the Tyablikov de- coupling of (4.2) to solve the general spin S problem. To find Gi(§,E) from (5.39) one uses the well-known Cramer's Rule. G1(k,E) is the ratio of two determinants: (H ~ ~ Gi(E,E) = det(gi)/det(E (5.41) where the matrix M1 is formed from the matrix B I - HE by replacing the i th column with the column vector g /2n. Explicitly the matrix M1 is, ~ 28 gi= E-2§(J(0)-J(k)) zl/Zn z2J(5) E-28J(0) 1 \) 22/2“ //w\ . --/ . 23J(5) . . Q . \‘ ° . . O \\ O . 2253(5) ~f 225/2n ° E-ZSJ(O) (5.42) The determinant may easily be obtained from (5.42) by expanding in minors about the first column. All the minors of the elements in the first column are zero except for the first and i th elements, as one can readily see by inspection. Then, det(gi) = (E—2§(J(0)—J(5))) (Zi/Zn) (E-2§J(0))25‘2 1+1 (5.43) + (-1) ziJ(E)X(Minor of the element ZiJ(E)). The minor of the element ZiJ(h) turns out to be (-l)l(Zl/2W)(B-2§J(O))2$_2. Hence, det(¥i) = (E-2§(J(0)-J(g)>)(zi/2n)(E-2§J(0))25‘2 (5.44) - ZiJ(E)(Zl/2n)(E-2§J(O))ZS-2. Since Zl=2§ this reduces to, det(gi) = (zi/zn)(E-2§J(0))25‘l. (5.45) 29 From (5.39) we can also quickly obtain, 1 det(E E - H = (E-2§(J(0)-J(g)))(E-2§J(0))25‘.'(5.46) ~E) Therefore, Zi . (5.47) E-2§(J(0)—J(§)> i l G (151E) - E' The wavevector dependent correlation function (5.34) may be gotten from (5.47) by means of the corollary to the Spectral Theorem (3.6), Z. - i _ i <8 (§)A > — exp(E(E)/KT)-l , (5.48) where E(§)=2§(J(O)—J(§)). Then using (5.37) the Spatially- dependent correlation function becomes, ° 1): Z' (5.49) -l__ l _'. _ (S A9> "N k exp(E(k)/KT)—1 eXp( 15 (59 Bh))° For g=h this reduces to, = Zi\%) E l/(exp(E(§)/KT)—l), (5.50) If we define, ~ 4 =(%)g l/(exp(E(E)/KT)-l)’ (5.51) then (5.50) becomes, = <[Ai,s‘]>¢ for i = 1,2,...,2s. (5.52) / _ . . _ As we have mentioned before, both f8 Al> and <[A1,S ]> . . ' l are linear combinations of the moments <(Sz)n> for n = 0,1,2,...,25—l. Hence (5.52) is a set of 2S equations 30 in the 28 unknown moments. Therefore (5.52) completely determines the moments. This set of ZS equations is com- pletely equivalent to the set of ZS equations contained in equation (3.13) in the paper by Tahir-Kheli and ter Haar. All that we have done here is to recast their solution into one that uses the Operators {A:}” Therefore the decoupling that we have used in (5.28)-(5.30) is completely equivalent to the RPA as stated in (4.2). C. The Complete Solution The matrix equation (5.21) is applicable for a sys- tem with no exchange, but just an external field and crys- tal field anisotrOpy. The matrix equation (5.39) is for a system which only has exchange coupling between the spins. The summarized form of each of these equations is, Crystal field: (El-8CF)'G = g/Zfl, (5.53) Exchange: (El-hE)'G = g/Zn. (5.54) Then if the system is described by a Hamiltonian H=H +H CF E’ the equations of motion become, (EI~H z :CF—g E)-§ = g/ZW, (5.55) or, if we define g = g + H (EE-§)-§ = g/ZN. (5.56) 31 Writing this out explicitly, we have, E-2§(J(0)-J(k)) -D Z2J(k) E-2§J(0) -D -Z3J(k) E-2§J(0) . . . -D ZZSJ(k)-alD -a2D -a3D . . . —aZS-1D E-28J(0) 1 G (5,3) 21 62(5,E) 22 (5.57) 6} w .5? ['11 II N'H :1 coo-N w G28(k,E) z This set of ZS equations completely determines the moments <(Sz)n> of the system that includes in its Hamiltonian exchange interactions and crystal field anisotrOpy. The only approximation that we have had to make is in decoup- ling the exchange terms via the RPA. The spectrum of excitations, i.e. the location of the poles of the Green functions is obtained by equating the determinant of the matrix EI-H to zero, ~ det(ET-H) = 0. (5.58) Since the determinant is a polynomial of degree 25 in the 32 variable B, there are ZS possible roots of this equation. Denoting the roots as Em(k), we can write, ZS aet(E§-H> = 1?] (E-Em(k)). (5.59) ~ 3 m=l ~ As before, the Green functions Gi(k,E) are obtained by divid- ing (5-59) into the determinant of the matrix which is formed by replacing the i th column of the matrix ET—H with the column vector Z/Zn. The determinant of this newly formed matrix is a polynomial of degree 28-1 in the variable E. Hence a partial fraction expansion will yield a Gi(k,E) which is a sum of simple poles located at the points Em(k) on the real axis in the complex E plane. Notationally we may write, i 1 2% Rg‘k’ (5 60) G (k,E) = —-' 'fi:-—:——'p - ~ 2n m=l E Em(k) where R;(k) is the residue of the m th pole of the i th Green function. The wavevector dependent correlation func- tion is then, i 25 Rmflg) .. i .- - mil exp(Em(k)/KT)-l ' (5.61) Then the spatial correlation function is, . zs 125L (k) =(l) X m ~ (5.62) N =1 l 2 . _ exp(eik°(R - ). g k m exp(Em(k)/KT) 1 ~ ~g §h and in particular, i - ZS R (k) - i l m ~ =(—)2 E _ . (5.63) h N k m=l exp(Em(k)/KT) 1 This set of ZS equations contains all the necessary 33 thermodynamic information. Unfortunately, (5.63) cannot be written in any simpler way than it is for a given size spin system. The complete solutions can only be obtained by a high Speed computer. D. The Transition Temperature The transition temperature is that temperature above which the system is no longer ordered when there is no external magnetic field present. A system which is not ordered has all of its odd moments ((Sz)2n+l> = 0, since the unordered state is a time reversible state. We may use this fact to obtain Tc’ the transition temperature. We cannot, though, merely set the odd moments in (5.63) to zero and expect to obtain a complete set of equations that are useful for solving for TC. For the particular cases worked out by this author, this process yielded redundant and sometimes trivial equations. The correct procedure for obtaining a complete set of equations, one unknown of which is Tc, is as follows. First define the quantities, _ Lim 82n+1 ‘ s+o (Z 2n+l/S)' (5.64) Then insert for Z in (5.63) the quantities 82n+1§' 2n+l and then expand the functions on the right hand sides of the equations in the variable 5. Comparing the coeffi- cients of the zeroth and first powers of S in each equa- tion yields 28 independent equations. The ZS unknown quantities will then be the quantities Zi for i even, the 34 quantities Bi for i odd, and TC(B1 is already known to be just the number 2). These equations must also be solved on a high speed computer. Unfortunately, again, no general expression for Tc has been found for an arbitrary size spin system. Chapter 7 of this thesis will present results for some special cases that have been worked out. E. A Decoupling Approximation for the Molecular Field Theory Suppose we had chosen to make the following decoupling approximation, <> = <52><>dg,h' (5,55) LOP HAN <> ;s'>>6 (5.66) + ,<< Zi Sh HT+ g,h' LQP instead of the ones we had made in (5.28) and (5.30). The equations of motion would then have been, E-2§J(0) G1(E) zl E-2§J(0) G2(E) 22 E-2§J(0) G3(E) z3 1 O 0 =5? 0 -a D -a D -a D E-2§J(0) G25(D) z 1 2 ° ’ ' ‘ 25-1 25 ' (5.67) where G1(E) = <>. Note that we could have gotten 35 (5-67) from (5.57) just by setting J(k) = 0. One can also see that (5.67) could have been obtained from (5-56) if we chose the Hamiltonian to be the MFT Hamiltonian, HMFT = E Hi’ (5.68) where Hi is as defined in (2.4). Therefore we are saying that the decoupling of (5.65) and (5.66) is completely equivalent to the MFT approximation, which, in turn, is equivalent to setting J(k) = O in (5-57)- The fact that setting J(k)=0 in the RPA equations of motion yields the MFT results when there is no crystal— field present (D=O) has been known for some time(15). Only until now is it known also to hold when D is finite. None of the elementary theories presented in Chapter 4 gives the MFT result when J(k) is set equal to zero. The decoupling of (5-65) and (5-66) is certainly consistent with what we already know about the MFT. This is a theory which neglects all correlations between spins at different lattice sites, as indicated by (2.2) where we have replaced neighboring spin operators by their en— semble averages, instead of considering their full quantum mechanical effect. 6. THE MANY-EXCITATION THEORY FOR AN ANTIFERROMAGNETIC SYSTEM All of our previous discussion has been concerned solely with ferromagnetic systems, i.e. systems in which the aligned spins all point in the same direction. More frequently though in insulating magnetic salts the Spins are arranged in parallel-antiparallel pairs so that the net magnetic moment of the crystal is zero. Such a con- figuration is called antiferromagnetic ordering. The simp- lest kind of such a system is when the nearest neighbors of each spin are antiparallel to it. Clearly such a lat- tice is not translationally invariant, but it is possible to break the lattice down into translationally invariant sublattices. We may define a magnetic sublattice as a Bravais lattice, which contains only spins that point in the same direction. We can always perform the subdivision of a lattice into many sublattices, though for complicated spin structures this may not be an easy thing to do. For convenience we will only consider a two sublattice model in this thesis. The extension to a larger number of sub— lattices is relatively easy to deduce from the two sub- lattice case. (7) (8) first solved the Anderson and Callen and Lines 36 37 antiferromagnetic problem within the Green function form- alism. We can use their technique for our problem also. Following their technique, we will consider two distinct types of Green functions. The first kind of Green func- tions will be concerned with correlations between ions that lie on the same sublattice, the second kind of Green functions will be concerned only with correlations between ions that lie on different sublattices. The second kind of Green functions are not really Green functions at all, since their equations of motion do not have the inhomo- geneous term. We will label the sublattices 1 and 2 respectively. Sublattice 1 will be the sublattice containing the "up" spins, and sublattice 2 will contain the "down" spins. Lattice site index g1 will refer to a point on sublattice 1, and the index g2 will refer to a point on sublat- tice 2. Rewriting equation (5.31) for an antiferromagnetic system yields, E<> = a Z./2w + 22 <5 2 i > g1 hl gl,hl 1 fl fl Jf1,g1<> z i - + ZEZJf2,gl<> i + — Jf1,g1<> . + _ gl>Jf2,gl<>’ (6'1) 38 Z > <> ‘ 22 > z >J <> f2 i Jfl,gZ<> 6.2 f2 g2 ( ) + _ Jf1,g2<>° We must be more careful now in how we use the ensemble averages and . Since the direction of the spin alignment is different on each sublattice, these quantities are no longer strictly positive nor are they independent of the index f, as they were before. We shall eliminate some of the confusion by adopting the following convention: = +8, (6.3) <5z > = -§ (6 4) f2 I o since 1 is the "up" sublattice and 2 is the "down" sub- lattice. The ensemble averages of the odd powers of S2 must also possess the same prOperty of changing sign from sublattice to sublattice. Recalling that Ci depends only on the odd powers of S2 if i is odd and only on the even powers of S2 if i is even, we adopt the following: (Cfl> = zi’ (6.5) 1 _ _ i — ( l) Zi' (6.6) and therefore the quantities S and Zi will then be strictly positive. 39 By again assuming that the Green functions and the exchange integrals depend only on the difference in their two spatial coordinates, we may perform Fourier trans- formations as follows: i i _ Gl(k,E) _ R ER <>exp(ik (R 91- Rhl)), (6.7) ~gl ~h1 i _ i . - - . _ G2(k,E) — R E <>exp(ik (R92 Rhl)), (6.8) ~92 Bhl Jll(k) = R R ng’hlexp(ik° (R gl -Rhl)), (6.9) ~gl ~hl (121(5) = R ER ng,hlexp““f°(8g2’§hl”' (6.10) ~g2 ~hl J 2(15) = R E ng'hZeXpHIE (R g-l Rh2)), (6.11) ~gl BhZ ~q2 ~h2 - i _ - i . _ _ R )— exp(ik-(Rgl Rhl). (6.13) gl Rhl The inverses of these transformations are as follows: <> _(Nl£ G2(k,E) exp( 15 (392 Bhl))' ~ 40 ng’hl =ifi‘é Jll(k) eXp(1R-(Rgl-Rhl)), (6.16) JgZ,hl = $1; J21(k) exp(ik°(§g2-§hl))p (6.17) ng’hz = é E J12(k) exp(ik'(§gl-Bh2)): (6.18) J92,h2 = g g J22= § 3 exp(-ik°(Bgl-Bhl))- (6.20) The sums over k are to be performed over the first Brillouin zone of the sublattice. The factor Z/N corresponds to the fact that there are N/Z sites on each sublattice, and hence N/Z k vectors in the first zone. Using the transformed quantities in (6.7)-(6.1Z) the equations (6.1) and (6.2) can be written concisely as, i _ - i _ - i EGl(k,E) — Zi/Zn + ZSJ11(O)Gl(k,E) ZSJ21(0)Gl(k,E) l l 1 _ — 1 _ - 1 - (-1)iz J (k)Gl(k E) - (-1)iz J (O)Gl(k E) 1 12 2 1 ~' 1 22 2 ~' ° Since the crystal field anisotrOpy acts on each spin at each lattice site separately, rather than coupling the spins to each other, it is easy to reintroduce the aniso- trOpy at this point, 41 1 _ - 1 _ - 1 EGl(k,E) _ Zi/Zn + 2SJll(O)Gl(k,E) ZSJ21(O)Gl(k,E) 1 1 ‘21311‘E’G1‘E'E’ - ziJ21(§)62(5.E) -DGi+l(k,E), for 1 s 2s-1, (6.23) 1 _ - 1 _ — 1 EG2(k,E) _ 28J12(0)G2(k,E) 28J22(0)Gz(k,E) i l i l - (-1) ziJ12(g)Gl(5.E) - (—1) ZiJ22(E)Gz‘E'E’ - ch+l(R,E), for 1 s 2s-1, (6.24) 1 _ — 1 _ - 1 EGl(k,E) _ Zi/Zn + 2SJll(O)Gl(k,E) ZSJ21(O)Gl(k,E) l l ‘21J11(E’G1(E'E’ - ziJ21(5)62(5.E) ZS-l . - ) Da.Gi(k,E), (6.25) j=1 3 ~ 1 _ - 1 _ - 1 EGZ(k,E) _ ZSJ12(0)G2(k,E) 2SJ22(0)G2(k,E) i l i l - (-1) ziJ12(5)61(r,E) - (-1) ZiJ22(k)G2(EIE) 25- ‘ ,l i=1 1 . J DajG2(k,E). (6.26) This set of equations may be written in block matrix form as, E3"§11 _§21 91 g N. NH :1 EI-H G (6.27) “E12 z 222 ~2 (O ‘ 42 where the elements are matrices themselves, which are defined as follows: 31-511 = E-28(J(0)-Jll(k)) -D 22 Jll(k) E-ZSJ(0). . Z3J11(§) ° . o . . —D Z28 Jll(k) -alD -a2D . . . °_aZS-1D E-ZSJ(0) Eg-gzz = E+2S(J(O)-Jll(k)) -D zlel (15) E+ZSJ(0). . . . z3 J11‘k> / . ' . ; \ ‘ -D (-1)252 J (k)-a D -a D —a D E+2§J(0) 2S 11 ~ 1 2 ° ° ' ° 23-1 '221 = 2 J 21(k) 21J21(k) 22J21(E) ‘ . \ ZZSJ(E) I (6.30) 22J12‘5) -z3J12(k) (-1)lz 25J12(¥’ 43 6% (15$) 2 G2 (15.12:) 3 62 (15.12) 28 62 (15.13) 0 o o o O (6.31) , (6.32) (6.33) The quantities Jll(k) and J22(k) are equal for a system such as ours which has two completely equivalent sublat- tices. By the same reasoning, J12(k) and J21(k) must be complex conjugates of each other, and for k=0 they are both real, which implies J12(0)=J21(0). The symbol J(O) 44 here represents either the quantity Jll(0)-J21(O) or J22(0)-J12(0). Note that the matrix Eg-gll matrix of (5.57) with J(k) replaced by Jll(k). The matrix ET-HZZ is the same as EI-Hll except that the signs before the quantities Zi are changed. is just the Equation (6.27) is really a set of 48 coupled equations in the 48 Green functions. There are, however, only ZS unknown thermodynamic quantities, the Zi’ or equivalently the moments of 82. All 48 Green functions depend on these ZS quantities, but it is only the Green functions of the form Gi(k,E) which can be related back to these moments via the Spectral Theorem, and hence they provide a set of ZS equations which completely determine the moments. The Green functions G§(k,E) certainly depend on these moments, but they cannot be used to find these quantities. The Green functions Gi(k,E) can be found by the same technique that we have used previously. Each of the Green functions is the ratio of two determinants, Gi(k,E) = det<§i)/aet(§i), (6.34) where, 5i = EE'EII ‘§21 -§12 51-522 (6.35) The matrix K1 is just the coefficient matrix of (6.27), and ~ the matrix N1 is formed by replacing the i th column of El ~ ~ 45 with the column vector, 1. 2n . (6.36) (GIN One can show that the determinant det(§l) is a polynomial 2 ~ in the variable E rather than the variable E. Hence the Green functions Gl(k,E) will have poles that come in posi- tive and negative pairs. Clearly then, Gi(k,E) has the form, i i Gi(k E) = l. gs .EE£§1_ + 13355:.., (5 37) 1 -' 21 _ E-E (k) 3+5 (k) ' ' m- l m ~ In ~ _. where the Em(k) are the ZS possible positive poles, the R;(k) are the residues of the positive poles, and the T;(k) are the residues of the negative poles. Then we have, . 25 r Ri(k) -1 _’-2_ ‘ m~ ‘ (R) g mil: exp(Em(k)/KT)-l T;(5) + exp(-Em(k)/KT>-l}" (6°38) It is the ZS equations of (6.38) that can then be used to determine the ZS moments <(Sz)n>. We can find the transition temperature for the anti- ferromagnet in the same way that we found it for the ferro- magnet. This method was described in section D of Chapter 5. The transition (or critical) temperature for an antifer- romagnet is called the Néel temperature, and is frequently denoted as TN by other authors, but we will simply con- tinue to use Tc here. The technique that we have presented here for an antiferromagnetic system which is divisible into two 46 sublattices is also useful in working with a ferromagnetic system which has two sublattices. The only change we need to make is in the matrix EI-HZZ. Recall that for an anti- ferromagnet we had to be careful about the signs of the quantities . For the ferromagnetic case these are all strictly positive. Then for a two sublattice ferro- magnet the equations for the Green functions are as in (6.27) only now the matrix EI-H22 is the same as the matrix EI-H and H i.e. there are no sign changes in z :11 212::21' the coefficients of the quantities Zi. 7. THE RESULTS FOR SOME SPECIAL CASES A. The Spin S=l Ferromagnet For S=l we need only the operators A1 and A2, which are S+ and SZS++S+Sz respectively. Then (5.56) becomes, —E-2S(J(0)-J(k)) -D TG1(R,E)W (2S — 2 Z2J(k)-D E-ZSJ(0)G _6 (3,5)- _224 (7.1) The poles of the Green functions are gotten in the usual way by equating the determinant of the coefficient matrix in (7.1) to zero, and solving for the roots of this equation. The two poles are, 2) 61(5) S(2J(0)-J(k)) + (DZ-22DJ(R)+§2J(§> .1 2 1 2 2 DJ(§)+§23(R) E2(k) S(2J(0)-J(k)) - (DZ-z2 ) Equations (5.63) become, Ri(E) ‘Lexp(El(k)/KT)-l (S-S+> = (1%) 1 R2(5) + exp(52(§)/KT)—1f ' - z + + z -1 RT(E) (S (S S +8 S )> = (E); .eXp(El(E)7kT)-l R§(k) 5 exp(E2(k)/ T)-l 3' 47 (7.2) (7.3) (7.4) (7.5) 48 where the residues are, Ri(k) = (2681(§)-4§2J(0)+z2D)/(El(5)-EZ(5)), (7.6) R;(k) = —(2§E2(§)—4§2J(0)+22D)/(E1(E)-EZ(E))’ (7.7) Ri‘?’ = (2281(5)-2522J(0)+2§D)/(El(§)-Ez(5)). (7.8) R§(k) = -(Z2E2(k)-2S22J(O)+2SD)/(El(k)-E2(k)). (7.9) By writing out the angular momentum matrices explicitly one finds that, Z = 6<(Sz)2> -4, (s s > = (8-S-22)/6, = 5-22/2. Then (7.4) and (7.5) become, _ Ri(E) (8-6S-ZZ)/6 = /KT)'1 1 R2(k) + exp(E2 (JE)/KT)"l/k' _ ( Ri(k) 8’22/2 = (exp(El(k)/KT)-l R§(k) + exp (E2 (k)/KT) -l>k' where we have used the notation, (7.10) (7.11) (7.12) (7.13) (7.14) 49 \/ >15 =§E (7.15) For any given temperature (7.13) and (7.14)determine S and 22 completely in terms of the exchange integrals Jfg and the anisotropy D. Using the technique outlined in Section D of Chapter 5, the equations that determine Tc and Z2 evaluated at T=TC are as follows, _ D (8 - Z2)/6 — (22/2) (370$) coth(EO(k)/2KT)>k, (7.16) ~ 4(Z2-2)/3Z 2 _<(J(1§)-2J(0)) exp(Eo(1§)/KTC)> 2 _ KTC(exp(EO(k)/KTc)-l)2 ~ 1 J(k) _ 2'E’ (7.17) 1 where 80(5) = (Dz-ZZDJ(k))2. Figures 1 and 2 contain plots of TC(D)/TC(O) as a function of D/J a simple cubic (SC) and a body centered cubic (BCC) lattice respectively. We have assumed that the exchange J is between nearest neighbors only, hence the variable J(O) = zJ, where z is the number of nearest neighbors. The Curie temperature TC is finite for all possible values of the ratio D/J as it is supposed to be. At very low anisotropy our results agree very closely with those of Lines. This is not so surprising since his tech- nique is certainly the most sophisticated of all the single excitation theories, and should be reasonably accurate for small values of D/J. For extremely large 50 Figure l. The ratio TC(D)/TC(0) versus D/J for a simple cubic spin S=l ferromagnet. (a) Narath's results. (b) Anderson and Callen's results. (0) Lines's results. (d) The results of this thesis. (e) The molecular field theory results. I6 51 A U v a? A v 8 A II .o (D V 0 <0 .. q E l l I 4 ¢ N l.0 (o)°1/(a)°1 4.0 3.0 2.0 D/J n.0, 52 Figure 2. The ratio TC(D)/TC(0) versus D/J for a body centered cubic spin S=l ferromagnet. (a) Narath's re- sults. (b) Anderson and Callen's results. (c) Lines's results. (d) The results of this thesis. (e) The molecular field theory results. 53 |.6 Q E A 32 .0 .1 T U) -( U 0 m - ‘8 1 q n 1 1 1 ‘2 c\! 0. (0)°1/(0)°1 4.0 3.0 2.0 D/J l.0 III—r v— 1.2,. - 54 anisotrOpy, one may expand functions in (7.16) and (7.17) to obtain KTC=2 J(O) at D/J = m. This is the same result that the MFT gives in the limit D/J+m. Hence for large values of the anisotropy our RPA results approach the results of the MFT. For each value of D/J both of our simple cubic and body centered cubic calculations gave the same value to Z at T=Tc as did the MFT. Therefore the value of Z <(Sz)2> at T=Tc for a given value of D/J was the same in each of the three calculations. Figure 3 contains plots of the sublattice magnetiza- tion, S, for a simple cubic lattice as a function of T/Tc for three different values of D/J. Note that the larger D/J the "squarer" the curve becomes.. This is to be expect- ed, since the larger D/J the higher one must go in tempera- ture, i.e. closer to Tc' before the intermediate spins states, |S-l>, |S-Z>,...,|-S+l>, become populated enough to reduce the ensemble average significantly. Also plot- ted in Figure 3 is the MFT prediction of S as a function of T/Tc. For all values of T/Tc the MFT curve lies above our curves. This MFT curve could just as well represent D/J = 1.0, 0.1, or 0.01 since the MFT does not distinguish between these values within the linewidth we have chosen for our graph. Hence the Green function calculation pre- dicts greater changes in S versus T/TC for different D/J ratios than does the MFT. This is consistent with the (14) known fact that the MFT underestimates the effect of the crystal field anisotropy on all the statistical 55 Figure 3. The sublattice magnetization versus T/Tc for a simple cubic spin S=l ferromagnet. Curves (a), (b), and (c) are the Green function calculations of this thesis. Curve (d) is the molecular field theory for the same three values of D/J, but in this case the three lines are indistinguishable. 0.8 ” 0.4 — 0.2 ” 56 (C) (0) (b) 0 0.2 04 05 T/TC 0.8 57 mechanical quantities. In figure 4 we have plotted <(SZ)2> versus T/Tc:for a simple cubic lattice with D/J = 1.0. For T 5 TC the Green function curve lies beneath the MFT curve, but for T 2 Tc’ the reverse is true. This behavior is typical of all values of D/J. The ensemble average <(Sz)2> is greater than the isotrOpic value S(S+l)/3 for all values of D/J. Only for D/J = 0 will <(sz)2> = 5(s+1)/3. B. The Spin S=l Antiferromagnet For an antiferromagnetic system with S=l we consider the four Green functions Gi(k,E), Gi(k,E), G§(k,E), and G§(k,E). Using the formalism presented in Chapter 4, the equations of motion are, rE-ZS(J(O)-Jll(k)) -D 2§J21(k) 0 Z2Jll(k)-D E-2§J(0) Z2J21(k) 0 -2§J12(k) o E+ZS(J(0)-Jll(k)) -D LZ2J12(k) o 22J11(R)-D E+ZSJ(0)J FGl(k,E)q .55. 1 ~ Gi‘E'E) Z2 1 JX< G:(k,E) =5? 0 (7.18) .Gi‘E'E’. .01 . 58 Figure 4. The ensemble average <(Sz)2> versus T/TC for a simple cubic spin S=l ferromagnet. (a) The Green function calculation of this thesis for D/J=l.0. (b) The molecular field theory prediction for D/J=l.0. Note that at T=TC the curves cross each other. 59 0.N m.- oexe of 0.0 0.0 [ A3 at E 1 50 - 0.0 1 0.0 0.. Ava 60 Clearly, there are four simple poles to be associated with each of the Green functions, but as we noted in the previous chapter, they come in positive and negative pairs. The two positive poles are, 131(k) = +(a + a + (a3)2);5.);5, (7.19) l + (4ala 2 2 E2(k) = +(cxl + a2 — (4ala2 + (d3)2)%)%, (7.20 where, _ 2_ ‘2 _ a1 — D ZZDJ11(0)+S (Jll(§) J12(§))(Jll(g)+J12(5)), (7.21) _ —2 _ _ _ a2 - s (Jll‘E’ J12<5) 2J(0))(Jll(g)+J12(g) 2J(0)>. (7.22) _ _ —2 a _ J12(§)(zzn 4s J(O)). (7.23) In the above calculation we have assumed that we are work- ing with a lattice which has the prOperty that J12(E)=J21(g), i.e. the quantity J (E) is assumed to be real. This is 12 not necessarily so for all lattices. The form of Gi(E,E) is, , Rl(k) Tl(k) &mE)=L5 1~ + 1~ 1 ~' 2N gE-ElZES E+ 1Zk$ \ ~ ~ Rl(k) Tl(k) 7 -+._:E_:__ + 2 ~ .. (7 24) E-E2(E) E¥E2(KY ;' ° 61 where we may write the residues in the form, 1 R1(E) 1 R2(§) 1 T1(E) 1 T2(E) where, 82(5) 33(5) 34(5) 3 2 (E1(E>31+E1(E>32(E) + El(g>B3<§) 2 2 +B4<§)>/<2El(5)(El<5>-EZ<§>>>, _ _ 3 2 — (E2(E)Bl+E2(E)BZ(E) + E2(E)B3(E) 2 2 +B4(E))/(2E2(E)(E1(E)-E2(E)))p 3 2 (El‘E)Bl‘E1‘E’Bz(E’ + El‘E)Ba‘E’ 2 2 +B4<§))/<2El(5>(51(5>—Ez(§))). 3 2 -(EZ(E)Bl-E2(E)BZ(E) + 32(E)B3(E) 2 2 +B4(E))/(2E2(E)(31(E)-E2(E)))p 28, —2 z D + 45 (J(0)-Jll(§)). 2 3 2 4§ZZDJ(O) - 16§ J2(0) - 2§D , 2 —2 _ D Jll(E) + 45 ZZDJ(O)(J(O) 2Jll(§)) 2 Z2 - 4§2J(0)D2 - 2203 -4 2 -lGS J (0)(J(O)-J11(E)). (7. (7. (7. (7. (7. (7. (7. (7. 25) 26) 27) 28) 29) 30) 31) 32) 62 Similarly, one has for Gi(k,E) the expression, {‘ R2 (k) T2 (k) G2(k E) = i— 1 ~ + 1 ~ 1 ~’ 2n )E-E1(E) E+El(~) R§(k) T§(k) + ————l—— + ————1—— , E-E2(E) E+E2(§)J where the residues are of the form, R2(k> = (E3(k)c +E2(k)C (k)+E (k>c (k) 1 ~ 1 ~ 1 1 ~ 2 ~ 1 ~ 3 ~ 2 2 + c4cg)>/(2El<5>(31(5>-E2(g)>>, 2 _ _ 3 2 R2 - <32c2(g)+32<5>c3(g) + c >/<2E >> 4 ~ 2 ~ 1 ~ 2 ~ ' T2 = (E3(k)C +E2(k)C (k)+E c (k) 1 ~ 1 ~ 1 1 ~ 2 ~ 1 ~ 3 ~ + c )/(2E (k)(E2(k)-E2(k))) 4 ~ 1 ~ 1 ~ 2 ~ ' T2(k) = -(E3(k)C +E2(k)C (k)+E (k)C (k) 2 ~ 2 ~ 1 2 ~ 2 ~ 2 ~ 3 ~ + c (k>>/(2E (k)(E2(k)_E2(k))) 4 ~ 2 ~ 1 ~ 2 ~ ’ where, C1 = 22, c2(5> = 23D - 25(J(0>-J11. c (k) = -z D2+ZZDJ (k)-4§2(J(0)-J (k))(Z (J(O) 3 ~ 2 2 11 ~ 11 ~ 2 + Jll(]i))—D)’ (7. (7. (7. (7. (7. (7. (7. (7. 33) 34) 35) 36) 37) 38) 39) 40) 63 2 2 3+4§z D2 c4(§) = -232 2 2 _ _ DJ11(E) - 28D Jll(E)-ZSD(J(O) -3 ‘J11)+8S J(o>(J(o>-Jll<§))(22J(o>-D). (7.41) Then using the corollary to the Spectral Theorem, (7.24), (7.33), and (7.10)-(7.12), we have, l l R (k) T (k) _‘_ _ 1 ~ 1 - (8 S ZZ)/6 - <;Xp(El(E)/KT)-l + exp(-E1(E)/KT)-l R§(5> R§(5) + exp(E2(§)7ET)-l + eXP(-E2(§)/KT)-;Zk’ (7.42) _ Ri T§ S - Z /2 = ” + ” 2 <;xp(El(§)/KT)-l exp(-El(§)/KT)-l R§<§> 23(5) 'exp(E2(E)/KT)-l + exp(-E2(E)/KT)-;>k° (7.43) The above two equations are ones that determine S and Z2 for all temperatures as a function of D and the Jfg Using the technique outlined in Section D of Chapter 5, we obtain the following two equations which determine the Neel temperature, Tc’ and Z evaluated at T=TC, 2 ./22D Z2D x + W COthEI (7.45) where, 2 F (k) = (2J<0)-Jll(g))(n 'ZzDJ11(E)’+ZzDJ12‘E’(J(O"J11‘E)), 1 . zZDJ12(E)EO(5) (7.46) 2 F (k) = —(2J(0)-Jll(g))(D -22DJ111E))-ZZDJ12(§)(J(0)-Jll(g))' 2 “ z29J12(}E’Eoo(}f) (7.47) E (k) = (DZ-Z D(J (k)-J (k)));5 (7 48) 0 ~ 2 11 ~ 12 ~ ' ° E (k) = (DZ-Z D(J (k)+J (k)))* (7 49) 00 ~ 2 11 ~ 12 ~ ° - Solutions of (7.44) and (7.45) have been obtained for simple cubic and body centered cubic lattices assuming only nearest neighbor exchange. The numerical results obtained for the Neel temperature T0 are equal to the results obtained for the Curie temperature of a ferromag- netic system under the same conditions. a priori there is no reason to expect that the Curie temperature should be equal to the Néel temperature for equivalent systems. For systems which have no crystal field anisotropy, Rushbrook and Wood(l6) obtained a larger transition temp- erature for the antiferromagnet than for the ferromagnet. The RPA, though, does predict rigorously the equivalence 65 of the transition temperatures in the case of D=0. For our calculation it is not easy to see at first that the temperatures are equal; it is only after the lengthy com- puter calculations that we find this out. Evidently, the decoupling of (5.28)-(5.30) neglects the fine distinction between antiferromagnetic ordering.and ferromagnetic ordering in the region of the phase transition. The ensemble averages <(Sz)2> also agree with the predictions obtained in the ferromagnetic case for each value of D/J. Figure 5 shows S versus T/TC for the same three values of D/J as we used in Figure 3, which is the same plot but only for a ferromagnetic system. This curve is quite different than the one for the ferromagnet, especially at low temperatures. The sublattice magnetization, S, does not reach its full value of unity for any of the D/J ratios, whereas in the ferromagnetic case S=l at T=O for all values of D/J. Our results are consistent with (17) in that they the antiferromagnetic spin wave results both have S less than unity, even at T=O. For larger values of D/J the value of S at T=0 is larger, and only for D/J = w will S=l. . The average <(Sz)2> is plotted in Figure 6, versus T/TC, and it possesses the same characteristics as the sublattice magnetization curve. For T=O, the average <(Sz)2 > does not approach its ferromagnetic value of unity, but some other value, which is slightly less than unity. The larger D/J is, the larger the T=O value of 66 Figure 5. The sublattice magnetization versus T/Tc for a simple cubic spin S=l antiferromagnet, (a) D/J=0.0l. (b) D/J=O.l. (C) D/J=l.0. 68 Figure 6. The ensemble average <(Sz)2> versus T/Tc for a simple cubic spin S=l antiferromagnet. (a) D/J = 0.01. (b) D/J = 0.1. (C) D/J = 1.0. 69 0.0 0.0 o»\» v.0 «.0 0 0.0 NO 0.0 0.0 0.. ~ A umv 70 <(Sz)2>. The results for T 2 TC are the same for the average <(Sz)2> in the antiferromagnetic case as for the ferromagnetic case which we have plotted in Figure 4. C. The Spin S=3/2 Ferromagnet For a spin S=3/2 the equations (5.56) become, E—2§(J(0)—J(g)) -D o 31(g,E) 2§ Z2J(E) E-2§J(0) -D 62(g,E) =%F 22 Z3J(E) -4D E-2§J(0) G3(g,E) z3 . (7.50) The poles of these Green functions are located at, 31(3) = A + 2(13/3);5 cos(Q/3), (7.51) 22(5) = A + 2(13/3);5 cos(Q/3 + 2n/3), (7.52) 53(3) = A + Hrs/3);5 cos(Q/3 + 4w/3), (7.53) where, A = 2§(J(0)—J(g)/3), (7.54) B = 4D2-Z2DJ(E) + 4§2J2(g)/3, (7.55) c = 222§DJ2(§)/3 + 16§DZJ(g)/3 - Z2D2J(E) - 16S3J3(K)/27, (7.56) cos 9 = (3C/2B)(3/B)% (7.57) The residue R;(h) of the m th pole of the i th Green 71 function is of the form, i _ _ i 2 i i R1 — (33(5) 82(8))(c181(g)+c281(g)+c3)/c4, (7.58) i _ _ 1 2 i i 82(5) — (81(5) 83(5))(0182(5)+c282(§)+c3)/c4. (7.59) i _ _ i 2 i i ' 83(5) — (82(5) 81(8))(clE3(§)+c283(5)+c3>/c4. (7.60) where, c4 = 28(8)15 sin 0, (7.61) l _. cl = 28, (7.62) 1 -2 c2 = -BS J(O) + Z2D, (7.63) c; = 2§(4§2J2(0)-4D2)-2§ZZDJ(O)+23D2, (7.64) 2- cl - 22, (7.65) c2 = -4SZ J(O)+Z D (7 66) 2 2 3’ ° c2 = 4§2z J2(O)-2§Z DJ(O), (7.67) 3 2 3 3 = 68 c1 Z3, (7.) 3_ -_ c2 _ 4220 4sz3J(0), (7.69) cg = ~8S22DJ(0)+4§2Z3J2(0). (7.70) Using the technique outlined in section D of Chapter 5 the equations which determine Tc are, 72 Tc = (4-222/3)/]5 , (7.71) _ D 6 — ZZZ}§+KTCE -( 22N()§)exp(Eo(l§)/KT) > KTCE§(§)(exp(EO(5)/KT)-1)2) E 2 “<)Tc/Eo‘5”‘zzL(E’ - 422J(0) + 81D 2 + (4Z2J (O)-281DJ(0))/L(E)£>E7 (7.73) where, 8(5) = 4(8-Bl)D2/(2(J(0)-J(§))(4D2-ZZDJ(§)) 22 2 2 - 2DJ (E) + 81D J(k)), (7.74) _ 3 _ 2 8(5) — (4z20 J(E)(8 Bl)/EO(E))(2(J(O) 2 2 2 — J(5))(4D -Z2DJ(E))-222DJ (5)+ 10 J(5)). (7.75) _ _ _ 3 M(g) - 0(22J(5) 80)(22J(5) 810/2)/(280(§)). (7.76) N(g) = 1602J(0)-4ZZDJ(0)8(g)+2220J2(5)—81028(5), (7.77) L(5) = (2(J(0)-J(5))(402-220J(8))-2220J2(5) + 8102J(k))/E§(5), (7.78) 73 EO(E) = (4D2-Z2DJ(E))%, (7.79) 81 = gig (23/§). (7.80) We have solved these equations numerically for a simple cubic lattice. The results are plotted in Figure 7. As with the spin S=l ferromagnet, our results closely approxi- mate Lines's results for small anisotropy. For large anisotropy the value of KTC approaches the value predicted by the MFT which is 9J(0)/2. The values of Z2 and 81 at TC agree with the MFT values for each value of the ratio D/J, just as they did for the S=l case. 74 Figure 7. The ratio TC(D)/TC(O) versus D/J for a simple cubic spin S=3/2 ferromagnet. (a) Narath's results. (b) Anderson and Callen's results. (c) Lines's results. (d) The results of this thesis. (e) The molecular field theory results. 75 (e) LG 9’: N "‘ \ 8 r0 II (n 3 0 (r) E L l l V N. (0)°1/(0)°1 |.O 4.0 3.0 2.0 D/J |.O 8. SUMMARY AND CONCLUSIONS We have demonstrated two significant ideas concern- ing the RPA decoupling of a magnetically ordered system which contains crystal field anisotrOpy. First in Chapter 5 we presented a convenient set of spin Operators which are appropriate for the solving of a system of non-interacting spins in a crystalline field environment. The inclusion of exchange coupling into the Hamiltonian then introduced higher order Green functions corresponding to the increased complexity of the system. These higher order functions, though,‘were only higher order exchange Green functions, and not higher order aniso- tropy Green functions. Hence, unlike previous authors(7’8'l3) we never had to approximate any of the anisotrOpy functions. The problem then was reduced to the decoupling of the ex- change Green functions. Second we found a decoupling approximation which is appropriate for each of the exchange Green functions, and which reproduced completely the results of the usual RPA for D=O. Therefore, we have successfully obtained re- sults from a system which contains crystal field anisotrOpy within the confines of the RPA. We have demonstrated that this technique is valid for any size spin system, 76 77 and also for antiferromagnetic as well as ferromagnetic lattices. The decoupling scheme we have presented here does remarkably well in comparison with the earlier theories on this subject. The transition temperature Tc(D) as a function of the crystal field parameter D is finite for all values of D. Our particular RPA decoupling shows that TC depends more strongly on D, for low anisotropy, than does the MFT. For anisotropies of the order of the ex- change, J, or smaller our results are equivalent to those of Lines. At the transition temperature the moments of S2 as predicted by our RPA calculation and as predicted by the MFT are the same. The temperature dependence of the moments, below Tc’ shows a greater variation with changes in D/J than does the MFT prediction. For very large anisotropies our results approach those of the MFT, and at D/J = m both theories agree exactly. This is to be contrasted with the earlier single excitation schemes which all predict TC = m at D/J = m. In Chapter 7 we presented the numerical solutions for a spin S=l ferromagnetic, a spin S=l antiferromagnet, and a spin S=3/2 ferromagnet. At the present time there is no insulating magnetic salt with S=l, or 3/2 for which the parameters D and J are known accurately. Therefore, there is no way to directly check.our calculations with experiments. Then, at best, the calculations presented here can only infer one or the other of the parameters 78 by using the known value of Tc. It is generally not pos- sible to infer these parameters from the sublattice magneti- zation curves since large changes.in D really produce relatively small changes in S(T). REFERENCES 10. 11. 12. 13. 14. 15. 16. 17. REFERENCES K. W. H. Stevens in Magnetism, edited by G. T. Rado and H. Suhl, (Academic Press, New York and London, 1964), Vol. I. N. N. Bogolyubov and S. V. Tyablikov, Dokl. Akad. Nauk SSSR 126, 53 (l959)[English transl.: Soviet Phys.--Dok1ady 4, 589 (1959)]. S. V. Tyablikov, Ukrain. Mat. Zh. 11, 287 (1959). R. A. Tahir-Kheli and D. ter Haar, Phys. Rev. 127, 88 (1962). H. B. Callen, Phys. Rev. 30, 890 (1963). R. A. Tahir-Kheli, Phys. Rev. 132, 689 (1963). F. B. Anderson and H. B. Callen, Phys. Rev. 136, A1068 (1964). M. E. Lines, Phys. Rev. 135, A1336 (1964). R. E. Mills, R. P. Kenan, and F. J. Milford, J. Appl. Phys. 36, 1131 (1965). D. N. Zubarev, Usp. Fiz. Nauk 71, 71 (1960) [English I transl.: Soviet Phys.--Uspekhi- 320 (1960)]. F. J. Dyson, Phys. Rev. 102, 1217, 1230 (1956). G. S. Rushbrooke and P. J. Wood, Mol. Phys. 1, 257 (1958). A. Narath, Phys. Rev. 140, A854 (1965). M. E. Lines, Phys. Rev. 156, 534 (1967). H. B. Callen and S. Shtrikman, Solid State Commun. 3, 5 (1965). G. S. Rushbrooke and P. J. Wood, Mol. Phys. 6, 409 (1963). P. W. Anderson, Phys. Rev. 86, 694 (1952). 79 18. 19. 20. 80 T. Murao and T. Matsubara, J. Phys. Soc. Japan 25, 352 (1968). J. H. Hetherington (unpublished). Copies of the rou- tine ANCR are available upon request. Handbook of Mathematical Functions, edited by M. Abramowitz and I. Stegan (NatiEnal Bureau of Standards, Washington, D. C., 1966) p. 886. APPENDIX A APPENDIX A COMMENTS ON MURAO AND MATSUBARA'S DECOUPLING APPROXIMATIONS (18) attempted to solve the spin Murao and Matsubara S=l ferromagnet by a technique similar to ours. Their results, however, are unsatisfactory in a number of import- ant respects. First, they predict a ratio of TC(D)/TC(O) which is smaller than the MFT prediction for each value of D. Since the MFT is known(l4) to underestimate the effect Of the crystal field anisotrOpy, especially when D is small, Murao and Matsubara's results are an even further underestimation. Second, when D is set equal to zero their results do not reduce to the RPA, but to the MFT. Both of these defects in their theory resulted from the particular Green functions they chose to decouple. For the Operators that appear on the left hand side of the Green function, they chose 5; = A; and A:’ just as we had. In our problem we had fixed the Operator on the right hand side of the Green function to be SQ. Murao and Matsubara chose the adjoints of the operators A; and AE to be inserted on the right hand side of the Green function. Hence they obtained twice as many equations 81 82 of motion as we did because they not only had our set of equations but also the set of equations with Sh replaced by the adjoint of Ag. turned out to be equivalent, leaving three independent Two Of their-four equations equations to be decoupled. But they therefore over— determined the problem as only two moments were unknown. Their resolution of this overdeterminancy was to drop the one equation which they felt, for heuristic reasons, neglected the important correlations. The remaining two equations produced the inadequacies mentioned above. Murao and Matsubara's generation of too many equa- tions of motion points out an important fact about the decoupling of Green function equations. Suppose we wish to approximate a Green function of the form, <§ (A.1) where A, B, and C are spin operators. If A is the only spin Operator with a finite ensemble average, and if one believes A and B to be relatively uncorrelated, then one is tempted to try the following, <>= <>. (A.2) The question now is whether this is automatically a valid thing to do for all possibilities for the Operator C. The point we wish to make here is that (A.2) could be a very good approximation for some choices of the Operator C, and quite bad for some other choices for C. 83 Since Murao and Matsubara performed approximations of the form of (A.2) for two different choices for C, they had no guarantee that both would be equally valid approx- imations. Clearly their two choices were mutually con- flicting, since they did Obtain thermodynamic incon- sistencies. All of the approximations that we made in (5.28)- (5.30) are internally consistent with each other since h ator C. This fact alone does not insure that we have we have only made the single choice of S for the Oper- decoupled correctly, it only says that we have been con- sistent. The final justification.for the choice C=S£ is that it does give physically meaningful results. APPENDIX B APPENDIX B THE WAVEVECTOR SUMS All of the wavevector sums <>k were performed by replacing the sum by an integral, <>k = fi- X + :1;- fffdkxdkydkz, (13.1) ~ 5 ~ lst Zone where N is the number of points in the first zone and vk is the volume of the first zone. For the antiferro- magnet, the zone referred to here is the first zone of the sublattice. The replacement in (B.1) is valid since N is Of the order of Avagadro's Number. The variable k appears in all the sums as the argu- ment of the Fourier transform of the exchange integrals, J(k), and nowhere else. For a simple cubic lattice with nearest neighbor exchange only, J(k) is as follows: J(k) = J(0)(cos(kx) + cos(ky) + cos(kz))/3, (B.2) where we have chosen the lattice parameter a=l for con- venience. For a body centered cubic lattice with nearest neighbor exchange only, J(k) is, J(k) = J(O) cos(kx) cos(ky) cos (kz), (8.3) 84 85 where we have chosen the lattice parameter, for this lat- tice, to be a=2, for convenience. For an antiferromag- net, J12(k)=J(k), and Jll(k) = 0, since neighbors on the same sublattice are really next nearest neighbors. The functions defined in (B.2) and (8.3) both have cubic symmetry, that is, they are both invariant under the 48 crystallographic Operations of a cube. Since k appears only in the functions J(k), the integrands themselves are invariant under the 48 Operations of the cube. Therefore we need not perform the integrals over the whole zone, but simply over an appropriate 1/48 of it. The boundaries of the regions of integration and the volumes, Of the regions of integration we used are as follows: Simple Cubic Ferromagnet, Boundaries: k = 0, k =k , k =k , k =0, 2 x y y z x Volume: v 03/6, k Simple Cubic Antiferromagnet (sublattice is a face centered cubic lattice), Boundaries: kz=0, kx=ky, ky=kz, kx+ky+kz=30/2,kx=n, Volume: vk=n3/12, Body Centered Cubic Ferromagnet, Boundaries: k =0, k =k , k =k , k =71—k , z x y y z x y Volume: vk=TP/24, 86 Body Centered Cubic Antiferromagnet (sublattice is a simple cubic lattice), Boundaries: k =0, k =k , k =k , k =n/2, z x y y z x Volume: vk=w3/48. The integrals were performed on a CDC 6500 high (19) speed computer, using Hetherington's integration routine ANCR. This routine determines the integral by fitting the function with appropriate Newton-Cotes‘zo) quadratures. The particular choice of quadratures de- pends on the desired error and the variations of the function over the interval of integration. Four place accuracy was generally obtained with approximately a 2000 point mesh, though this figure varied greatly depend- ing on the particular integrand. Calculation time for each integral was of the order of one second.