This preview shows pages 1–2. Sign up to view the full content.
This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: Chapter 5 Magnetic Systems
c 2010 by Harvey Gould and Jan Tobochnik 12 May 2010 We apply the general formalism of statistical mechanics developed in Chapter 4 to the Ising model, a model for which the interactions between the magnetic moments are important. We will ﬁnd that these interactions lead to a wide range of interesting phenomena, including the existence of phase transitions. Computer simulations will be used extensively and a simple, but powerful approximation method known as meanﬁeld theory will be introduced. 5.1 Paramagnetism The most familiar magnetic system in our everyday experience is probably the magnet on a refrigerator door. This magnet likely consists of iron ions localized on sites of a lattice with conduction electrons that are free to move throughout the crystal. The iron ions each have a magnetic moment and, due to a complicated interaction with each other and with the conduction electrons, they tend to line up with each other. At suﬃciently low temperatures, the moments can be aligned by an external magnetic ﬁeld and produce a net magnetic moment or magnetization which remains even if the magnetic ﬁeld is removed. Materials that retain a nonzero magnetization in zero magnetic ﬁeld are called ferromagnetic. At high enough temperatures there is enough energy to destroy the magnetization, and the iron is said to be in the paramagnetic phase. One of the key goals of this chapter is to understand the transition between the ferromagnetic and paramagnetic phases. In the simplest model of magnetism the magnetic moment can be in one of two states as discussed in Section 4.3.1. The next level of complexity is to introduce an interaction between neighboring magnetic moments. A model that includes such an interaction is discussed in Section 5.4. 229 CHAPTER 5. MAGNETIC SYSTEMS 230 5.2 Noninteracting Magnetic Moments We ﬁrst review the behavior of a system of noninteracting magnetic moments with spin 1/2 in equilibrium with a heat bath at temperature T . We discussed this system in Section 4.3.1 and in Example 4.1 using the microcanonical ensemble. The energy of interaction of a magnetic moment µ in a magnetic ﬁeld B is given by E = −µ · B = −µz B, (5.1) where µz is the component of the magnetic moment in the direction of the magnetic ﬁeld B. Because the magnetic moment has spin 1/2, it has two possible orientations. We write µz = sµ, where s = ±1. The association of the magnetic moment of a particle with its spin is an intrinsic quantum mechanical eﬀect (see Section 5.10.1). We will refer to the magnetic moment or the spin of a particle interchangeably. What would we like to know about the properties of a system of noninteracting spins? In the absence of an external magnetic ﬁeld, there is little of interest. The spins point randomly up or down because there is no preferred direction, and the mean internal energy is zero. In contrast, in the presence of an external magnetic ﬁeld, the net magnetic moment and the energy of the system are nonzero. In the following we will calculate their mean values as a function of the external magnetic ﬁeld B and the temperature T . We assume that the spins are ﬁxed on a lattice so that they are distinguishable even though the spins are intrinsically quantum mechanical. Hence the only quantum mechanical property of the system is that the spins are restricted to two values. As we will learn, the usual choice for determining the thermal properties of systems deﬁned on a lattice is the canonical ensemble. Because each spin is independent of the others and distinguishable, we can ﬁnd the partition N function for one spin, Z1 , and use the relation ZN = Z1 to obtain ZN , the partition function for N spins. (We reached a similar conclusion in Example 4.2.) We can derive this relation between N Z1 and ZN by writing the energy of the N spins as E = −µB i=1 si and expressing the partition function ZN for the N spin system as ZN =
s1 =±1 s2 =±1 ...
sN =±1 eβµB Σi=1 si eβµBs1 eβµBs2 . . . eβµBsN
sN =±1 N (5.2a) (5.2b) (5.2c) (5.2d) =
s1 =±1 s2 =±1 ... eβµBs1
s1 =±1 s2 =±1 = = eβµBs2 . . .
sN =±1 N N = Z1 . eβµBsN eβµBs1
s1 =±1 To ﬁnd Z1 we write Z1 =
s=±1 e−βµBs = eβµB (−1) + eβµB (+1) = 2 cosh βµB. (5.3) CHAPTER 5. MAGNETIC SYSTEMS Hence, the partition function for N spins is ZN = (2 cosh βµB )N . 231 (5.4) We now use the canonical ensemble formalism that we developed in Section 4.6 to ﬁnd the thermodynamic properties of the system for a given T and B . The free energy is given by F = −kT ln ZN = −N kT ln Z1 = −N kT ln(2 cosh βµB ). The mean energy E is E=− ∂ (βF ) ∂ ln ZN = = −N µB tanh βµB. ∂β ∂β (5.6) (5.5) From (5.6) we see that E → 0 as T → ∞ (β → 0). In the following we will frequently omit the mean value notation when it is clear from the context that an average is implied. Problem 5.1. Comparison of the results for two ensembles (a) Compare the result (5.6) for the mean energy E (T ) of a system of noninteracting spins in the canonical ensemble to the result that you found in Problem 4.21 for T (E ) using the microcanonical ensemble. (b) Why is it much easier to treat a system of noninteracting spins in the canonical ensemble? (c) What is the probability p that a spin is parallel to the magnetic ﬁeld B given that the system is in equilibrium with a heat bath at temperature T ? Compare your result to the result in (4.74) using the microcanonical ensemble. (d) What is the relation of the results that we have found for a system of noninteracting spins to the results obtained in Example 4.2? The heat capacity C is a measure of the change of the temperature due to the addition of energy at constant magnetic ﬁeld. The heat capacity at constant magnetic ﬁeld can be expressed as ∂E ∂E = −kβ 2 . (5.7) C= ∂T B ∂β (We will write C rather than CB .) From (5.6) and (5.7) we ﬁnd that the heat capacity of a system of N noninteracting spins is given by C = kN (βµB )2 sech2 βµB. (5.8) Note that the heat capacity is always positive, goes to zero at high T , and goes to zero as T → 0, consistent with the third law of thermodynamics. Magnetization and susceptibility . Two other macroscopic quantities of interest for magnetic systems are the mean magnetization M (in the z direction) given by
N M =µ
i=1 si , (5.9) CHAPTER 5. MAGNETIC SYSTEMS and the isothermal susceptibility χ, which is deﬁned as χ= ∂M ∂B
T 232 . (5.10) The susceptibility χ is a measure of the change of the magnetization due to the change in the external magnetic ﬁeld and is another example of a response function. We will frequently omit the factor of µ in (5.9) so that M becomes the number of spins pointing in a given direction minus the number pointing in the opposite direction. Often it is more convenient to work with the mean magnetization per spin m, an intensive variable, which is deﬁned as 1 (5.11) m = M. N As for the discussion of the heat capacity and the speciﬁc heat, the meaning of M and m will be clear from the context. We can express M and χ in terms of derivatives of ln Z by noting that the total energy can be expressed as E = E0 − M B, (5.12) where E0 is the energy of interaction of the spins with each other (the energy of the system when B = 0) and −M B is the energy of interaction of the spins with the external magnetic ﬁeld. (For noninteracting spins E0 = 0.) The form of E in (5.12) implies that we can write Z in the form Z=
s e−β (E0,s −Ms B ) , (5.13) where Ms and E0,s are the values of M and E0 in microstate s. From (5.13) we have ∂Z = ∂B βMs e−β (E0,s −Ms B ) ,
s (5.14) and hence the mean magnetization is given by M= = 1 Z Ms e−β (E0,s −Ms B )
s (5.15a) (5.15b) 1 ∂Z ∂ ln ZN = kT . βZ ∂B ∂B If we substitute the relation F = −kT ln Z , we obtain M =− ∂F . ∂B (5.16) Problem 5.2. Relation of the susceptibility to the magnetization ﬂuctuations Use considerations similar to that used to derive (5.15b) to show that the isothermal susceptibility can be written as 1 2 χ= (5.17) [M 2 − M ] . kT Note the similarity of the form (5.17) with the form (4.88) for the heat capacity CV . CHAPTER 5. MAGNETIC SYSTEMS 233 The relation of the response functions CV and χ to the equilibrium ﬂuctuations of the energy and magnetization, respectively, are special cases of a general result known as the ﬂuctuationdissipation theorem. Example 5.1. Magnetization and susceptibility of a noninteracting system of spins From (5.5) and (5.16) we ﬁnd that the mean magnetization of a system of noninteracting spins is M = N µ tanh(βµB ). The susceptibility can be calculated using (5.10) and (5.18) and is given by χ = N µ2 β sech2 (βµB ). (5.19) ♦ (5.18) Note that the arguments of the hyperbolic functions in (5.18) and (5.19) must be dimensionless and be proportional to the ratio µB/kT . Because there are only two energy scales in the system, µB , the energy of interaction of a spin with the magnetic ﬁeld, and kT , the arguments must depend only on the dimensionless ratio µB/kT . For high temperatures kT ≫ µB (βµB ≪ 1), sech(βµB ) → 1, and the leading behavior of χ is given by N µ2 (kT ≫ µB ). (5.20) χ → N µ2 β = kT The result (5.20) is known as the Curie form of the isothermal susceptibility and is commonly observed for magnetic materials at high temperatures. From (5.18) we see that M is zero at B = 0 for all T > 0, which implies that the system is paramagnetic. Because a system of noninteracting spins is paramagnetic, such a model is not applicable to materials such as iron which can have a nonzero magnetization even when the magnetic ﬁeld is zero. Ferromagnetism is due to the interactions between the spins. Problem 5.3. Thermodynamics of noninteracting spins (a) Plot the magnetization given by (5.18) and the heat capacity C given in (5.8) as a function of T for a given external magnetic ﬁeld B . Give a simple argument why C must have a broad maximum somewhere between T = 0 and T = ∞. (b) Plot the isothermal susceptibility χ versus T for ﬁxed B and describe its limiting behavior for low temperatures. (c) Calculate the entropy of a system of N noninteracting spins and discuss its limiting behavior at low (kT ≪ µB ) and high temperatures (kT ≫ µB ). Does S depend on kT and µB separately? Problem 5.4. Adiabatic demagnetization Consider a solid containing N noninteracting paramagnetic atoms whose magnetic moments can be aligned either parallel or antiparallel to the magnetic ﬁeld B . The system is in equilibrium with a heat bath at temperature T . The magnetic moment is µ = 9.274 × 10−24 J/tesla. CHAPTER 5. MAGNETIC SYSTEMS (a) If B = 4 tesla, at what temperature is 75% of the spins oriented in the +z direction? 234 (b) Assume that N = 1023 , T = 1 K, and B is increased quasistatically from 1 tesla to 10 tesla. What is the magnitude of the energy transfer from the heat bath? (c) If the system is now thermally isolated at T = 1 K and B is quasistatically decreased from 10 tesla to 1 tesla, what is the ﬁnal temperature of the system? This process is known as adiabatic demagnetization. (This problem can be solved without elaborate calculations.) 5.3 Thermodynamics of Magnetism The fundamental magnetic ﬁeld is B. However, we can usually control only the part of B due to currents in wires, and cannot directly control that part of the ﬁeld due to the magnetic dipoles in a material. Thus, we deﬁne a new ﬁeld H by H= 1 M B− , µ0 V (5.21) where M is the magnetization and V is the volume of the system. In this section we use V instead of N to make contact with standard notation in electromagnetism. Our goal in this section is to ﬁnd the magnetic equivalent of the thermodynamic relation dW = −P dV in terms of H, which we can control, and M, which we can measure. To gain insight into how to do so we consider a solenoid of length L and n turns per unit length with a magnetic material inside. When a current I ﬂows in the solenoid, there is an emf E generated in the solenoid wires. The power or rate at which work is done on the magnetic substance is −E I . By Faraday’s law we know that E =− dΦ dB = −ALn , dt dt (5.22) where the crosssectional area of the solenoid is A, the magnetic ﬂux through each turn is Φ = BA, and there are Ln turns. The work done on the system is dW = −E Idt = ALnIdB. (5.23) Ampere’s law can be used to ﬁnd that the ﬁeld H within the ideal solenoid is uniform and is given by H = nI. (5.24) Hence, (5.23) becomes dW = ALHdB = V HdB. We use (5.21) to express (5.25) as dW = µ0 V HdH + µ0 HdM. (5.26) (5.25) The ﬁrst term on the righthand side of (5.26) refers only to the ﬁeld energy, which would be there even if there were no magnetic material inside the solenoid. Thus, for the purpose of understanding the thermodynamics of the magnetic material, we can neglect the ﬁrst term and write dW = µ0 HdM. (5.27) CHAPTER 5. MAGNETIC SYSTEMS The form of (5.27) leads us to introduce the magnetic free energy G(T, M ), given by dG(T, M ) = −SdT + µ0 HdM. 235 (5.28) We use the notation G for the free energy as a function of T and M and reserve F for the free energy F (T, H ). We deﬁne F = G − µ0 HM, (5.29) and ﬁnd dF (T, H ) = dG − µ0 HdM − µ0 M dH = −SdT + µ0 HdM − µ0 HdM − µ0 M dH = −SdT − µ0 M dH. µ0 M = − Thus, we have (5.30a) (5.30b) (5.30c) ∂F . ∂H (5.31) The factor of µ0 is usually incorporated into H , so that we will usually write F (T, H ) = G(T, M ) − HM, (5.32) as well as dW = HdM and dG = −SdT + HdM . Similarly, we will write dF = −SdT − M dH , and thus ∂F M =− . (5.33) ∂H T We also write χ= ∂M ∂H
T . (5.34) The free energy F (T, H ) is frequently more useful because the quantities T and H are the easiest to control experimentally as well as in computer simulations. 5.4 The Ising Model As we saw in Section 5.1 the absence of interactions between spins implies that the system can only be paramagnetic. The most important and simplest system that exhibits a phase transition is the Ising model.1 The model was proposed by Wilhelm Lenz (1888–1957) in 1920 and was solved exactly for one dimension by his student Ernst Ising2 (1900–1998) in 1925. Ising was disappointed because the onedimensional model does not have a phase transition. Lars Onsager (1903–1976) solved the Ising model exactly in 1944 for two dimensions in the absence of an external magnetic ﬁeld and showed that there was a phase transition in two dimensions.3
1 Each year hundreds of papers are published that apply the Ising model to problems in ﬁelds as diverse as neural networks, protein folding, biological membranes, and social behavior. For this reason the Ising model is sometimes known as the “fruit ﬂy” of statistical mechanics. 2 A biographical note about Ernst Ising can be found at <www.bradley.edu/las/phy/personnel/ising.html> . 3 The model is sometimes known as the LenzIsing model. The history of the Ising model is discussed by Brush (1967). CHAPTER 5. MAGNETIC SYSTEMS 236 –J +J Figure 5.1: Two nearest neighbor spins (in any dimension) have an interaction energy −J if they are parallel and interaction energy +J if they are antiparallel. In the Ising model the spin at every site is either up (+1) or down (−1). Unless otherwise stated, the interaction is between nearest neighbors only and is given by −J if the spins are parallel and +J if the spins are antiparallel. The total energy can be expressed in the form4
N N E = −J i,j =nn(i) si sj − H si
i=1 (Ising model), (5.35) where si = ±1 and J is known as the exchange constant. We will assume that J > 0 unless otherwise stated and that the external magnetic ﬁeld is in the up or positive z direction. In the following we will refer to s as the spin.5 The ﬁrst sum in (5.35) is over all pairs of spins that are nearest neighbors. The interaction between two nearest neighbor spins is counted only once. A factor of µ has been incorporated into H , which we will refer to as the magnetic ﬁeld. In the same spirit the magnetization becomes the net number of positive spins, that is, the number of up spins minus the number of down spins. Because the number of spins is ﬁxed, we will choose the canonical ensemble and evaluate the partition function. In spite of the apparent simplicity of the Ising model it is possible to obtain exact solutions only in one dimension and in two dimensions in the absence of a magnetic ﬁeld.6 In other cases we need to use various approximation methods and computer simulations. There is no general recipe for how to perform the sums and integrals needed to calculate thermodynamic quantities. 5.5 The Ising Chain In the following we obtain an exact solution of the onedimensional Ising model and introduce an additional physical quantity of interest.
4 If we interpret the spin as an operator, then the energy is really a Hamiltonian. The distinction is unimportant here. 5 Because the spin S is a quantum mechanical object, we might expect that the commutator of the spin operator ˆ with the Hamiltonian is nonzero. However, because the Ising model retains only the component of the spin along ˆ the direction of the magnetic ﬁeld, the commutator of the spin S with the Hamiltonian is zero, and we can treat the spins in the Ising model as if they were classical. 6 It has been shown that the threedimensional Ising model (and the twodimensional Ising model with nearest neighbor and next nearest neighbor interactions) is computationally intractable and falls into the same class as other problems such as the traveling salesman problem. See <www.sandia.gov/LabNews/LN042100/sorin_story.html> and <www.siam.org/siamnews/0700/ising.pdf> . The Ising model is of interest to computer scientists in part for this reason. CHAPTER 5. MAGNETIC SYSTEMS 237 (a) (b) Figure 5.2: (a) Example of free boundary conditions for N = 9 spins. The spins at each end interact with only one spin. In contrast, all the other spins interact with two spins. (b) Example of toroidal boundary conditions. The N th spin interacts with the ﬁrst spin so that the chain forms a ring. As a result, all the spins have the same number of neighbors and the chain does not have a surface. 5.5.1 Exact enumeration As we mentioned, the canonical ensemble is the natural choice for calculating the thermodynamic N properties of systems deﬁned on a lattice. Because the spins are interacting, the relation ZN = Z1 is not applicable. The calculation of the partition function ZN is straightforward in principle. The goal is to enumerate all the microstates of the system and their corresponding energies, calculate ZN for ﬁnite N , and then take the limit N → ∞. The diﬃculty is that the total number of states, 2N , is too many for N ≫ 1. However, for the onedimensional Ising model (Ising chain) we can calculate ZN for small N and easily see how to generalize to arbitrary N . For a ﬁnite chain we need to specify the boundary conditions. One possibility is to choose free ends so that the spin at each end has only one neighbor instead of two [see Figure 5.2(a)]. Another choice is toroidal boundary conditions as shown in Figure 5.2(b). This choice implies that the N th spin is connected to the ﬁrst spin so that the chain forms a ring. In this case every spin is equivalent, and there is no boundary or surface. The choice of boundary conditions does not matter in the thermodynamic limit, N → ∞. In the absence of an external magnetic ﬁeld it is more convenient to choose free boundary conditions when calculating Z directly. (We will choose toroidal boundary conditions when doing simulations.) The energy of the Ising chain in the absence of an external magnetic ﬁeld with free boundary conditions is given explicitly by
N −1 E = −J si si+1
i=1 (free boundary conditions). (5.36) We begin by calculating the partition function for two spins. There are four possible states: both spins up with energy −J , both spins down with energy −J , and two states with one spin up CHAPTER 5. MAGNETIC SYSTEMS 238 –J –J +J +J Figure 5.3: The four possible microstates of the N = 2 Ising chain. and one spin down with energy +J (see Figure 5.3). Thus Z2 is given by Z2 = 2eβJ + 2e−βJ = 4 cosh βJ. In the same way we can enumerate the eight microstates for N = 3. We ﬁnd that Z3 = 2 e2βJ + 4 + 2 e−2βJ = 2(eβJ + e−βJ )2 = (e
βJ (5.37) (5.38a) (5.38b) +e −βJ )Z2 = (2 cosh βJ )Z2 . The relation (5.38b) between Z3 and Z2 suggests a general relation between ZN and ZN −1 : ZN = (2 cosh βJ )ZN −1 = 2 2 cosh βJ
N −1 . (5.39) We can derive the recursion relation (5.39) directly by writing ZN for the Ising chain in the form P N −1 ZN = ··· (5.40) eβJ i=1 si si+1 .
s1 =±1 sN =±1 The sum over the two possible states for each spin yields 2N microstates. To understand the meaning of the sums in (5.40), we write (5.40) for N = 3: Z3 =
s1 =±1 s2 =±1 s3 =±1 eβJs1 s2 +βJs2 s3 . (5.41) The sum over s3 can be done independently of s1 and s2 , and we have Z3 =
s1 =±1 s2 =±1 eβJs1 s2 eβJs2 + e−βJs2 eβJs1 s2 2 cosh βJs2 = 2 cosh βJ
s1 =±1 s2 =±1 s1 =±1 s2 =±1 (5.42a) eβJs1 s2 . (5.42b) = We have used the fact that the cosh function is even and hence cosh βJs2 = cosh βJ , independently of the sign of s2 . The sum over s1 and s2 in (5.42b) is Z2 . Thus Z3 is given by Z3 = (2 cosh βJ )Z2 , in agreement with (5.38b). The analysis of (5.40) for ZN proceeds similarly. We note that spin sN occurs only once in the exponential, and we have, independently of the value of sN −1 , eβJsN −1 sN = 2 cosh βJ.
sN =±1 (5.43) (5.44) CHAPTER 5. MAGNETIC SYSTEMS
0.5 239 0.4 C Nk 0.3 0.2 0.1 0.0 0 1 2 3 4 5 6 7 8 kT/J
Figure 5.4: The temperature dependence of the speciﬁc heat (in units of k ) of an Ising chain in the absence of an external magnetic ﬁeld. At what value of kT /J does C exhibit a maximum? Hence we can write ZN as ZN = (2 cosh βJ )ZN −1 . We can continue this process to ﬁnd ZN = (2 cosh βJ )2 ZN −2 = (2 cosh βJ ) ZN −3 . . . = (2 cosh βJ )N −1 Z1 = 2(2 cosh βJ )N −1 ,
3 (5.45) (5.46a) (5.46b) (5.46c) where we have used the fact that Z1 = s1 =±1 1 = 2. No Boltzmann factor appears in Z1 because there are no interactions with one spin. We can use the general result (5.39) for ZN to ﬁnd the Helmholtz free energy: F = −kT ln ZN = −kT ln 2 + (N − 1) ln(2 cosh βJ ) . (5.47) In the thermodynamic limit N → ∞, the term proportional to N in (5.47) dominates, and we have the desired result: F = −N kT ln 2 cosh βJ . (5.48) Problem 5.5. Exact enumeration Enumerate the 2N microstates for the N = 4 Ising chain and ﬁnd the corresponding contributions to Z4 for free boundary conditions. Then show that Z4 and Z3 satisfy the recursion relation (5.45) for free boundary conditions. CHAPTER 5. MAGNETIC SYSTEMS Problem 5.6. Thermodynamics of the Ising chain (a) What is the ground state of the Ising chain? 240 (b) What is the entropy S in the limits T → 0 and T → ∞? The answers can be found without doing an explicit calculation. (c) Use (5.48) for the free energy F to verify the following results for the entropy S , the mean energy E , and the heat capacity C of the Ising chain: S = N k ln(e2βJ + 1) − E = −N J tanh βJ,
2 2 2βJ , 1 + e−2βJ (5.49) (5.50) (5.51) C = N k (βJ ) (sech βJ ) . Verify that the results in (5.49)–(5.51) reduce to the appropriate behavior for low and high temperatures. (d) A plot of the T dependence of the heat capacity in the absence of a magnetic ﬁeld is given in Figure 5.4. Explain why it has a maximum. 5.5.2 Spinspin correlation function We can gain further insight into the properties of the Ising model by calculating the spinspin correlation function G(r) deﬁned as G(r) = sk sk+r − sk sk+r . (5.52) Because the average of sk is independent of the choice of the site k (for toroidal boundary conditions) and equals m = M/N , G(r) can be written as G(r) = sk sk+r − m2 . (5.53) The average is over all microstates. Because all lattice sites are equivalent, G(r) is independent of the choice of k and depends only on the separation r (for a given T and H ), where r is the separation between the two spins in units of the lattice constant. Note that G(r = 0) = m2 − m2 ∝ χ [see (5.17)]. The spinspin correlation function G(r) is a measure of the degree to which a spin at one site is correlated with a spin at another site. If the spins are not correlated, then G(r) = 0. At high temperatures the interaction between spins is unimportant, and hence the spins are randomly oriented in the absence of an external magnetic ﬁeld. Thus in the limit kT ≫ J , we expect that G(r) → 0 for any r. For ﬁxed T and H , we expect that, if spin k is up, then the two adjacent spins will have a greater probability of being up than down. For spins further away from spin k , we expect that the probability that spin k + r is up or correlated will decrease. Hence, we expect that G(r) → 0 as r → ∞. CHAPTER 5. MAGNETIC SYSTEMS
1.0 241 0.8 G(r)
0.6 0.4 0.2 0.0 0 25 50 75 100 r 125 150 Figure 5.5: Plot of the spinspin correlation function G(r) as given by (5.54) for the Ising chain for βJ = 2. Problem 5.7. Calculation of G(r) for three spins Consider an Ising chain of N = 3 spins with free boundary conditions in equilibrium with a heat bath at temperature T and in zero magnetic ﬁeld. Enumerate the 23 microstates and calculate G(r = 1) and G(r = 2) for k = 1, the ﬁrst spin on the left. We show in the following that G(r) can be calculated exactly for the Ising chain and is given by G(r) = tanh βJ .
r (5.54) A plot of G(r) for βJ = 2 is shown in Figure 5.5. Note that G(r) → 0 for r ≫ 1 as expected. We also see from Figure 5.5 that we can associate a length with the decrease of G(r). We will deﬁne the correlation length ξ by writing G(r) in the form G(r) = e−r/ξ where ξ=− 1 . ln(tanh βJ ) (5.56) (r ≫ 1), (5.55) At low temperatures, tanh βJ ≈ 1 − 2e−2βJ , and ln tanh βJ ≈ −2e−2βJ . Hence ξ= 1 2βJ e 2 (βJ ≫ 1). (5.58) (5.57) The correlation length is a measure of the distance over which the spins are correlated. From (5.58) we see that the correlation length becomes very large for low temperatures (βJ ≫ 1). CHAPTER 5. MAGNETIC SYSTEMS 242 Problem 5.8. What is the maximum value of tanh βJ ? Show that for ﬁnite values of βJ , G(r) given by (5.54) decays with increasing r.
∗ General calculation of G(r ) in one dimension. To calculate G(r) in the absence of an external magnetic ﬁeld we assume free boundary conditions. It is useful to generalize the Ising model and assume that the magnitude of each of the nearest neighbor interactions is arbitrary so that the total energy E is given by
N −1 E=− Ji si si+1 ,
i=1 (5.59) where Ji is the interaction energy between spin i and spin i + 1. At the end of the calculation we will set Ji = J . We will ﬁnd in Section 5.5.4 that m = 0 for T > 0 for the onedimensional Ising model. Hence, we can write G(r) = sk sk+r . For the form (5.59) of the energy, sk sk+r is given by s k s k +r = where
N −1 1 ZN s N −1 1 =± 1 ··· sk sk+r exp
sN =±1 i=1 βJi si si+1 , (5.60) ZN = 2
i=1 2 cosh βJi . (5.61) The righthand side of (5.60) is the value of the product of two spins separated by a distance r in a particular microstate times the probability of that microstate. We now use a trick similar to that used in Section 3.5 and the Appendix to calculate various sums and integrals. If we take the derivative of the exponential in (5.60) with respect to Jk , we bring down a factor of βsk sk+1 . Hence, the spinspin correlation function G(r = 1) = sk sk+1 for the Ising model with Ji = J can be expressed as sk sk+1 = 1 ZN s
N −1 1 =± 1 ··· sk sk+1 exp
sN =±1 i=1 N −1 βJi si si+1 (5.62a) = = = 11∂ ZN β ∂Jk s1 =±1 ··· exp
sN =±1 i=1 βJi si si+1 (5.62b) (5.62c) 1 1 ∂ZN (J1 , · · · , JN −1 ) ZN β ∂Jk sinh βJ = tanh βJ, cosh βJ J i =J (5.62d) where we have used the form (5.61) for ZN . To obtain G(r = 2), we use the fact that s2 +1 = 1 to k write sk sk+2 = sk (sk+1 sk+1 )sk+2 = (sk sk+1 )(sk+1 sk+2 ). We write G(r = 2) = = 1 ZN
N −1 sk sk+1 sk+1 sk+2 exp
{sj } i=1 βJi si si+1 (5.63a) (5.63b) 1 1 ∂ 2 ZN (J1 , · · · , JN −1 ) = [tanh βJ ]2 . ZN β 2 ∂Jk ∂Jk+1 CHAPTER 5. MAGNETIC SYSTEMS 243 The method used to obtain G(r = 1) and G(r = 2) can be generalized to arbitrary r. We write G(r) = and use (5.61) for ZN to ﬁnd that G(r) = tanh βJk tanh βJk+1 · · · tanh βJk+r−1
r ∂ ∂ 11∂ ··· ZN , ZN β r ∂Jk Jk+1 Jk+r−1 (5.64) (5.65a) (5.65b) =
k=1 tanh βJk+r−1 . For a uniform interaction, Ji = J , and (5.65b) reduces to the result for G(r) in (5.54). 5.5.3 Simulations of the Ising chain Although we have found an exact solution for the onedimensional Ising model in the absence of an external magnetic ﬁeld, we can gain additional physical insight by doing simulations. As we will see, simulations are essential for the Ising model in higher dimensions. As we discussed in Section 4.11, page 215, the Metropolis algorithm is the simplest and most common Monte Carlo algorithm for a system in equilibrium with a heat bath at temperature T . In the context of the Ising model, the Metropolis algorithm can be implemented as follows: 1. Choose an initial microstate of N spins. The two most common initial states are the ground state with all spins parallel or the T = ∞ state where each spin is chosen to be ±1 at random. 2. Choose a spin at random and make a trial ﬂip. Compute the change in energy of the system, ∆E , corresponding to the ﬂip. The calculation is straightforward because the change in energy is determined by only the nearest neighbor spins. If ∆E < 0, then accept the change. If ∆E > 0, accept the change with probability p = e−β ∆E . To do so, generate a random number r uniformly distributed in the unit interval. If r ≤ p, accept the new microstate; otherwise, retain the previous microstate. 3. Repeat step 2 many times choosing spins at random. 4. Compute averages of the quantities of interest such as E , M , C , and χ after the system has reached equilibrium. In the following problem we explore some of the qualitative properties of the Ising chain. Problem 5.9. Computer simulation of the Ising chain Use Program Ising1D to simulate the onedimensional Ising model. It is convenient to measure the temperature in units such that J/k = 1. For example, a temperature of T = 2 means that T = 2J/k . The “time” is measured in terms of Monte Carlo steps per spin (mcs), where in one Monte Carlo step per spin, N spins are chosen at random for trial changes. (On the average each spin will be chosen equally, but during any ﬁnite interval, some spins might be chosen more than others.) Choose H = 0. The thermodynamic quantities of interest for the Ising model include the mean energy E , the heat capacity C , and the isothermal susceptibility χ. CHAPTER 5. MAGNETIC SYSTEMS 244 (a) Determine the heat capacity C and susceptibility χ for diﬀerent temperatures, and discuss the qualitative temperature dependence of χ and C . Choose N ≥ 200. (b) Why is the mean value of the magnetization of little interest for the onedimensional Ising model? Why does the simulation usually give M = 0? (c) Estimate the mean size of the domains at T = 1.0 and T = 0.5. By how much does the mean size of the domains increase when T is decreased? Compare your estimates with the correlation length given by (5.56). What is the qualitative temperature dependence of the mean domain size? (d) Why does the Metropolis algorithm become ineﬃcient at low temperatures? 5.5.4 *Transfer matrix So far we have considered the Ising chain only in zero external magnetic ﬁeld. The solution for nonzero magnetic ﬁeld requires a diﬀerent approach. We now apply the transfer matrix method to solve for the thermodynamic properties of the Ising chain in nonzero magnetic ﬁeld. The transfer matrix method is powerful and can be applied to various magnetic systems and to seemingly unrelated quantum mechanical systems. The transfer matrix method also is of historical interest because it led to the exact solution of the twodimensional Ising model in the absence of a magnetic ﬁeld. A background in matrix algebra is important for understanding the following discussion. To apply the transfer matrix method to the onedimensional Ising model, it is necessary to adopt toroidal boundary conditions so that the chain becomes a ring with sN +1 = s1 . This boundary condition enables us to write the energy as E = −J 1 si si+1 − H 2 i=1
N N (si + si+1 ) (toroidal boundary conditions).
i=1 (5.66) The use of toroidal boundary conditions implies that each spin is equivalent. The transfer matrix T is deﬁned by its four matrix elements, which are given by Ts,s′ = eβ [Jss + 2 H (s+s )] . The explicit form of the matrix elements is T++ = eβ (J +H ) , T−− = e
β (J −H )
′ 1 ′ (5.67) (5.68a) (5.68b) , (5.68c) ,
−βJ T − + = T +− = e or T= T++ T+− T−+ T−− = eβ (J +H ) e−βJ . e−βJ e β (J −H ) (5.69) The deﬁnition (5.67) of T allows us to write ZN in the form ZN (T, H ) =
s1 s2 ··· sN Ts1 ,s2 Ts2 ,s3 · · · TsN ,s1 . (5.70) CHAPTER 5. MAGNETIC SYSTEMS The form of (5.70) is suggestive of the interpretation of T as a transfer function. 245 Problem 5.10. Transfer matrix method in zero magnetic ﬁeld Show that the partition function for a system of N = 3 spins with toroidal boundary conditions can be expressed as the trace (the sum of the diagonal elements) of the product of three matrices: eβJ e−βJ e−βJ eβJ eβJ e−βJ e−βJ eβJ eβJ e−βJ e−βJ eβJ . (5.71) The rule for matrix multiplication that we need for the transfer matrix method is (T2 )s1 ,s3 =
s2 Ts1 ,s2 Ts2 ,s3 , (5.72) or (T)2 = If we multiply N matrices, we obtain (TN )s1 ,sN +1 =
s2 T++ T++ T − + T +− T +− T − + . T−− T−− (5.73) s3 ··· sN Ts1 ,s2 Ts2 ,s3 · · · TsN ,sN +1 . (5.74) This result is very close to the form of ZN in (5.70). To make it identical, we use toroidal boundary conditions and set sN +1 = s1 , and sum over s1 : (TN )s1 ,s1 =
s1 s1 s2 s3 ··· sN Ts1 ,s2 Ts2 ,s3 · · · TsN ,s1 = ZN . (5.75) Because we have N s1 (T )s1 ,s1 is the deﬁnition of the trace (the sum of the diagonal elements) of (TN ), ZN = Tr (TN ). (5.76) The fact that ZN is the trace of the N th power of a matrix is a consequence of our use of toroidal boundary conditions. Because the trace of a matrix is independent of the representation of the matrix, the trace in (5.76) may be evaluated by bringing T into diagonal form: T= λ+ 0 . 0 λ− (5.77) The matrix TN is diagonal with the diagonal matrix elements λN , λN . In the diagonal represen− + tation for T in (5.77), we have ZN = Tr (TN ) = λN + λN , (5.78) + − where λ+ and λ− are the eigenvalues of T. The eigenvalues λ± are given by the solution of the determinant equation e β ( J +H ) − λ e−βJ = 0. −βJ β (J −H ) e e −λ (5.79) CHAPTER 5. MAGNETIC SYSTEMS The roots of (5.79) are λ± = eβJ cosh βH ± e−2βJ + e2βJ sinh2 βH
1/2 246 . (5.80) It is easy to show that λ+ > λ− for all β and H , and consequently (λ−/λ+ )N → 0 as N → ∞. In the thermodynamic limit N → ∞ we obtain from (5.78) and (5.80) 1 λ− ln ZN (T, H ) = ln λ+ + ln 1 + N →∞ N λ+ lim and the free energy per spin is given by f (T, H ) = 1 F (T, H ) = −kT ln eβJ cosh βH + e2βJ sinh2 βH + e−2βJ N
1/2 N = ln λ+ , (5.81) . (5.82) We can use (5.82) and (5.31) and some algebraic manipulations to ﬁnd the magnetization per spin m at nonzero T and H : m=− ∂f sinh βH . = ∂H (sinh2 βH + e−4βJ )1/2 (5.83) A system is paramagnetic if m = 0 only when H = 0, and is ferromagnetic if m = 0 when H = 0. From (5.83) we see that m = 0 for H = 0 because sinh x ≈ x for small x. Thus for H = 0, sinh βH = 0 and thus m = 0. The onedimensional Ising model becomes a ferromagnet only at T = 0 where e−4βJ → 0, and thus from (5.83) m → 1 at T = 0. Problem 5.11. Isothermal susceptibility of the Ising chain More insight into the properties of the Ising chain can be found by understanding the temperature dependence of the isothermal susceptibility χ. (a) Calculate χ using (5.83). (b) What is the limiting behavior of χ in the limit T → 0 for H > 0? (c) Show that the limiting behavior of the zeroﬁeld susceptibility in the limit T → 0 is χ ∼ e2βJ . (The zeroﬁeld susceptibility is found by calculating the susceptibility for H = 0 and then taking the limit H → 0 before other limits such as T → 0 are taken.) Express the limiting behavior in terms of the correlation length ξ . Why does χ diverge as T → 0? Because the zeroﬁeld susceptibility diverges as T → 0, the ﬂuctuations of the magnetization also diverge in this limit. As we will see in Section 5.6, the divergence of the magnetization ﬂuctuations is one of the characteristics of the critical point of the Ising model. That is, the phase transition from a paramagnet (m = 0 for H = 0) to a ferromagnet (m = 0 for H = 0) occurs at zero temperature for the onedimensional Ising model. We will see that the critical point occurs at T > 0 for the Ising model in two and higher dimensions. CHAPTER 5. MAGNETIC SYSTEMS 247 (a) (b) Figure 5.6: A domain wall in one dimension for a system of N = 8 spins with free boundary conditions. In (a) the energy of the system is E = −5J (H = 0). The energy cost for forming a domain wall is 2J (recall that the ground state energy is −7J ). In (b) the domain wall has moved with no cost in energy. 5.5.5 Absence of a phase transition in one dimension We found by direct calculations that the onedimensional Ising model does not have a phase transition for T > 0. We now argue that a phase transition in one dimension is impossible if the interaction is shortrange, that is, if a given spin interacts with only a ﬁnite number of spins. At T = 0 the energy is a minimum with E = −(N − 1)J (for free boundary conditions), and the entropy S = 0.7 Consider all the excitations at T > 0 obtained by ﬂipping all the spins to the right of some site [see Figure 5.6(a)]. The energy cost of creating such a domain wall is 2J . Because there are N − 1 sites where the domain wall may be placed, the entropy increases by ∆S = k ln(N − 1). Hence, the free energy cost associated with creating one domain wall is ∆F = 2J − kT ln(N − 1). (5.84) We see from (5.84) that for T > 0 and N → ∞, the creation of a domain wall lowers the free energy. Hence, more domain walls will be created until the spins are completely randomized and the net magnetization is zero. We conclude that M = 0 for T > 0 in the limit N → ∞. Problem 5.12. Energy cost of a single domain Compare the energy of the microstate in Figure 5.6(a) with the energy of the microstate shown in Figure 5.6(b) and discuss why the number of spins in a domain in one dimension can be changed without any energy cost. 5.6 The TwoDimensional Ising Model We ﬁrst give an argument similar to the one that was given in Section 5.5.5 to suggest the existence of a paramagnetic to ferromagnetism phase transition in the twodimensional Ising model at a nonzero temperature. We will show that the mean value of the magnetization is nonzero at low, but nonzero temperatures and in zero magnetic ﬁeld. The key diﬀerence between one and two dimensions is that in the former the existence of one domain wall allows the system to have regions of up and down spins whose size can be changed without any cost of energy. So on the average the number of up and down spins is the same. In two dimensions the existence of one domain does not make the magnetization zero. The regions of
7 The ground state for H = 0 corresponds to all spins up or all spins down. It is convenient to break this symmetry by assuming that H = 0+ and letting T → 0 before letting H → 0+ . CHAPTER 5. MAGNETIC SYSTEMS 248 (a) (b) Figure 5.7: (a) The ground state of a 5 × 5 Ising model. (b) Example of a domain wall. The energy cost of the domain is 10J , assuming free boundary conditions. down spins cannot grow at low temperature because their growth requires longer boundaries and hence more energy. From Figure 5.7 we see that the energy cost of creating a rectangular domain in two dimensions is given by 2JL (for an L × L lattice with free boundary conditions). Because the domain wall can be at any of the L columns, the entropy is at least order ln L. Hence the free energy cost of creating one domain is ∆F ∼ 2JL − T ln L. hence, we see that ∆F > 0 in the limit L → ∞. Therefore creating one domain increases the free energy and thus most of the spins will remain up, and the magnetization remains positive. Hence M > 0 for T > 0, and the system is ferromagnetic. This argument suggests why it is possible for the magnetization to be nonzero for T > 0. M becomes zero at a critical temperature Tc > 0, because there are many other ways of creating domains, thus increasing the entropy and leading to a disordered phase. 5.6.1 Onsager solution As mentioned, the twodimensional Ising model was solved exactly in zero magnetic ﬁeld for a rectangular lattice by Lars Onsager in 1944. Onsager’s calculation was the ﬁrst exact solution that exhibited a phase transition in a model with shortrange interactions. Before his calculation, some people believed that statistical mechanics was not capable of yielding a phase transition. Although Onsager’s solution is of much historical interest, the mathematical manipulations are very involved. Moreover, the manipulations are special to the Ising model and cannot be generalized to other systems. For these reasons few workers in statistical mechanics have gone through the Onsager solution in great detail. In the following, we summarize some of the results of the twodimensional solution for a square lattice. The critical temperature Tc is given by sinh 2J = 1, kTc (5.85) CHAPTER 5. MAGNETIC SYSTEMS
1.2 1.0 249 κ
0.8 0.6 0.4 0.2 0.0 0.0 0.5 1.0 1.5 2.0 2.5 βJ 3.0 3.5 4.0 Figure 5.8: Plot of the function κ deﬁned in (5.87) as a function of J/kT . or 2 kTc √ ≈ 2.269. = J ln(1 + 2) (5.86) It is convenient to express the mean energy in terms of the dimensionless parameter κ deﬁned as κ=2 sinh 2βJ . (cosh 2βJ )2 (5.87) A plot of the function κ versus βJ is given in Figure 5.8. Note that κ is zero at low and high temperatures and has a maximum of 1 at T = Tc . The exact solution for the energy E can be written in the form E = −2N J tanh 2βJ − N J where K1 (κ) =
0 sinh2 2βJ − 1 2 K1 (κ) − 1 , sinh 2βJ cosh 2βJ π dφ 1 − κ2 sin2 φ (5.88) π /2 . (5.89) K1 is known as the complete elliptic integral of the ﬁrst kind. The ﬁrst term in (5.88) is similar to the result (5.50) for the energy of the onedimensional Ising model with a doubling of the exchange interaction J for two dimensions. The second term in (5.88) vanishes at low and high temperatures (because of the term in brackets) and at T = Tc because of the vanishing of the term sinh2 2βJ − 1. The function K1 (κ) has a logarithmic singularity at T = Tc at which κ = 1. Hence, the second term behaves as (T − Tc ) ln T − Tc  in the vicinity of Tc . We conclude that E (T ) is continuous at T = Tc and at all other temperatures [see Figure 5.9(a)]. CHAPTER 5. MAGNETIC SYSTEMS
–0.4 4.0 250 –0.8 3.0 E NJ
–1.2 C Nk
2.0 –1.6 1.0 –2.0 0.0 1.0 2.0
(a) 3.0 4.0 5.0 0.0 0.0 1.0 2.0
(b) 3.0 4.0 5.0 kT/J kT/J Figure 5.9: (a) Temperature dependence of the energy of the Ising model on the square lattice according to (5.88). Note that E (T ) is a continuous function of kT /J . (b) Temperature dependence of the speciﬁc heat of the Ising model on the square lattice according to (5.90). Note the divergence of the speciﬁc heat at the critical temperature. The heat capacity can be obtained by diﬀerentiating E (T ) with respect to temperature. It can be shown after some tedious algebra that 4 C (T ) = N k (βJ coth 2βJ )2 K1 (κ) − E1 (κ) π π + (2 tanh2 2βJ − 1)K1 (κ) − (1 − tanh2 2βJ ) 2 where E1 (κ) =
0 π /2 , (5.90) dφ 1 − κ2 sin2 φ. (5.91) E1 is called the complete elliptic integral of the second kind. A plot of C (T ) is given in Figure 5.9(b). The behavior of C near Tc is given by C ≈ −N k 2 2J π kTc
2 ln 1 − T + constant Tc (T near Tc ). (5.92) An important property of the Onsager solution is that the heat capacity diverges logarithmically at T = Tc : C (T ) ∼ − ln ǫ, (5.93) where the reduced temperature diﬀerence is given by ǫ = (Tc − T )/Tc . (5.94) CHAPTER 5. MAGNETIC SYSTEMS 251 A major test of the approximate treatments that we will develop in Section 5.7 and in Chapter 9 is whether they can yield a heat capacity that diverges as in (5.93). The power law divergence of C (T ) can be written in general as C (T ) ∼ ǫ−α , (5.95) Because the divergence of C in (5.93) is logarithmic, which is slower than any power of ǫ, the critical exponent α equals zero for the twodimensional Ising model. To know whether the logarithmic divergence of the heat capacity in the Ising model at T = Tc is associated with a phase transition, we need to know if there is a spontaneous magnetization. That is, is there a range of T > 0 such that M = 0 for H = 0? (Onsager’s solution is limited to zero magnetic ﬁeld.) To calculate the spontaneous magnetization we need to calculate the derivative of the free energy with respect to H for nonzero H and then let H = 0. In 1952 C. N. Yang calculated the magnetization for T < Tc and the zeroﬁeld susceptibility.8 Yang’s exact result for the magnetization per spin can be expressed as m(T ) = 1 − [sinh 2βJ ]−4 0
1/8 (T < Tc ), (T > Tc ). (5.96) A plot of m is shown in Figure 5.10. We see that m vanishes near Tc as m ∼ ǫβ (T < Tc ), (5.97) where β is a critical exponent and should not be confused with the inverse temperature. For the twodimensional Ising model β = 1/8. The magnetization m is an example of an order parameter. For the Ising model m = 0 for T > Tc (paramagnetic phase), and m = 0 for T ≤ Tc (ferromagnetic phase). The word “order” in the magnetic context is used to denote that below Tc the spins are mostly aligned in the same direction; in contrast, the spins point randomly in both directions for T above Tc . The behavior of the zeroﬁeld susceptibility for T near Tc was found by Yang to be χ ∼ ǫ−7/4 ∼ ǫ−γ , (5.98) where γ is another critical exponent. We see that γ = 7/4 for the twodimensional Ising model. The most important results of the exact solution of the twodimensional Ising model are that the energy (and the free energy and the entropy) are continuous functions for all T , m − vanishes continuously at T = Tc , the heat capacity diverges logarithmically at T = Tc , and the zeroﬁeld susceptibility and other quantities show power law behavior which can be described by critical exponents. We say that the paramagnetic ↔ ferromagnetic transition in the twodimensional Ising model is continuous because the order parameter m vanishes continuously rather
8 C. N. Yang, “The spontaneous magnetization of a twodimensional Ising model,” Phys. Rev. 85, 808– 816 (1952). The result (5.96) was ﬁrst announced by Onsager at a conference in 1944 but not published. C. N. Yang and T. D. Lee shared the 1957 Nobel Prize in Physics for work on parity violation. See <nobelprize.org/physics/laureates/1957/> . CHAPTER 5. MAGNETIC SYSTEMS 252 1.0 m
0.8 0.6 0.4 0.2 0.0 0.0 Tc 0.5 1.0 1.5 2.0 2.5 3.0 kT/J
Figure 5.10: The temperature dependence of the spontaneous magnetization m(T ) of the twodimensional Ising model. than discontinuously. Because the transition occurs only at T = Tc and H = 0, the transition occurs at a critical point. So far we have introduced the critical exponents α, β , and γ to describe the behavior of the speciﬁc heat, magnetization, and susceptibility near the critical point. We now introduce three more critical exponents: η , ν , and δ (see Table 5.1). The notation χ ∼ ǫ−γ means that χ has a singular contribution proportional to ǫ−γ . The deﬁnitions of the critical exponents given in Table 5.1 implicitly assume that the singularities are the same whether the critical point is approached from above or below Tc . The exception is m, which is zero for T > Tc . The critical exponent δ characterizes the dependence of m on the magnetic ﬁeld at T = Tc : m ∼ H 1/15 ∼ H 1/δ We see that δ = 15 for the twodimensional Ising model. The behavior of the spinspin correlation function G(r) for T near Tc and large r is given by G(r) ∼ 1 e−r/ξ rd−2+η (r ≫ 1 and ǫ ≪ 1), (5.100) (T = Tc ). (5.99) where d is the spatial dimension and η is another critical exponent. The correlation length ξ diverges as ξ ∼ ǫ−ν . (5.101) The exact result for the critical exponent ν for the twodimensional (d = 2) Ising model is ν = 1. At T = Tc , G(r) decays as a power law for large r: G(r) = 1 rd−2+η (T = Tc , r ≫ 1). (5.102) For the twodimensional Ising model η = 1/4. The values of the various critical exponents for the Ising model in two and three dimensions are summarized in Table 5.1. CHAPTER 5. MAGNETIC SYSTEMS values d = 2 (exact) 0 (logarithmic) 1/8 7/4 15 1 1/4 253 of the exponents d = 3 meanﬁeld theory 0.113 0 (jump) 0.324 1/2 1.238 1 4.82 3 0.629 1/2 0.031 0 Quantity speciﬁc heat order parameter susceptibility equation of state (ǫ = 0) correlation length correlation function ǫ = 0 Critical behavior C ∼ ǫ −α m ∼ ǫβ χ ∼ ǫ −γ m ∼ H 1/δ ξ ∼ ǫ −ν G(r) ∼ 1/rd−2+η Table 5.1: Values of the critical exponents for the Ising model in two and three dimensions. The values of the critical exponents of the Ising model are known exactly in two dimensions and are ratios of integers. The results in three dimensions are not ratios of integers and are approximate. The exponents predicted by meanﬁeld theory are discussed in Sections 5.7, and 9.1, pages 256 and 434. There is a fundamental diﬀerence between the exponential behavior of G(r) for T = Tc in (5.100) and the power law behavior of G(r) for T = Tc in (5.102). Systems with correlation functions that decay as a power law are said to be scale invariant. That is, power laws look the same on all scales. The replacement x → ax in the function f (x) = Ax−η yields a function that is indistinguishable from f (x) except for a change in the amplitude A by the factor a−η . In contrast, this invariance does not hold for functions that decay exponentially because making the replacement x → ax in the function e−x/ξ changes the correlation length ξ by the factor a. The fact that the critical point is scale invariant is the basis for the renormalization group method (see Chapter 9). Scale invariance means that at the critical point there will be domains of spins of the same sign of all sizes. We stress that the phase transition in the Ising model is the result of the cooperative interactions between the spins. Although phase transitions are commonplace, they are remarkable from a microscopic point of view. For example, the behavior of the system changes dramatically with a small change in the temperature even though the interactions between the spins remain unchanged and shortrange. The study of phase transitions in relatively simple systems such as the Ising model has helped us begin to understand phenomena as diverse as the distribution of earthquake sizes, the shape of snowﬂakes, and the transition from a boom economy to a recession. 5.6.2 Computer simulation of the twodimensional Ising model The implementation of the Metropolis algorithm for the twodimensional Ising model proceeds as in one dimension. The only diﬀerence is that an individual spin interacts with four nearest neighbors on a square lattice rather than two nearest neighbors in one dimension. Simulations of the Ising model in two dimensions allow us to test approximate theories and determine properties that cannot be calculated analytically. We explore some of the properties of the twodimensional Ising model in Problem 5.13. Problem 5.13. Simulation of the twodimensional Ising model Use Program Ising2D to simulate the Ising model on a square lattice at a given temperature T and external magnetic ﬁeld H . (Remember that T is given in terms of J/k .) First choose N = L2 = 322 CHAPTER 5. MAGNETIC SYSTEMS and set H = 0. For simplicity, the initial orientation of the spins is all spins up. 254 (a) Choose T = 10 and run until equilibrium has been established. Is the orientation of the spins random such that the mean magnetization is approximately equal to zero? What is a typical size of a domain, a region of parallel spins? (b) Choose a low temperature such as T = 0.5. Are the spins still random or do a majority choose a preferred direction? You will notice that M ≈ 0 for suﬃcient high T and M = 0 for suﬃciently low T . Hence, there is an intermediate value of T at which M ﬁrst becomes nonzero. (c) Start at T = 4 and determine the temperature dependence of the magnetization per spin m, the zeroﬁeld susceptibility χ, the mean energy E , and the speciﬁc heat C . (Note that we have used the same notation for the speciﬁc heat and the heat capacity.) Decrease the temperature in intervals of 0.2 until T ≈ 1.6, equilibrating for at least 1000 mcs before collecting data at each value of T . Describe the qualitative temperature dependence of these quantities. Note that when the simulation is stopped, the mean magnetization and the mean of the absolute value of the magnetization is returned. At low temperatures the magnetization can sometimes ﬂip for small systems so that the value of M  is a more accurate representation of the magnetization. For the same reason the susceptibility is given by χ= 1 kT M 2 − M 
2 , (5.103) rather than by (5.17). A method for estimating the critical exponents is discussed in Problem 5.41. (d) Set T = Tc ≈ 2.269 and choose L ≥ 128. Obtain M for H = 0.01, 0.02, 0.04, 0.08, and 0.16. Make sure you equilibrate the system at each value of H before collecting data. Make a loglog plot of m versus H and estimate the critical exponent δ using (5.99). (e) Choose L = 4 and T = 2.0. Does the sign of the magnetization change during the simulation? Choose a larger value of L and observe if the sign of the magnetization changes. Will the sign of M change for L ≫ 1? Should a theoretical calculation of M yield M = 0 or M = 0 for T < Tc ? Problem 5.14. Ising antiferromagnet So far we have considered the ferromagnetic Ising model for which the energy of interaction between two nearest neighbor spins is J > 0. Hence, all spins are parallel in the ground state of the ferromagnetic Ising model. In contrast, if J < 0, nearest neighbor spins must be antiparallel to minimize their energy of interaction. (a) Sketch the ground state of the onedimensional antiferromagnetic Ising model. Then do the same for the antiferromagnetic Ising model on a square lattice. What is the value of M for the ground state of an Ising antiferromagnet?
∗ CHAPTER 5. MAGNETIC SYSTEMS 255 √3 a 2 a Figure 5.11: Each spin has six nearest neighbors on a hexagonal lattice. This lattice structure is sometimes called a triangular lattice. ? Figure 5.12: The six nearest neighbors of the central spin on a hexagonal lattice are successively antiparallel, corresponding to the lowest energy of interaction for an Ising antiferromagnet. The central spin cannot be antiparallel to all its neighbors and is said to be frustrated. (b) Use Program IsingAnitferromagnetSquareLattice to simulate the antiferromagnetic Ising model on a square lattice at various temperatures and describe its qualitative behavior. Does the system have a phase transition at T > 0? Does the value of M show evidence of a phase transition? (c) In addition to the usual thermodynamic quantities the program calculates the staggered magnetization and the staggered susceptibility. The staggered magnetization is calculated by considering the square lattice as a checkerboard with black and red sites so that each black site has four red sites as nearest neighbors and vice versa. The staggered magnetization is calculated from ci si where ci = +1 for a black site and ci = −1 for a red site. Describe the behavior of these quantities and compare them to the behavior of M and χ for the ferromagnetic Ising model. (d) *Consider the Ising antiferromagnetic model on a hexagonal lattice (see Figure 5.11), for which each spin has six nearest neighbors. The ground state in this case is not unique because of frustration (see Figure 5.12). Convince yourself that there are multiple ground states. Is the entropy zero or nonzero at T = 0?9 Use Program IsingAntiferromagnetHexagonalLattice
9 The entropy at zero temperature is S (T = 0) = 0.3383kN . See G. H. Wannier, “Antiferromagnetism. The CHAPTER 5. MAGNETIC SYSTEMS 256 to simulate the antiferromagnetic Ising model on a hexagonal lattice at various temperatures and describe its qualitative behavior. This system does not have a phase transition for T > 0. Are your results consistent with this behavior? 5.7 MeanField Theory Because it is not possible to solve the thermodynamics of the Ising model exactly in three dimensions and the twodimensional Ising model in the presence of a magnetic ﬁeld, we need to develop approximate theories. In this section we develop an approximation known as meanﬁeld theory. Meanﬁeld theories are relatively easy to treat and usually yield qualitatively correct results, but are not usually quantitatively correct. In Section 5.10.4 we consider a more sophisticated version of meanﬁeld theory for Ising models that yields more accurate values of Tc , and in Section 9.1 we consider a more general formulation of meanﬁeld theory. In Section 8.6 we discuss how to apply similar ideas to gases and liquids. In its simplest form meanﬁeld theory assumes that each spin interacts with the same eﬀective magnetic ﬁeld. The eﬀective ﬁeld is due to the external magnetic ﬁeld plus the internal ﬁeld due to all the neighboring spins. That is, spin i “feels” an eﬀective ﬁeld Heﬀ given by
q Heﬀ = J
j =1 sj + H, (5.104) where the sum over j in (5.104) is over the q nearest neighbors of i. (Recall that we have incorporated a factor of µ into H so that H in (5.104) has units of energy.) Because the orientation of the neighboring spins depends on the orientation of spin i, Heﬀ ﬂuctuates from its mean value, which is given by
q H eﬀ = J
j =1 sj + H = Jqm + H, (5.105) where sj = m. In meanﬁeld theory we ignore the deviations of Heﬀ from H eﬀ and assume that the ﬁeld at i is H eﬀ , independent of the orientation of si . This assumption is an approximation because if si is up, then its neighbors are more likely to be up. This correlation is ignored in meanﬁeld theory. The form of the mean eﬀective ﬁeld in (5.105) is the same throughout the system. The result of this approximation is that the system of N interacting spins has been reduced to a system of one spin interacting with an eﬀective ﬁeld which depends on all the other spins. The partition function for one spin in the eﬀective ﬁeld H eﬀ is Z1 =
s1 =±1 es1 H eff /kT = 2 cosh[(Jqm + H )/kT ]. (5.106) The free energy per spin is f = −kT ln Z1 = −kT ln 2 cosh[(Jqm + H )/kT ] ,
triangular Ising net,” Phys. Rev. 79, 357–364 (1950), erratum, Phys. Rev. B 7, 5017 (1973). (5.107) CHAPTER 5. MAGNETIC SYSTEMS
1.0 0.5 0.0 –0.5 –1.0 –1.5 –1.0 –0.5 0.0
(a) 257
1.0 0.5 0.0 βJq = 0.8 βJq = 2.0 stable stable –0.5 –1.0 stable unstable 0.5 m 1.0 1.5 –1.5 –1.0 –0.5 0.0
(b) 0.5 m 1.0 1.5 Figure 5.13: Graphical solution of the selfconsistent equation (5.108) for H = 0. The line y1 (m) = m represents the lefthand side of (5.108), and the function y2 (m) = tanh Jqm/kT represents the righthand side. The intersection of y1 and y2 gives a possible solution for m. The solution m = 0 exists for all T . Stable solutions m = ±m0 (m0 > 0) exist only for T suﬃciently small such that the slope Jq/kT of y2 near m = 0 is greater than 1. and the magnetization is ∂f = tanh[(Jqm + H )/kT ]. (5.108) ∂H Equation (5.108) is a selfconsistent transcendental equation whose solution yields m. The meanﬁeld that inﬂuences the mean value of m in turn depends on the mean value of m. From Figure 5.13 we see that nonzero solutions for m exist for H = 0 when qJ/kT ≥ 1. The critical temperature satisﬁes the condition that m = 0 for T ≤ Tc and m = 0 for T > Tc . Thus the critical temperature Tc is given by kTc = Jq. (5.109) m=− Problem 5.15. Numerical solutions of (5.108) Use Program IsingMeanField to ﬁnd numerical solutions of (5.108). (a) Set H = 0 and q = 4 and determine the value of the meanﬁeld approximation to the critical temperature Tc of the Ising model on a square lattice. Start with kT /Jq = 10 and then proceed to lower temperatures. Plot the temperature dependence of m. The equilibrium value of m is the solution with the lowest free energy (see Problem 5.18). (b) Determine m(T ) for the onedimensional Ising model (q = 2) and H = 0 and H = 1 and compare your values for m(T ) with the exact solution in one dimension [see (5.83)]. For T near Tc the magnetization is small, and we can expand tanh Jqm/kT (tanh x ≈ x − x3 /3 for x ≪ 1) to ﬁnd 1 (5.110) m = Jqm/kT − (Jqm/kT )3 + · · · . 3 CHAPTER 5. MAGNETIC SYSTEMS Equation (5.110) has two solutions: m(T > Tc ) = 0 and m(T < Tc ) = ± 31/2 ((Jq/kT ) − 1)1/2 . (Jq/kT )3/2 258 (5.111a) (5.111b) The solution m = 0 in (5.111a) corresponds to the high temperature paramagnetic state. The solution in (5.111b) corresponds to the low temperature ferromagnetic state (m = 0). How do we know which solution to choose? The answer can be found by calculating the free energy for both solutions and choosing the solution that gives the smaller free energy (see Problems 5.15 and 5.17). If we let Jq = kTc in (5.111b), we can write the spontaneous magnetization (the nonzero magnetization for T < Tc ) as m(T < Tc ) = 31/2 T Tc Tc − T Tc
1/2 . (5.112) We see from (5.112) that m approaches zero as a power law as T approaches Tc from below. It is convenient to express the temperature dependence of m near the critical temperature in terms of the reduced temperature ǫ = (Tc − T )/Tc [see (5.94)] and write (5.112) as m(T ) ∼ ǫβ . (5.113) From (5.112) we see that meanﬁeld theory predicts that β = 1/2. Compare this prediction to the value of β for the twodimensional Ising model (see Table 5.1). We now ﬁnd the behavior of other important physical properties near Tc . The zeroﬁeld isothermal susceptibility (per spin) is given by χ = lim For T 1 − tanh2 Jqm/kT ∂m 1 − m2 = . = 2 ∂H kT − Jq (1 − m2 ) kT − Jq (1 − tanh Jqm/kT ) 1 k (T − Tc ) (5.114) H →0 Tc we have m = 0 and χ in (5.114) reduces to χ= (T > Tc , H = 0), (5.115) where we have used the relation (5.108) with H = 0. The result (5.115) for χ is known as the CurieWeiss law. For T Tc we have from (5.112) that m2 ≈ 3(Tc − T )/Tc , 1 − m2 = (3T − 2Tc)/Tc , and χ≈ 1 1 = k [T − Tc (1 − m2 )] k [T − 3T + 2Tc ] 1 = (T Tc , H = 0). 2k (Tc − T ) (5.116a) (5.116b) We see that we can characterize the divergence of the zeroﬁeld susceptibility as the critical point is approached from either the low or high temperature side by χ ∼ ǫ−γ [see (5.98)]. The meanﬁeld prediction for the critical exponent γ is γ = 1. CHAPTER 5. MAGNETIC SYSTEMS 259 The magnetization at Tc as a function of H can be calculated by expanding (5.108) to third order in H with kT = kTc = qJ : 1 m = m + H/kTc − (m + H/kTc )3 + · · · . 3 For H/kTc ≪ m we ﬁnd m = (3H/kTc )1/3 ∝ H 1/3 (T = Tc ). (5.117) (5.118) The result (5.118) is consistent with our assumption that H/kTc ≪ m. If we use the power law dependence given in (5.99), we see that meanﬁeld theory predicts that the critical exponent δ is δ = 3, which compares poorly with the exact result for the twodimensional Ising model given by δ = 15. The easiest way to obtain the energy per spin for H = 0 in the meanﬁeld approximation is to write 1 E = − Jqm2 , (5.119) N 2 which is the average value of the interaction energy divided by 2 to account for double counting. Because m = 0 for T > Tc , the energy vanishes for all T > Tc , and thus the speciﬁc heat also vanishes. Below Tc the energy per spin is given by E 1 2 = − Jq tanh(Jqm/kT ) . N 2 Problem 5.16. Behavior of the speciﬁc heat near Tc Use (5.120) and the fact that m2 ≈ 3(Tc − T )/Tc for T Tc to show that the speciﬁc heat according to meanﬁeld theory is − C (T → Tc ) = 3k/2. (5.121)
∗ (5.120) Hence, meanﬁeld theory predicts that there is a jump (discontinuity) in the speciﬁc heat. Problem 5.17. Improved meanﬁeld theory approximation for the energy We write si and sj in terms of their deviation from the mean as si = m + ∆i and sj = m + ∆j , and write the product si sj as si sj = (m + ∆i )(m + ∆j ) = m + m(∆i + ∆j ) + ∆i ∆j .
2 (5.122a) (5.122b) We have ordered the terms in (5.122b) in powers of their deviation from the mean. If we neglect the last term, which is quadratic in the ﬂuctuations from the mean, we obtain si sj ≈ m2 + m(si − m) + m(sj − m) = −m2 + m(si + sj ). (a) Show that we can approximate the energy of interaction in the Ising model as −J si sj = + J
i,j =nn(i) (5.123) i,j =nn(i) m2 − Jm
N (si + sj )
i,j =nn(i) (5.124a) = JqN m2 si . − Jqm 2 i=1 (5.124b) CHAPTER 5. MAGNETIC SYSTEMS (b) Show that the partition function Z (T, H, N ) can be expressed as Z (T, H, N ) = e−N qJm = e−N qJm = e−N qJm
2 260 /2kT s1 =±1 ··· e(Jqm+H )
sN =±1 N P i si /kT (5.125a) (5.125b) 2 /2kT s=±1 e(qJm+H )s/kT 2 /2kT 2 cosh(qJm + H )/kT N . (5.125c) Show that the free energy per spin f (T, H ) = −(kT /N ) ln Z (T, H, N ) is given by f (T, H ) = 1 Jqm2 − kT ln 2 cosh(qJm + H )/kT . 2 (5.126) The expressions for the free energy in (5.107) and (5.126) contain both m and H rather than H only. In this case m represents a parameter. For arbitrary values of m these expressions do not give the equilibrium free energy, which is determined by minimizing f treated as a function of m. Problem 5.18. Minima of the free energy (a) To understand the meaning of the various solutions of (5.108), expand the free energy in (5.126) about m = 0 with H = 0 and show that the form of f (m) near the critical point (small m) is given by f (m) = a + b(1 − βqJ )m2 + cm4 . (5.127) Determine a, b, and c. (b) If H is nonzero but small, show that there is an additional term −mH in (5.127). (c) Show that the minimum free energy for T > Tc and H = 0 is at m = 0, and that m = ±m0 corresponds to a lower free energy for T < Tc . (d) Use Program IsingMeanField to plot f (m) as a function of m for T > Tc and H = 0. For what value of m does f (m) have a minimum? (e) Plot f (m) for T = 1 and H = 0. Where are the minima of f (m)? Do they have the same depth? If so, what is the meaning of this result? (f) Choose H = 0.5 and T = 1. Do the two minima have the same depth? The global minimum corresponds to the equilibrium or stable phase. If we quickly “ﬂip” the ﬁeld and let H → −0.5, the minimum at m ≈ 1 will become a local minimum. The system will remain in this local minimum for some time before it switches to the global minimum (see Section 5.10.6). We now compare the results of meanﬁeld theory near the critical point with the exact results for the one and twodimensional Ising models. The fact that the meanﬁeld result (5.109) for Tc depends only on q , the number of nearest neighbors, and not the spatial dimension d is one of the inadequacies of the simple version of meanﬁeld theory that we have discussed. The simple meanﬁeld theory even predicts a phase transition in one dimension, which we know is incorrect. In Table 5.2 the meanﬁeld predictions for Tc are compared to the bestknown estimates of the CHAPTER 5. MAGNETIC SYSTEMS lattice square hexagonal diamond simple cubic bcc fcc d 2 2 3 3 3 3 q 4 6 4 6 8 12 Tmf /Tc 1.763 1.648 1.479 1.330 1.260 1.225 261 Table 5.2: Comparison of the meanﬁeld predictions for the critical temperature of the Ising model with exact results and the bestknown estimates for diﬀerent spatial dimensions d and lattice symmetries. critical temperatures for the Ising model for two and threedimensional lattices. We see that the meanﬁeld theory predictions for Tc improve as the number of neighbors increases for a given dimension. The inaccuracy of meanﬁeld theory is due to the fact that it ignores correlations between the spins. In Section 5.10.4 we discuss more sophisticated treatments of meanﬁeld theory that include shortrange correlations between the spins and yield better estimates of the critical temperature, but not the critical exponents. Meanﬁeld theory predicts that various thermodynamic properties exhibit power law behavior near Tc , in qualitative agreement with the exact solution of the twodimensional Ising model and the known properties of the threedimensional Ising model. This qualitative behavior is characterized by critical exponents. The meanﬁeld predictions for the critical exponents β , γ , and δ are β = 1/2, γ = 1, and δ = 3 respectively (see Table 5.1). (The meanﬁeld theory predictions of the other critical exponents are given in Section 9.1.) These values of the critical exponents do not agree with the results of the Onsager solution of the twodimensional Ising model (see Table 5.1). Also meanﬁeld theory predicts a jump in the speciﬁc heat, whereas the Onsager solution predicts a logarithmic divergence. Similar disagreements between the predictions of meanﬁeld theory and the known critical exponents are found in three dimensions, but the discrepancies are not as large. We also note that the meanﬁeld results for the critical exponents are independent of dimension. In Section 9.1 we discuss a more general version of meanﬁeld theory, which is applicable to a wide variety of systems, and shows why all meanﬁeld theories predict the same values for the critical exponents independent of dimension. Problem 5.19. Improvement of meanﬁeld theory with dimension From Table 5.1 we see that the predictions of meanﬁeld theory increase in accuracy with increasing dimension. Why is this trend reasonable? Why meanﬁeld theory fails . The main assumption of meanﬁeld theory is that each spin feels the same eﬀective magnetic ﬁeld due to all the other spins. That is, meanﬁeld theory ignores the ﬂuctuations in the eﬀective ﬁeld. But if meanﬁeld theory ignores ﬂuctuations, why does the susceptibility diverge near the critical point? (Recall that the susceptibility is a measure of the ﬂuctuations of the magnetization.) Because the ﬂuctuations are ignored in one context, but not another, we see that meanﬁeld theory carries with it the seeds of its own destruction. That is, meanﬁeld theory does not treat the ﬂuctuations consistently. This inconsistency is unimportant if the ﬂuctuations are weak. CHAPTER 5. MAGNETIC SYSTEMS 262 Because the ﬂuctuations become more important as the system approaches the critical point, we expect that meanﬁeld theory breaks down for T close to Tc . A criterion for the range of temperatures for which meanﬁeld theory theory is applicable is discussed in Section 9.1, page 434, where it is shown that the ﬂuctuations can be ignored if
− ξ0 d ǫ(d/2)−2 ≪ 1, (5.128) where ξ0 is the correlation length at T = 0 and is proportional to the eﬀective range of interaction. The inequality in (5.128) is always satisﬁed for d > 4 near the critical point where ǫ ≪ 1. That is, meanﬁeld theory yields the correct results for the critical exponents in higher than four dimensions. (In four dimensions the power law behavior is modiﬁed by logarithmic factors.) In a conventional superconductor such as tin, ξ0 ≈ 2300 ˚, and the meanﬁeld theory of a superconductor (known A as BCS theory) is applicable near the superconducting transition for ǫ as small as 10−14 . 5.7.1 *Phase diagram of the Ising model Nature exhibits two qualitatively diﬀerent kinds of phase transitions. The more familiar kind, which we observe when ice freezes or water boils, involves a discontinuous change in various thermodynamic quantities such as the energy and the entropy. For example, the density as well as the energy and the entropy change discontinuously when water boils and when ice freezes. This type of phase transition is called a discontinuous or ﬁrstorder transition. We will discuss ﬁrstorder transitions in the context of gases and liquids in Chapter 7. The other type of phase transition is more subtle. In this case thermodynamic quantities such as the energy and the entropy are continuous, but various derivatives such as the speciﬁc heat and the compressibility of a ﬂuid and the susceptibility of a magnetic system show divergent behavior at the phase transition. Such transitions are called continuous phase transitions. We have seen that the Ising model in two dimensions has a continuous phase transition in zero magnetic ﬁeld such that below the critical temperature Tc there is a nonzero spontaneous magnetization, and above Tc the mean magnetization vanishes as shown by the solid curve in the phase diagram in Figure 5.14. The threedimensional Ising model has the same qualitative behavior, but the values of the critical temperature and the critical exponents are diﬀerent. The behavior of the Ising model is qualitatively diﬀerent if we apply an external magnetic ﬁeld H . If H = 0, the magnetization m is nonzero at all temperatures and has the same sign as H (see Figure 5.15). The same information is shown in a diﬀerent way in Figure 5.14. Each point in the unshaded region corresponds to an equilibrium value of m for a particular value of T and H .10 At a phase transition at least one thermodynamic quantity diverges or has a discontinuity. For example, both the speciﬁc heat and the susceptibility diverge at Tc for a ferromagnetic phase
10 In a ferromagnetic material such as iron, nickel, and cobalt, the net magnetization frequently vanishes even below Tc . In these materials there are several magnetic domains within which the magnetization is nonzero. These domains are usually oriented at random, leading to zero net magnetization for the entire sample. Such a state is located in the shaded region of Figure 5.14. The randomization of the orientation of the domains occurs when the metal is formed and cooled below Tc and is facilitated by crystal defects. When a piece of iron or similar material is subject to an external magnetic ﬁeld, the domains align, and the iron becomes “magnetized.” When the ﬁeld is removed the iron remains magnetized. If the iron is subject to external forces such as being banged with a hammer, the domains can be randomized again, and the iron loses its net magnetization. The Ising model is an example of a singledomain ferromagnet. CHAPTER 5. MAGNETIC SYSTEMS 263 1.0 m
0.5 H>0 0.0 –0.5 H<0 –1.0 0.0 0.5 1.0 T/Tc 1.5 2.0 Figure 5.14: Sketch of the phase diagram for an Ising ferromagnet. The bold line represents the magnetization for H = 0. Equilibrium magnetization values for H = 0 are possible at any point in the unshaded region. Points in the shaded region represent nonequilibrium values of m. transition with H = 0. For H = 0 and T > Tc there is no phase transition as the magnetic ﬁeld is decreased to zero and then to negative values because no quantity diverges or has a discontinuity. In contrast, for T < Tc there is a transition because as we change the ﬁeld from H = 0+ to H = 0− m changes discontinuously from a positive value to a negative value (see Figure 5.15). Problem 5.20. Ising model in an external magnetic ﬁeld Use Program Ising2D to simulate the Ising model on a square lattice at a given temperature. Choose N = L2 = 322 . Run for at least 200 Monte Carlo steps per spin at each value of the ﬁeld. (a) Set H = 0.2 and T = 3 and estimate the approximate value of the magnetization. Then change the ﬁeld to H = 0.1 and continue updating the spins (do not press New) so that the simulation is continued from the last microstate. Note the value of the magnetization. Continue this procedure with H = 0, then H = −0.1 and then H = −0.2. Do your values of the magnetization change abruptly as you change the ﬁeld? Is there any indication of a phase transition as you change H ? (b) Repeat the same procedure as in part (a) at T = 1.8 which is below the critical temperature. What happens now? Is there evidence of a sudden change in the magnetization as the direction of the ﬁeld is changed? (c) Use Program Ising2DHysteresis with T = 1.8, the initial magnetic ﬁeld H = 1, ∆H = 0.01, and 10 mcs for each value of H . The program plots the mean magnetization for each value of H , and changes H by ∆H until H reaches H = −1, when it changes ∆H to −∆H . Describe what you obtain and why it occurred. The resulting curve is called a hysteresis loop, and is CHAPTER 5. MAGNETIC SYSTEMS
1.0 264 m
0.5 T < Tc T = Tc T > Tc 0.0 –0.5 –1.0 –2.0 –1.0 0.0 1.0 H 2.0 Figure 5.15: The equilibrium values of m as a function of the magnetic ﬁeld H for T > Tc , T = Tc , and T < Tc . The plot for T > Tc is smooth in contrast to the plot for T < Tc which has a discontinuity at H = 0. For T = Tc there is no discontinuity, and the function m(H ) has an inﬁnite slope at H = 0. characteristic of discontinuous phase transitions. An example of the data you can obtain in this way is shown in Figure 5.16. (d) Change the number of mcs per ﬁeld value to 1 and view the resulting plot for m versus H . Repeat for mcs per ﬁeld value equal to 100. Explain the diﬀerences you see. 5.8 *Simulation of the Density of States The probability that a system in equilibrium with a heat bath at a temperature T has energy E is given by Ω(E ) −βE P (E ) = e , (5.129) Z where Ω(E ) is the number of states with energy E ,11 and the partition function Z = E Ω(E )e−βE . If Ω(E ) is known, we can calculate the mean energy (and other thermodynamic quantities) at any temperature from the relation 1 E Ω(E )e−βE . (5.130) E= Z
E Hence, the quantity Ω(E ) is of much interest. In the following we discuss an algorithm that directly computes Ω(E ) for the Ising model. In this case the energy is a discrete variable and hence the quantity we wish to compute is the number of spin microstates with the same energy.
11 The quantity Ω(E ) is the number of states with energy E for a system such as the Ising model which has discrete values of the energy. It is common to refer to Ω(E ) as the density of states even when the values of E are discrete. CHAPTER 5. MAGNETIC SYSTEMS 265 1.0 m
0.5 0.0 –0.5 –1.0 –1.0 –0.5 0.0 0.5 H 1.0 Figure 5.16: Hysteresis curve obtained from the simulation of the twodimensional Ising model at T = 1.8. The ﬁeld was reduced by ∆H = 0.01 every 10 mcs. The arrows indicate the direction in which the magnetization is changing. Note that starting from saturation at m = 1, a “coercive ﬁeld” of about H = −0.25 is needed to reduce the magnetization to zero. The magnetization at zero magnetic ﬁeld is called the “remnant” magnetization. The equilibrium value of m drops discontinuously from a positive value to a negative value as H decreases from H = 0+ to H = 0− . Hysteresis is a nonequilibrium phenomenon and is a dynamic manifestation of a system that remains in a local minimum for some time before it switches to the global minimum. Suppose that we were to try to compute Ω(E ) by doing a random walk in energy space by ﬂipping the spins at random and accepting all microstates that we obtain in this way. The histogram of the energy, H (E ), the number of visits to each possible energy E of the system, would become proportional to Ω(E ) if the walk visited all possible microstates many times. In practice, it would be impossible to realize such a long random walk given the extremely large number of microstates. For example, an Ising model with N = 100 spins has 2100 ≈ 1.3 × 1030 microstates. An even more important limitation on doing a simple random walk to determine Ω(E ) is that the walk would spend most of its time visiting the same energy values over and over again and would not reach the values of E that are less probable. The idea of the WangLandau algorithm is to do a random walk in energy space by ﬂipping single spins at random and accepting the changes with a probability that is proportional to the reciprocal of the density of states. In this way energy values that would be visited often using a simple random walk would be visited less often because they have a larger density of states. There is only one problem – we don’t know Ω(E ). We will see that the WangLandau12 algorithm estimates Ω(E ) at the same time that it does a random walk. To implement the algorithm we begin with an arbitrary microstate and a guess for the density of states. The simplest guess is to set Ω(E ) = 1 for all possible energies E . The algorithm can be
12 We have mentioned the name Landau several times in the text. This Landau is not Lev D. Landau, but David Landau, a wellknown physicist at the University of Georgia. CHAPTER 5. MAGNETIC SYSTEMS summarized by the follow steps. 266 1. Choose a spin at random and make a trial ﬂip. Compute the energy before, E1 , and after the ﬂip, E2 , and accept the change with probability ˜ ˜ p(E1 → E2 ) = min(Ω(E1 )/Ω(E2 ), 1), (5.131) ˜ ˜ where Ω(E ) is the current estimate of Ω(E ). Equation (5.131) implies that if Ω(E2 ) ≤ ˜ E1 ), the state with energy E2 is always accepted; otherwise, it is accepted with probability Ω( ˜ ˜ Ω(E1 )/Ω(E2 ). That is, the state with energy E2 is accepted if a random number r satisﬁes ˜ ˜ the condition r ≤ Ω(E1 )/Ω(E2 ). After the trial ﬂip the energy of the system is E2 if the change is accepted or remains at E1 if the change is not accepted. ˜ 2. To estimate Ω(E ) multiply the current value of Ω(E ) by a modiﬁcation factor f > 1: ˜ ˜ Ω(E ) = f Ω(E ). We also update the existing entry for H (E ) in the energy histogram H (E ) → H (E ) + 1. (5.133) (5.132) ˜ Because Ω(E ) becomes very large, we must work with the logarithm of the density of states, ˜ so that ln Ω(E ) will ﬁt into double precision numbers. Therefore, each update of the density ˜ ˜ of states is implemented as ln Ω(E ) → ln Ω(E ) + ln f , and the ratio of the density of states is ˜ E1 ) − ln Ω(E2 )]. A reasonable choice of the initial modiﬁcation factor ˜ computed as exp[ln Ω( is f = f0 = e ≈ 2.71828 . . .. If f0 is too small, the random walk will need a very long time to reach all possible energies. Too large a choice of f0 will lead to large statistical errors. 3. Proceed with the random walk in energy space until an approximately ﬂat histogram H (E ) is obtained, that is, until all the possible energy values are visited an approximately equal number of times. Because it is impossible to obtain a perfectly ﬂat histogram, we will say that H (E ) is “ﬂat” when H (E ) for all possible E is not less than ∆ of the average histogram H (E ); ∆ is chosen according to the size and the complexity of the system and the desired accuracy of the density of states. For the twodimensional Ising model on small lattices, ∆ can be chosen to be as high as 0.95, but for large systems the criterion for ﬂatness might never be satisﬁed if ∆ is too close to 1. 4. Once the ﬂatness criterion has been satisﬁed, reduce the modiﬁcation factor f using a function √ such as f1 = f0 , reset the histogram to H (E ) = 0 for all values of E , and begin the next iteration of the random walk during which the density of states is modiﬁed by f1 at each trial ﬂip. The density of states is not reset during the simulation. We continue performing the random walk until the histogram H (E ) is again ﬂat. We then reduce the modiﬁcation √ factor, fi+1 = fi , reset the histogram to H (E ) = 0 for all values of E , and continue the random walk. 5. The simulation is stopped when f is smaller than a predeﬁned value [such as f = exp(10−8 ) ≈ 1.00000001]. The modiﬁcation factor acts as a control parameter for the accuracy of the density of states during the simulation and also determines how many Monte Carlo sweeps are necessary for the entire simulation. CHAPTER 5. MAGNETIC SYSTEMS 267 The algorithm provides an estimate of the density of states because if the current estimate ˜ of Ω(E ) is too low, then the acceptance criterion (5.131) pushes the system to states with lower ˜ ˜ ˜ Ω(E ), thus increasing Ω(E ). If Ω(E ) is too high, the reverse happens. Gradually, the calculation tends to the true value of Ω(E ). At the end of the simulation, the algorithm provides only a relative density of states. To determine the normalized density of states Ω(E ), we can either use the fact that the total number of states for the Ising model is Ω(E ) = 2N , (5.134)
E or that the number of ground states (for which E = −2N J ) is 2. The latter normalization ensures the accuracy of the density of states at low energies which is important in the calculation of thermodynamic quantities at low temperatures. If we apply (5.134), we cannot guarantee the accuracy of Ω(E ) for energies at or near the ground state, because the rescaling factor is dominated by the maximum density of states. We may use one of these normalization conditions to obtain the absolute density of states and the other normalization condition to check the accuracy of our result.
∗ Problem 5.21. WangLandau algorithm for the Ising model Program IsingDensityOfStates implements the WangLandau algorithm for the Ising model on a square lattice. (a) Calculate the exact values of Ω(E ) for the 2 × 2 Ising model. Run the simulation for L = 2 and verify that the computed density of states is close to your exact answer. (b) Choose larger values of L, for example, L = 16, and describe the qualitative energy dependence of Ω(E ). (c) The program also computes the speciﬁc heat as a function of temperature using the estimated ˜ value of Ω(E ). Describe the qualitative temperature dependence of the speciﬁc heat. The Potts model . The Potts model is a generalization of the Ising model in which each lattice site contains an entity (a spin) that can be in one of q states. If two nearest neighbor sites are in the same state, then the interaction energy is −K . The interaction energy is zero if they are in diﬀerent states. Potts models are useful for describing the absorption of molecules on crystalline surfaces and the behavior of foams, for example, and exhibit a discontinuous or continuous phase transition depending on the value of q . The Potts model exhibits a phase transition between a high temperature phase where the q states equally populate the sites and a low temperature phase where one of the q states is more common than the others. In two dimensions the transition between these two phases is ﬁrstorder (see Section 5.7.1) for q > 4 and is continuous otherwise. In Problem 5.22 we explore how the WangLandau algorithm can provide some insight into the nature of the Potts model and its phase transitions.
∗ Problem 5.22. Application of WangLandau algorithm to the Potts model CHAPTER 5. MAGNETIC SYSTEMS 268 (a) (b) Figure 5.17: (a) A typical microstate of the Ising model. (b) The same microstate in the lattice gas picture with spin up replaced by a particle and spin down replaced by an empty cell. (a) What is the relation of the q = 2 Potts model to the usual Ising model? In particular, what is the relation of the interaction energy K deﬁned in the Potts model and the interaction energy J deﬁned in the Ising model? (b) Program PottsModel implements the WangLandau algorithm for the Potts model on a square lattice. Run the program for L = 16 and various values of q and verify that the peak in the heat capacity occurs near the known exact value for the transition temperature given by √ Tc = (ln (1 + q ))−1 . You will need more than 100,000 Monte Carlo steps to obtain reliable data. (c) Choose q = 2 and q = 3 and observe the energy distribution P (E ) = Ωe−βE at T = Tc . Do you see one peak or two? (d) Choose q = 10 and observe the energy distribution P (E ) at T = Tc . You should notice two peaks in this distribution. Discuss why the occurrence of two peaks is appropriate for a ﬁrstorder transition. 5.9 *Lattice Gas The Ising model is useful not only because it is the simplest model of magnetism, but also because it can be applied to many other systems. Two common applications are to ﬂuids and to binary alloys. In the ﬂuid model a spin of +1 represents a particle and a spin of −1 represents a void (see Figure 5.17). The hard core repulsion between particles at short distances is modeled by the restriction that there is at most one particle per site. The shortrange attractive interaction between particles is modeled by the nearestneighbor Ising interaction. Binary alloys are modeled in a similar way with +1 representing one type of atom and −1 representing a second type. The Ising interaction models the tendency of like atoms to be near each other because the attraction between like atoms is stronger than that between unlike atoms in binary alloys. In this section we will focus on the ﬂuid model, which is called the lattice gas. CHAPTER 5. MAGNETIC SYSTEMS 269 We could proceed by taking the form of the Ising energy given in (5.35) and converting all our previous results in the language of magnetism to the language of ﬂuids. For example, when most spins are up the system is mostly a liquid, and when most spins are down the system is a gas. However, it is useful to change variables so that we can more directly describe the behavior of the system in the language of particles. For this purpose we deﬁne a new variable ni ≡ (si + 1)/2, so that ni = 0 or 1. The energy function (5.35) becomes
N N E = −J i,j =nn(i) (2ni − 1)(2nj − 1) − H i=1 (2ni − 1). (5.135) We expand and rearrange terms and ﬁnd that on a square lattice
N N E = −4J i,j =nn(i) ni nj − (2H − 8J ) i=1 ni + N (H − 2J ).
N i (5.136) To obtain the factor of 8 in (5.136) note that N i,j =nn(i) 2 ni = 4 2ni , where the factor of 4
N arises from the sum over j and is the number of nearest neighbors of i. The sum i,j =nn(i) 2nj gives the same contribution. To avoid double counting we need to divide the sums by a factor of 2. We deﬁne the energy u0 ≡ 4J and the chemical potential µ ≡ 2H − 8J , and express the energy of the lattice gas as
N N E = −u0 i,j =nn(i) ni nj − µ i=1 ni + N (H − 2J ). (5.137) The constant term N (H − 2J ) can be eliminated by redeﬁning the zero of energy. It is natural to ﬁx the temperature T and external magnetic ﬁeld H of the Ising model because we can control these quantities experimentally and in simulations. Hence, we usually simulate and calculate the properties of the Ising model in the canonical ensemble. The Metropolis algorithm for simulating the Ising model ﬂips individual spins, which causes the magnetization to ﬂuctuate. Because the magnetization is not conserved, the number of particles in the lattice gas context is not conserved, and hence the same Metropolis algorithm is equivalent to the grand canonical ensemble for a lattice gas. We can modify the Metropolis algorithm to simulate a lattice gas in the canonical ensemble with the number of particles ﬁxed. Instead of ﬂipping individual spins (single spin ﬂip dynamics), we have to interchange two spins. The algorithm proceeds by choosing a pair of nearest neighbor spins at random. If the two spins are parallel, we include the unchanged microstate in various averages. If the two spins are antiparallel, we interchange the two spins and compute the trial change in the energy ∆E as before and accept the trial change with the usual Boltzmann probability. (This algorithm is called spin exchange or Kawasaki dynamics.) Although the Ising model and the lattice gas are equivalent and all the critical exponents are the same, the interpretation of the phase diagram diﬀers. Suppose that the number of occupied sites equals the number of unoccupied sites. In this case the transition from high to low temperature is continuous. For T > Tc the particles exist in small droplets and the voids exist in small bubbles. In this case the system is neither a liquid or a gas. Below Tc the particles coalesce into a macroscopically large cluster, and the bubbles coalesce into a large region of unoccupied sites. CHAPTER 5. MAGNETIC SYSTEMS 270 This change is an example of phase separation, and the simultaneous existence of both a gas and liquid is referred to as gasliquid coexistence.13 The order parameter of the lattice gas is taken to be ρ∗ ≡ (ρL − ρG )/ρL , where ρL is the particle density of the liquid region and ρG is the density of the gas region. Above Tc there is no phase separation (there are no separate liquid and gas regions) and thus ρ∗ = 0. At T = 0, ρL = 1 and ρG = 0.14 The power law behavior of ρ∗ as T → Tc from below Tc is described by the critical exponent β , which has the same value as that for the magnetization in the Ising model. The equality of the critical exponents for the Ising and lattice gas models as well as more realistic models of liquids and gases is an example of universality of critical phenomena for qualitatively diﬀerent systems (see Sections 9.5 and 9.6). If the number of occupied and unoccupied sites is unequal, then the transition from a ﬂuid to twophase coexistence as the temperature is lowered is discontinuous. For particle systems it is easier to analyze the transition by varying the pressure rather than the temperature. As we will discuss in Chapter 7, there is a jump in the density as the pressure is changed along an isotherm in a real ﬂuid. This density change can occur either for a ﬁxed number of particles with a change in the volume or as a change in the number of particles for a ﬁxed volume. We will consider the latter by discussing the lattice gas in the grand canonical ensemble. From (2.168) we know that the thermodynamic potential Ω associated with the grand canonical ensemble is given by Ω = −P V . For the lattice gas the volume V is equal to the number of sites N , a ﬁxed quantity. We know from (4.144) that Ω = −kT ln ZG , where ZG is the grand partition function. The grand partition function for the lattice gas with the energy given by (5.137) is identical to the partition function for the Ising model in a magnetic ﬁeld with the energy given by (5.35), because the diﬀerence is only a change of variables. The free energy for the Ising model in a magnetic ﬁeld is F = −kT ln Z . Because Z = ZG , we have that F = Ω, and thus we conclude that −P V for the lattice gas equals F for the Ising model. This identiﬁcation will allow us to understand what happens to the density as we change the pressure along an isotherm. In a lattice gas the density is the number of particles divided by the number of sites or ρ = ni . In Ising language ρ is ρ= 1 (si + 1) = (m + 1)/2, 2 (5.138) where m is the magnetization. Because −P V = F , changing the pressure at ﬁxed temperature and volume in the lattice gas corresponds to changing the free energy F (T, V, H ) by changing the ﬁeld H in the Ising model. We know that when H changes from 0+ to 0− for T < Tc there is a jump in the magnetization from a positive to a negative value. From (5.138) we see that this jump corresponds to a jump in the density in the lattice gas model, corresponding to a change in the density from a liquid to a gas. Problem 5.23. Simulation of the twodimensional lattice gas (a) What is the value of the critical temperature Tc for the lattice gas in two dimensions (in units of u0 /k )? (b) Program LatticeGas simulates the lattice gas on a square lattice of linear dimension L. The initial state has all the particles at the bottom of the simulation cell. Choose L = 32 and set
13 If the system were subject to gravity, the liquid region would be at the bottom of a container and the gas would be at the top. 14 For T > 0 the cluster representing the liquid region would have some unoccupied sites and thus ρ < 1, and L the gas region would have some particles so that ρG > 0. CHAPTER 5. MAGNETIC SYSTEMS 271 the gravitational ﬁeld equal to zero. Do a simulation at T = 0.4 with N = 600 particles. After a few Monte Carlo steps you should see the bottom region of particles (green sites) develop a few small holes or bubbles and the unoccupied region contain a few isolated particles or small clusters of particles. This system represents a liquid (the predominately green region) in equilibrium with its vapor (the mostly white region). Record the energy. To speed up the simulation set steps per display equal to 100. (c) Increase the temperature in steps of 0.05 until T = 0.7. At each temperature run for at least 10,000 mcs to reach equilibrium and then press the Zero Averages button. Run for at least 20,000 mcs before recording your estimate of the energy. Describe the visual appearance of the positions of the particle and empty sites at each temperature. At what temperature does the one large liquid region break up into many pieces, such that there is no longer a sharp distinction between the liquid and vapor region? At this temperature there is a single ﬂuid phase. Is there any evidence from your estimates of the energy that a transition from a twophase to a onephase system has occurred? Repeat your simulations with N = 200. (d) Repeat part (c) with N = 512. In this case the system will pass through a critical point. The change from a one phase to a twophase system occurs continuously in the thermodynamic limit. Can you detect this change or does the system look similar to the case in part (c)? (e) If we include a gravitational ﬁeld, the program removes the periodic boundary conditions in the vertical direction, and thus sites in the top and bottom rows have three neighbors instead of four. The ﬁeld should help deﬁne the liquid and gas regions. Choose g = 0.01 and repeat the above simulations. Describe the diﬀerences you see. (f) Simulate a lattice gas of N = 2048 particles on a L = 64 lattice at T = 2.0 with no gravitational ﬁeld for 5000 mcs. Then change the temperature to T = 0.2. This process is called a (temperature) quench, and the resulting behavior is called spinodal decomposition. The domains grow very slowly as a function of time. Discuss why it is diﬃcult for the system to reach its equilibrium state for which there is one domain of mostly occupied sites in equilibrium with one domain of mostly unoccupied sites. 5.10
5.10.1 Supplementary Notes
The Heisenberg model of magnetism Classical electromagnetic theory tells us that magnetic ﬁelds are due to electrical currents and changing electric ﬁelds, and that the magnetic ﬁelds far from the currents are described by a magnetic dipole. It is natural to assume that magnetic eﬀects in matter are due to microscopic current loops created by the motion of electrons in atoms. However, it was shown by Niels Bohr in his doctoral thesis of 1911 and independently by Johanna H. van Leeuwen in her 1919 doctoral thesis that diamagnetism does not exist in classical physics (see Mattis). Magnetism is a quantum phenomenon. In the context of magnetism the most obvious new physics due to quantum mechanics is the existence of an intrinsic magnetic moment. The intrinsic magnetic moment is proportional to the intrinsic spin, another quantum mechanical property. We will now derive an approximation CHAPTER 5. MAGNETIC SYSTEMS 272 for the interaction energy between two magnetic moments. Because the electrons responsible for magnetic behavior are localized near the atoms of a regular lattice in most magnetic materials, we consider the simple case of two localized electrons. Each electron has spin 1/2 and are aligned up or down along the axis speciﬁed by the applied magnetic ﬁeld. The electrons interact with each other and with nearby atoms and are described in part by the spatial wavefunction ψ (r1 , r2 ). This wavefunction must be multiplied by the spin eigenstates to obtain the actual state of the two electron system. We denote the basis for the spin eigenstates as ↑↑ , ↓↓ , ↑↓ , ↓↑ , (5.139) where the arrows correspond to the spin of the electrons. These states are eigenstates of the ˆ ˆ z component of the total spin angular momentum15 Sz such that Sz operating on any of the states in (5.139) has an eigenvalue equal to the sum of the spins in the z direction. For example, ˆ ˆ ˆ ˆ Sz ↑↑ = 1↑↑ and Sz ↑↓ = 0↑↓ . Similarly, Sx or Sy give zero if either operator acts on these eigenstates. Because electrons are fermions, the basis states in (5.139) are not acceptable because if two electrons are interchanged, the wave function must be antisymmetric. Thus, ψ (r1 , r2 ) = +ψ (r2 , r1 ) if the spin state is antisymmetric, and ψ (r1 , r2 ) = −ψ (r2 , r1 ) if the spin state is symmetric. The simplest normalized linear combinations of the eigenstates in (5.139) that satisfy this condition are 1 √ [↑↓ − ↓↑ ], 2 ↑↑ , 1 √ [↑↓ + ↓↑ ], 2 ↓↓ . (5.140a) (5.140b) (5.140c) (5.140d) The state in (5.140a) is antisymmetric, because interchanging the two electrons leads to minus the original state. This state has a total spin S = 0 and is called the singlet state. The collection of the last three states is called the triplet state and has S = 1. If the spins are in the triplet state, then ψ (r1 , r2 ) = −ψ (r2 , r1 ). Similarly, if the spins are in the singlet state, then ψ (r1 , r2 ) = +ψ (r2 , r1 ). Hence, when r1 = r2 , ψ is zero for the triplet state, and thus the electrons stay further apart, and their electrostatic energy is smaller. For the singlet state at r1 = r2 , ψ is nonzero, and thus the electrons can be closer to each other, with a larger electrostatic energy. To ﬁnd a relation between the energy and the spin operators we note that ˆˆ ˆˆ ˆ ˆ ˆ2 ˆ2 S · S = (S1 + S2 )2 = S1 + S2 + 2 S1 · S2 , (5.141) ˆ ˆ2 where the operator S is the total spin. Because both electrons have spin 1/2, the eigenvalues of S1 2 ˆ ˆ and S2 are equal and are given by (1/2)(1 + 1/2) = 3/4. We see that the eigenvalue S of S is zero ˆ for the singlet state and is one for the triplet state. Hence, the eigenvalue of S 2 is S (S + 1) = 0 ˆˆ for the singlet state and (1(1 + 1) = 2 for the triplet state. Similarly, the eigenvalue S12 of S1 · S2 equals −3/4 for the singlet state and 1/4 for the triplet state. These considerations allows us to write E = c − JS12 , (5.142)
15 We will denote operators with a caret symbol in this section. CHAPTER 5. MAGNETIC SYSTEMS 273 where c is a constant and J is known as the exchange constant. If we denote Etriplet and Esinglet as the triplet energy and the singlet energy, respectively, and let J = Esinglet − Etriplet , we can determine c and ﬁnd 1 (5.143) E = (Esinglet + 3Etriplet) − JS12 . 4 You can check (5.143) by showing that when S12 = −3/4, E = Esinglet and when S12 = 1/4, E = Etriplet . The term in parentheses in (5.143) is a constant and can be omitted by suitably deﬁning the zero of energy. The second term represents a convenient form of the interaction between two spins. The total energy of the Heisenberg model of magnetism is based on the form (5.143) for the spinspin interaction and is expressed as
N N ˆ H=− i<j =1 ˆˆ Jij Si · Sj − gµ0 H · ˆ Si ,
i=1 (5.144) where gµ0 is the magnetic moment of the electron. Usually we combine the factors of g and µ0 into H and write the Heisenberg Hamiltonian as
N N ˆ H=− i<j =1 ˆˆ Jij Si · Sj − H ˆ Sz,i
i=1 (Heisenberg model). (5.145) The form (5.145) of the interaction energy is known as the Heisenberg model. The exchange ˆ ˆ constant Jij can be either positive or negative. Note that S as well as the Hamiltonian H is an operator, and that the Heisenberg model is quantum mechanical in nature. The distinction ˆ between the operator H and the magnetic ﬁeld H will be clear from the context. The Heisenberg model assumes that we can treat all interactions in terms of pairs of spins. This assumption means that the magnetic ions in the crystal must be suﬃciently far apart that the overlap of their wavefunctions is small. We have also neglected any orbital contribution to the total angular momentum. In addition, dipolar interactions can be important and lead to a coupling between the spin degrees of freedom and the relative displacements of the magnetic ions. It is very diﬃcult to obtain the exact Hamiltonian from ﬁrst principles. The Heisenberg model is the starting point for most microscopic models of magnetism. We can go to the classical limit S → ∞, consider spins with one, two, or three components, place the spins on lattices of any dimension and any crystal structure, and take J to be positive, negative, random, nearestneighbor, longrange, and so on. In addition, we can include other interactions such as the interaction of an electron with an ion. The theoretical possibilities are very rich as are the types of magnetic materials of interest experimentally. 5.10.2 Low temperature expansion The existence of exact analytical solutions for systems with nontrivial interactions is the exception, and we usually must be satisﬁed with approximate solutions with limited ranges of applicability. If the ground state is known and if we can determine the excitations from the ground state, we can determine the behavior of a system at low temperatures. CHAPTER 5. MAGNETIC SYSTEMS 274 (a) (b) Figure 5.18: (a) The ground state of N = 5 Ising spins in an external magnetic ﬁeld H . For toroidal boundary conditions, the ground state energy is E0 = −5J − 5H . (b) The ﬂip of a single spin of N = 5 Ising spins. The corresponding energy cost is 4J + 2H . (a) (b) Figure 5.19: Microstates corresponding to two ﬂipped spins of a system of N = 5 spins in one dimension. In (a) the ﬂipped spins are nearest neighbors and in (b) the ﬂipped spins are not nearest neighbors. To understand the nature of this class of approximations we consider the onedimensional Ising model at low temperatures. We know that the ground state corresponds to all spins completely aligned. When we raise the temperature slightly above T = 0, the system can raise its energy by ﬂipping one or more spins. At a given temperature we can consider the excited states corresponding to 1, 2, . . . , f ﬂipped spins. These f spins may be connected or may consist of disconnected groups. Each number of ﬂipped spins corresponds to a term in the low temperature expansion of the partition function. As an example, consider a system of N = 5 spins with toroidal boundary conditions. The ground state is shown in Figure 5.18(a). The energy cost of ﬂipping a single spin is 4J + 2H . A typical microstate with one ﬂipped spin is shown in Figure 5.18(b). (The energy of interaction of the ﬂipped spin with its two neighbors changes from −2J to +2J .) Because the ﬂipped spin can be at N = 5 diﬀerent sites, we have Z = [1 + 5 e−β (4J +2H ) ]e−βE0 where E0 = −5(J + H ). (f = 1), (5.146) The next higher energy excitation consists of a pair of ﬂipped spins with one type of contribution arising from pairs that are nearest neighbors [see Figure 5.19(a)] and the other type arising from pairs that are not nearest neighbors [see Figure 5.19(b)]. We leave it as an exercise (see Problem 5.24) to determine the corresponding energies and the number of diﬀerent ways that this type of excitation occurs.
∗ Problem 5.24. Low temperature expansion for ﬁve spins (a) Determine the contribution to the partition function corresponding to two ﬂipped spins out of N = 5. (Use toroidal boundary conditions.) (b) Enumerate the 25 microstates of the N = 5 Ising model in one dimension and classify the microstates corresponding to the energy of the microstate. Then use your results to ﬁnd the CHAPTER 5. MAGNETIC SYSTEMS low temperature expansion of Z . Express your result for Z in terms of the variables u = e−2J/kT and w = e−2H/kT .
∗ 275 (5.147a) (5.147b) Problem 5.25. Low temperature behavior of N spins (a) Generalize your results to N spins in one dimension and calculate Z corresponding to f = 1 and f = 2 ﬂipped spins. Find the free energy, mean energy, heat capacity, magnetization, and susceptibility, assuming that u ≪ 1 and w ≪ 1 [see (5.147)] so that you can use the approximation ln (1 + ǫ) ≈ ǫ. (b) Generalize the low temperature expansion for the onedimensional Ising model by expanding in the number of domain walls, and show that the expansion can be summed exactly (see Problem 5.34). (The low temperature expansion of the Ising model can be summed only approximately in higher dimensions using what are known as Pad´ approximants.) e 5.10.3 High temperature expansion At high temperatures for which J/kT ≪ 1, the eﬀects of the interactions between the spins become small, and we can develop a perturbation method based on expanding Z in terms of the small parameter J/kT . For simplicity, we consider the Ising model in zero magnetic ﬁeld. We write (5.148) eβJsi sj , ZN =
{si =±1} i,j =nn(i) where the sum is over all states of the N spins, and the product is restricted to nearest neighbor pairs of sites ij in the lattice. We ﬁrst use the identity eβJsi sj = cosh βJ + si sj sinh βJ = (1 + vsi sj ) cosh βJ, where v = tanh βJ. (5.150) The identity (5.149) can be demonstrated by considering the various cases si , sj = ±1 (see Problem 5.43). The variable v approaches zero as T → ∞ and will be used as an expansion parameter instead of J/kT for reasons that will become clear. Equation (5.148) can now be written as ZN = (cosh βJ )p
{si } ij (5.149) (1 + vsi sj ), (5.151) where p is the total number of nearest neighbor pairs in the lattice, that is, the total number of interactions. For a lattice with toroidal boundary conditions p= 1 N q, 2 (5.152) CHAPTER 5. MAGNETIC SYSTEMS 276 where q is the number of nearest neighbor sites of a given site; q = 2 for one dimension and q = 4 for a square lattice. To make this procedure explicit, consider an Ising chain with toroidal boundary conditions for N = 3. For this case p = 3(2)/2 = 3, and there are three factors in the product in (5.151): (1 + vs1 s2 )(1 + vs2 s3 )(1 + vs3 s1 ). We can expand this product in powers of v to obtain the 2p = 8 terms in the partition function:
1 1 1 ZN =3 = (cosh βJ )3
s1 =−1 s2 =−1 s3 =−1 1 + v (s1 s2 + s2 s3 + s3 s1 ) (5.153) + v 2 (s1 s2 s2 s3 + s1 s2 s3 s1 + s2 s3 s3 s1 ) + v 3 s1 s2 s2 s3 s3 s1 . It is convenient to introduce a onetoone correspondence between each of the eight terms in the brackets in (5.153) and a diagram. The set of eight diagrams is shown in Figure 5.20. Because v enters into the product in (5.153) as vsi sj , a diagram of order v n has n v bonds. We can use the topology of the diagrams to help us to keep track of the terms in (5.153). The term of order v 0 is 2N =3 = 8. Because si =±1 si = 0, each of the terms of order v vanishes. Similarly, each of the three terms of order v 2 contains at least one of the spin variables raised to an odd power so that these terms also vanish. For example, s1 s2 s2 s3 = s1 s3 , and both s1 and s3 enter to ﬁrstorder. In general, we have 1 2 n even, si n = (5.154) 0 n odd. s =− 1
i From (5.154) we see that only terms of order v 0 and v 3 contribute so that ZN =3 = cosh3 βJ [8 + 8v 3 ] = 23 (cosh3 βJ + sinh3 βJ ). (5.155) We can now generalize our analysis to arbitrary N . We have observed that the diagrams that correspond to nonvanishing terms in Z are those that have an even number of bonds from each vertex; these diagrams are called closed. A bond from site i corresponds to a product of the form si sj . An even number of bonds from site i implies that si raised to an even power enters into the sum in (5.151). Hence, only diagrams with an even number of bonds from each vertex yield a nonzero contribution to ZN . For the Ising chain only two bonds can come from a given site. Hence, although there are 2N diagrams for a Ising chain of N spins with toroidal boundary conditions, only the diagrams of order v 0 (with no bonds) and of order v N contribute to ZN . We conclude that ZN = (cosh βJ )N [2N + 2N v N ]. (5.156) Problem 5.26. The form of ZN in (5.156) is not the same as the form of ZN given in (5.39). Use the fact that v < 1 and take the thermodynamic limit N → ∞ to show the equivalence of the two results for ZN . Problem 5.27. High temperature expansion for four spins Draw the diagrams that correspond to the nonvanishing terms in the high temperature expansion of the partition function for the N = 4 Ising chain. CHAPTER 5. MAGNETIC SYSTEMS
1 277 v0
2 3 v1 v2 v3
Figure 5.20: The eight diagrams that correspond to the eight terms in the partition function for the N = 3 Ising chain. The term si sj is represented by a line between the neighboring sites i and j [see Stanley (1971)]. It is possible to generalize the diagrammatic analysis to higher dimensions. The results of low temperature and high temperature expansions have been used to estimate the values of the various critical exponents (see Domb). An analogous diagrammatic expansion is discussed in Chapter 8 for particle systems. 5.10.4 *Bethe approximation In Section 5.7 we introduced a simple meanﬁeld theory of the Ising model. In the following we discuss how to improve this approximation.16 The idea is that instead of considering a single spin, we consider a group or cluster of spins and the eﬀective ﬁeld experienced by it. In particular, we will choose the group to be a spin and its q nearest neighbors (see Figure 5.21). The interactions of the nearest neighbors with the central spin are calculated exactly, and the rest of the spins in the system are assumed to act on the nearest neighbors through a selfconsistent eﬀective ﬁeld. The energy of the cluster is
q q Hc = −Js0 j =1 sj − Hs0 − Heﬀ
q sj
j =1 (5.157a) = −(Js0 + Heﬀ )
16 This j =1 sj − Hs0 . (5.157b) approach is due to Bethe, who received a Nobel Prize for his work on the theory of stellar nucleosynthesis. CHAPTER 5. MAGNETIC SYSTEMS
s4 278 s1 s0 s2 s3 Figure 5.21: The simplest cluster on the square lattice used in the Bethe approximation. The interaction of the central spin with its q = 4 nearest neighbors is treated exactly. For a square lattice q = 4. Note that the ﬂuctuating ﬁeld acting on the nearest neighbor spins s1 , . . . , sq has been replaced by the eﬀective ﬁeld Heﬀ . The cluster partition function Zc is given by Zc =
s0 =±1,sj =±1 e−βHc . (5.158) We ﬁrst do the sum over s0 = ±1 using (5.157b) and write Zc = eβH
sj =±1 eβ (J +Heff )( Pq j =1 sj ) + e−βH
sj =±1 eβ (−J +Heff )( Pq j =1 sj ) . (5.159) For simplicity, we will evaluate the partition function of the cluster for the onedimensional Ising model for which q = 2. Because the two neighboring cluster spins can take the values ↑↑, ↑↓, ↓↑, and ↓↓, the sums in (5.159) yield Zc = eβH e2β (J +Heff ) + 2 + e−2β (J +Heff ) + e−βH e2β (−J +Heff ) + 2 + e−2β (−J +Heff ) =4 e
βH (5.160a) (5.160b) cosh β (J + Heﬀ ) + e 2 −βH cosh β (J − Heﬀ ) . 2 The expectation value of the central spin is given by s0 = 1 ∂ ln Zc 4 βH = e cosh2 β (J + Heﬀ ) − e−βH cosh2 β (J − Heﬀ ) . β ∂H Zc (5.161) In the following we will set H = 0 to ﬁnd the critical temperature. We also want to calculate the expectation value of the spin of the nearest neighbors sj for j = 1, . . . , q . Because the system is translationally invariant, we require that s0 = sj and ﬁnd the eﬀective ﬁeld Heﬀ by requiring that this condition be satisﬁed. From (5.159) we see that sj = 1 ∂ ln Zc . q ∂ (βHeﬀ ) (5.162) CHAPTER 5. MAGNETIC SYSTEMS If we substitute (5.160b) for Zc in (5.162) with H = 0, we ﬁnd sj = 4 sinh β (J + Heﬀ ) cosh β (J + Heﬀ ) Zc − sinh β (J − Heﬀ ) cosh β (J − Heﬀ ) . 279 (5.163) The requirement s0 = sj yields the relation cosh2 β (J + Heﬀ ) − cosh2 β (J − Heﬀ ) = sinh β (J + Heﬀ ) cosh β (J + Heﬀ ) − sinh β (J − Heﬀ ) cosh β (J − Heﬀ ). (5.164) Equation (5.164) can be simpliﬁed by writing sinh x = cosh x − e−x with x = β (J ± Heﬀ ). The result is cosh β (J + Heﬀ ) = e2βHeff . (5.165) cosh β (J − Heﬀ ) We can follow a similar procedure to ﬁnd the generalization of (5.165) to arbitrary q . The result is coshq−1 β (J + Heﬀ ) = e2βHeff . (5.166) coshq−1 β (J − Heﬀ ) Equation (5.166) always has the solution Heﬀ = 0 corresponding to the high temperature phase. Is there a nonzero solution for Heﬀ for low temperatures? As Heﬀ → ∞, the lefthand side of (5.166) approaches e2βJ (q−1) , a constant independent of Heﬀ , and the righthand side diverges. Therefore, if the slope of the function on the left at Heﬀ = 0 is greater than 2β , the two functions must intersect again at ﬁnite Heﬀ . If we take the derivative of the lefthand side of (5.166) with respect to Heﬀ and set it equal to 2β , we ﬁnd that the condition for a solution to exist is coth βc J = q − 1, (5.167) On the square lattice (q = 4) the condition (5.167) yields kTc /J ≈ 2.885 in comparison to the Onsager solution kTc /J ≈ 2.269 (see [5.86)] and the result of simple meanﬁeld theory, kTc /J = 4. For the onedimensional Ising model (q = 2), the Bethe approximation predicts Tc = 0 in agreement with the exact result. That is, the Bethe approximation does not predict a phase transition in one dimension. Better results can be found by considering larger clusters. Although such an approach yields more accurate results for Tc , it yields the same meanﬁeld exponents because it depends on the truncation of correlations beyond a certain distance. Hence, this approximation must break down in the vicinity of a critical point where the correlation between spins becomes inﬁnite. Problem 5.28. The Bethe approximation (a) Work out the details of the Bethe approximation for the cluster in Figure 5.21 and derive (5.166). (b) Derive (5.167) for the critical temperature. (c) Show that kTc /J ≈ 2.885 for q = 4 and kTc /J = 0 for q = 2. where coth x = cosh x/ sinh x. Because (5.166) is invariant under Heﬀ → −Heﬀ , there will be two solutions for T ≤ Tc . CHAPTER 5. MAGNETIC SYSTEMS 280 5.10.5 Fully connected Ising model We expect that meanﬁeld theory becomes exact for a system in which every spin interacts equally strongly with every other spin in the system because the ﬂuctuations of the eﬀective ﬁeld go to zero in the limit N → ∞. We will refer to this model as the fully connected Ising model. For such a system the energy is given by (see Problem 5.29) E= JN (N − M 2 ), 2 (5.168) where M is the magnetization and JN is the interaction between any two spins. Note that E depends only on M . In Problem 5.29 we also ﬁnd that the number of states with magnetization M is given by N! Ω(M ) = , (5.169) n!(N − n)! where n is the number of up spins. As before, n = N/2 + M/2 and N − n = N/2 − M/2. Problem 5.29. Energy and density of states of the fully connected Ising model (a) Show that the energy of a system for which every spin interacts with every other spin is given by (5.168). A straightforward way to proceed is to consider a small system, say N = 9, and determine the energy of various microstates. As you do so, you will see how to generalize your results to arbitrary N . (b) Use similar considerations as in part (a) to ﬁnd the number of states as in (5.169). The energy of interaction JN of two spins has to scale as 1/N so that the total energy of N spins will be proportional to N . We will choose JN = q J . N (5.170) The factor of q is included so that we will obtain the usual meanﬁeld result for Tc . Given the form of the energy in (5.168) and the number of states in (5.169), we can write the partition function as ZN =
M N 2 + M 2 N! ! N− 2 M 2 ! e−βJN (N −M 2 )/2 βHM e . (5.171) We have included the interaction with an external magnetic ﬁeld. For N ≫ 1 we can convert the sum to an integral. We write
∞ ZN =
−∞ Z (M ) dM, (5.172) and Z (M ) = N! e−βE eβHM , n!(N − n)! (5.173) CHAPTER 5. MAGNETIC SYSTEMS 281 where n = (M + N )/2. A plot of Z (M ) shows that it is peaked about a particular value of M . So let us do our usual trick of expanding ln ZM about its maximum. We ﬁrst ﬁnd the value of M for which Z (M ) is a maximum. We write ln Z (M ) = ln N ! − ln n! − ln(N − n)! − βE + βHM. (5.174) We then use the weaker form of Stirling’s approximation (3.102), d(ln x!)/dx = ln x, dn/dM = 1/2, and d(N − n)/dM = −1/2. The result is 1 1 d ln Z (M ) = − ln n + ln(N − n) + βJN M + βH dM 2 2 1N 1N = − ln (1 + m) + ln (1 − m) + qβJm + βH 2 2 2 2 1 1 = − ln(1 + m) + ln(1 − m) + qβJm + βH = 0, 2 2 (5.175a) (5.175b) (5.175c) where m = M/N . We set d(ln Z (M ))/dM = 0 to ﬁnd the value of m that maximizes Z (M ). We have 1 1−m ln = −β (qJm + H ), 2 1+m so that 1−m = e−2β (qJm+H ) ≡ x. 1+m (5.177) (5.176) Finally we solve (5.177) for m in terms of x and obtain 1 − m = x(1 + m), m(−1 − x) = −1 + x. Hence m= 1−x 1 − e−2β (Jqm+H ) = −2β (Jqm+H ) 1+x e +1 eβ (Jqm+H ) − e−β (Jqm+H ) = −β (Jqm+H ) e + eβ (Jqm+H ) = tanh(β (Jqm + H ). (5.178a) (5.178b) (5.178c) Note that (5.178c) is identical to the meanﬁeld result in (5.108).17
∗ Problem 5.30. Fully connected Ising form of Z (a) Show that Z (M ) can be written as a Gaussian and then do the integral over M in (5.172) to ﬁnd the meanﬁeld form of Z . (b) Use the form of Z from part (a) to ﬁnd the meanﬁeld result for the free energy F . Compare your result to (5.126).
17 Meanﬁeld theory corresponds to taking the limit N → ∞ before letting the range of interaction go to inﬁnity. In contrast, the fully connected Ising model corresponds to taking both limits simultaneously. Although the Ising model gives the same results for the partition function as meanﬁeld theory, the fully connected Ising model can yield diﬀerent results in some other contexts. CHAPTER 5. MAGNETIC SYSTEMS 282 5.10.6 Metastability and nucleation To introduce the concepts of metastability and nucleation we ﬁrst consider the results of the simulations in Problem 5.31. Problem 5.31. Simulations of metastability (a) Use Program Ising2D to simulate the Ising model on a square lattice. Choose L = 64, T = 1, and H = 0.7. Run the simulation until the system reaches equilibrium. You will notice that most of the spins are aligned with the magnetic ﬁeld. (b) Pause the simulation and let H = −0.7; we say that we have “ﬂipped” the ﬁeld. Continue the simulation after the change of the ﬁeld and watch spins. Do the spins align themselves with the magnetic ﬁeld immediately after the ﬂip? Monitor the magnetization of the system as a function of time. Is there an interval of time for which the mean value of m does not change appreciably? At what time does m change sign? What is the equilibrium state of the system after the change of the ﬁeld? (c) Keep the temperature ﬁxed at T = 1 and decrease the ﬁeld to H = 0.6 and ﬂip the ﬁeld as in part (b). Does m become negative sooner or later compared to H  = 0.7? Is this time the same each time you do the simulation? (The program uses a diﬀerent random number seed each time it is run.) You probably found that the spins did not immediately ﬂip to align themselves with the magnetic ﬁeld. Instead most of the spins remained up and the mean values of the magnetization and energy did not change appreciably for many Monte Carlo steps per spin. We say that the system is in a metastable state. The reason that the spins do not ﬂip as soon as the ﬁeld is ﬂipped is that if the ﬁeld is not too large, it costs energy for a spin to ﬂip because it would likely no longer be parallel with its neighbors. If we wait long enough, we will see isolated “droplets” of spins pointing in the stable (down) direction. If a droplet is too small, it will likely shrink and vanish. In contrast, if the droplet is bigger than a critical size (see Figure 5.22), it will grow and the system will quickly reach its stable equilibrium state. If the droplet has a certain critical size, then it will grow with a probability of 50%. This droplet is called the nucleating droplet. The initial decay of the metastable state is called nucleation. Metastable states occur often in nature and in the laboratory. For example, if you take a container of distilled (very pure) water with no dirt, pollen, or other impurities, you can supercool it below the freezing temperature of 0◦ C. The supercooled water will remain a liquid unless there is a spontaneous density ﬂuctuation. More likely, an external disturbance will create the necessary ﬂuctuation. Search <youtube.com> for supercooled water to see some great demonstrations. Metastable states are important in forming crystalline metals from a molten liquid as well as in biological systems and the inﬂationary scenario of the early universe. To determine the size of the nucleating droplet consider nucleation at low temperatures so that we can ignore the entropy.18 A compact droplet (circular in two dimensions and spherical in three dimensions) minimizes the energy cost of creating a droplet of down spins. The energy is
18 At higher temperatures we would consider the free energy. CHAPTER 5. MAGNETIC SYSTEMS 283 Figure 5.22: Example of a nucleating droplet. The simulation was done for the Ising model on a square lattice with L = 64, T = 1, and H  = 0.7. The magnetic ﬁeld was originally up (in the direction of the lighter sites). The nucleating droplet appeared at t ≈ 50 mcs after the ﬂip of the ﬁeld and is the largest cluster of down (dark) sites. decreased by aligning the spins of the droplet with the ﬁeld. This energy decrease is proportional to the area (volume in three dimensions) of the droplet. Hence, Ebulk = −aHrd , (5.179) where r is the radius of the droplet, d is the spatial dimension, and a is a constant. Creating a surface costs energy. The associated energy cost is proportional to the circumference (surface area in three dimensions) of the droplet, and hence Esurf = σrd−1 , (5.180) where σ is the energy cost per unit length (per unit area in three dimensions). This quantity is known as the surface tension. The total energy cost of creating a droplet of radius r is E (r) = −aHrd + σrd−1 . (5.181) The energy cost of the droplet increases as a function of r until a critical radius rc (see Figure 5.23). The radius of the nucleating droplet can be obtained by determining the value of r for which E (r) has a maximum: dE d d = −adHrc −1 + (d − 1)σrc −2 = 0, (5.182) dr r=rc or − adHrc + (d − 1)σ = 0, (5.183) CHAPTER 5. MAGNETIC SYSTEMS 284 4 energy cost 2 0 –2 –4 0 1 2 3 r 4 5 Figure 5.23: The energy cost E (r) of a droplet of radius r from (5.181) for d = 3 (with a = 1, σ = 0.5, and H = 0.1 chosen for convenience). and (d − 1)σ . (5.184) adH The energy cost of creating the nucleating droplet is Ec = bσ d /(aH )d−1 , where b depends on d. The probability of creating the nucleating droplet is proportional to e−βEc . The lifetime of the metastable state, that is, the time before the nucleating droplet occurs, is proportional to the inverse of this probability. rc = Note that we used equilibrium considerations to estimate the radius of the nucleating droplet and the lifetime of the metastable state. This assumption of equilibrium is justiﬁed only if the lifetime of the metastable state is long. Hence, we must have Ec /kT ≫ 1, that is, small ﬁelds or low temperatures. Vocabulary
magnetization m, susceptibility χ free energies F (T, H ), G(T, M ) Ising model, exchange constant J , domain wall spinspin correlation function G(r), correlation length ξ order parameter, ﬁrstorder transition, continuous phase transition, critical point critical temperature Tc , critical exponents α, β , δ , γ , ν , η exact enumeration meanﬁeld theory, Bethe approximation low and high temperature expansions hysteresis, metastable state, nucleating droplet CHAPTER 5. MAGNETIC SYSTEMS
z 285 θ y φ x Figure 5.24: The direction of µ is determined by the angles θ and φ of a spherical coordinate system. Additional Problems
Problem 5.32. Classical paramagnet The energy of interaction of a classical magnetic dipole with the magnetic ﬁeld B is given by E = −µ · B. In the absence of an external ﬁeld the dipoles are randomly oriented so that the mean magnetization is zero. The goal of this problem is to ﬁnd the mean magnetization as a function of B and T . The direction of the magnetization is parallel to B. The sum over microstates becomes an integral over all directions of µ. The direction of µ in three dimensions is given by the angles θ and φ of a spherical coordinate system as shown in Figure 5.24. The integral is over the solid angle element dΩ = sin θdθdφ. In this coordinate system the energy of the dipole is given by E = −µB cos θ. (a) Choose spherical coordinates and show that the probability p(θ, φ)dθdφ that the dipole is between the angles θ and θ + dθ and φ and φ + dφ is given by p(θ, φ)dθdφ = where Z1 is given by Z1 =
0 0 eβµB cos θ sin θdθdφ, Z1 (5.185) 2π π eβµB cos θ sin θ dθ dφ. (5.186) (b) How is cos θ related to Z1 ? CHAPTER 5. MAGNETIC SYSTEMS (c) Show that the mean magnetization is given by M = N µL β µB , where the Langevin function L(x) is given by L(x) = 1 1 e x + e −x − = coth x − . e x − e −x x x 286 (5.187) (5.188) (d) For x < π , L(x) can be expanded as L(x) = x x3 22n B2n − + ...+ + ... 3 45 (2n)! (x ≪ 1), (5.189) where Bn is the Bernoulli number of order n (see the Appendix). What is M and the susceptibility in the limit of high T ? (e) For large x, L(x) is given by L(x) ≈ 1 − 1 + 2e−2x x (x ≫ 1). (5.190) What is the behavior of M in the limit of low T ? (f) What is the mean energy and the entropy of a system of N noninteracting magnetic dipoles? Is the behavior of the entropy at low temperatures consistent with the third law of thermodynamics? Problem 5.33. Arbitrary spin The magnetic moment of an atom or nucleus is associated with its angular momentum which is quantized. If the angular momentum is J , the magnetic moment along the direction of B is restricted to (2J + 1) orientations. We write the energy of an individual atom as E = −gµ0 J · B = −gµ0 Jz B. (5.191) The values of µ0 and g depend on whether we are considering a nucleus, an atom, or an electron. The values of Jz are restricted to −J , −J + 1, −J + 2, . . . , J − 1, J . Hence, the partition function for one atom contains (2J + 1) terms:
J Z1 =
m =− J e−β (−gµ0 mB ) . (5.192) The summation index m ranges from −J to J in steps of +1. To simplify the notation, we let α = βgµ0 B , and write Z1 as a ﬁnite geometrical series:
J Z1 =
m =− J emα (5.193a) (5.193b) = e−αJ (1 + eα + e2α + . . . + e2Jα ). CHAPTER 5. MAGNETIC SYSTEMS The sum of a ﬁnite geometrical series is given by
n 287 Sn =
p=0 xp = xn+1 − 1 . x−1 (5.194) Given that there are (2J + 1) terms in (5.193b), show that Z1 = e−αJ e(2J +1)α − 1 [1 − e(2J +1)α ] . = e−αJ eα − 1 1 − eα M = N gµ0 JBJ (α), where the Brillouin function BJ (α) is deﬁned as BJ (α) = 1 1 (J + 1/2) coth(J + 1/2)α − coth α/2 . J 2 (5.197) (5.195) Use the above relations to show that (5.196) What is the limiting behavior of M for high and low T for ﬁxed B ? What is the limiting behavior of M for J = 1/2 and J ≫ 1? Problem 5.34. Density of states In Problem 4.40 the density of states was given without proof for the onedimensional Ising model for even N and toroidal boundary conditions: Ω(E, N ) = 2 N i =2 N! i! (N − i)! (i = 0, 2, 4, . . . , N ), (4.18)
∗ with E = 2 i − N . Use this form of Ω and the relation ZN =
E Ω(E, N )e−βE (5.198) to ﬁnd the partition function for small values of (even) N . Problem 5.35. Sample microstates The ﬁve microstates shown in Figure 5.25 for the Ising chain were generated using the Metropolis algorithm (see Sections 4.11 and 5.5.3) at βJ = 2 using toroidal boundary conditions. On the basis of this limited sample, estimate the mean value of E/J and the magnetization. Calculate the spinspin correlation function G(r) for r = 1, 2, and 3 using the third spin as the origin, and then repeat the calculation using the sixth spin. Remember that the Metropolis algorithm simulates a system in equilibrium with a heat bath at temperature T with the correct weight. Explain why your results are not accurate. Problem 5.36. Enumeration of microstates of the twodimensional Ising model CHAPTER 5. MAGNETIC SYSTEMS 288 (a) (b) (c) (d) (e) Figure 5.25: Five microstates of the Ising chain with N = 10 spins with toroidal boundary conditions generated by the Metropolis algorithm at βJ = 2 and H = 0. (a) Calculate the partition function for the Ising model on a square lattice for N = 4 spins in the presence of an external magnetic ﬁeld H . Assume that the system is in equilibrium with a heat bath at temperature T . Choose toroidal boundary conditions. (b) Determine Ω(E ), the number of states with energy E , and discuss its dependence on E . Assume that H = 0. (c) The 42 = 16 microstates of the twodimensional Ising model for N = 4 can be grouped into four “ordered” states with energies ±J and 12 “disordered” states with zero energy. Test the hypothesis that the phase transition occurs when the partition function of the disordered states equals that of the ordered states. What is the resulting value of Tc ? This simple reasoning does not work as well for the Ising model in three dimensions. Problem 5.37. Form of P (E ) for the Ising model Consider the twodimensional Ising model in equilibrium with a heat bath at temperature T . (a) On the basis of general considerations, what is the form of the probability P (E )∆E that the system has energy between E and E + ∆E ? (b) Why is this form not applicable at the critical point? Problem 5.38. The Ising model and cooperative phenomena Explore the analogy between the behavior of the Ising model and the behavior of a large group of people. Under what conditions would a group of people act like a collection of individuals, each doing their “own thing?” Under what conditions might they act as a group? What factors could CHAPTER 5. MAGNETIC SYSTEMS 289 cause a transition from one behavior to the other? The relation of the Ising model to models of economic opinions, urban segregation, and language change is discussed by Stauﬀer. Problem 5.39. The demon and the Ising chain (a) Consider a demon that exchanges energy with the Ising chain by ﬂipping single spins. Show that the possible changes in the energy in zero magnetic ﬁeld are 0 and ±4J . Conﬁrm that the possible demon energies are Ed = 4nJ , where n = 0,, 1, 2, . . . (b) Derive an expression for the mean demon energy as a function of the temperature of the system.
∗ Problem 5.40. Applications of the transfer matrix method (a) Consider a onedimensional Isingtype model with si = 0, ±1. Use the transfer matrix method to calculate the dependence of the energy on T for H = 0. The solution requires the diﬀerentiation of the root of a cubic equation that you might wish to do numerically. (b) Use the transfer matrix method to ﬁnd the thermodynamic properties of the q = 3 Potts model in one dimension. Problem 5.41. Finitesize scaling and critical exponents Although a ﬁnite system cannot exhibit a true phase transition characterized by divergent physical quantities, we expect that if the correlation length ξ (T ) is less than the linear dimension L of the system, our simulations will yield results comparable to an inﬁnite system. However, if T is close to Tc , the results of simulations will be limited by ﬁnitesize eﬀects. Because we can only simulate ﬁnite lattices, it is diﬃcult to obtain estimates for the critical exponents α, β , and γ by using their deﬁnitions in (5.95), (5.97), and (5.98) directly. The eﬀects of ﬁnite system size can be made quantitative by the following argument, which is based on the fact that the only important length near the critical point is the correlation length. Consider, for example, the critical behavior of χ. If the correlation length ξ ≫ 1,19 but is much less than L, the power law behavior given by (5.98) is expected to hold. However, if ξ is comparable to L, ξ cannot change appreciably and (5.98) is no longer applicable. This qualitative change in the behavior of χ and other physical quantities occurs for ξ ∼ L ∼ T − Tc −ν . We invert (5.199) and write Hence, if ξ and L are approximately the same size, we can replace (5.98) by the relation χ(T = Tc ) ∼ [L−1/ν ]−γ ∼ Lγ/ν . (5.201) T − Tc  ∼ L−1/ν . (5.200) (5.199) The relation (5.201) between χ and L at T = Tc is consistent with the fact that a phase transition is deﬁned only for inﬁnite systems. We can use the relation (5.201) to determine the ratio γ/ν . This method of analysis is known as ﬁnite size scaling.
19 All lengths are measured in terms of the lattice spacing. CHAPTER 5. MAGNETIC SYSTEMS 290 (a) Use Program Ising2D to estimate χ at T = Tc for diﬀerent values of L. Make a loglog plot of χ versus L and use the scaling relation (5.201) to determine the ratio γ/ν . Use the exact result ν = 1 to estimate γ . Then use the same reasoning to determine the exponent β and compare your estimates for β and γ with the exact values given in Table 5.1. (b) Make a loglog plot of C versus L. If your data for C is suﬃciently accurate, you will ﬁnd that the loglog plot of C versus L is not a straight line but shows curvature. The reason is that the exponent α equals zero for the twodimensional Ising model, and C ∼ C0 ln L. Is your data for C consistent with this form? The constant C0 is approximately 0.4995. Problem 5.42. Low temperature behavior in meanﬁeld theory (a) Write (5.108) in the form βqJm = tanh−1 m = (1/2) ln[(1 + m)/(1 − m)] and show that m(T ) ≈ 1 − 2e−βqJ as T → 0. (b) Determine the low temperature behavior of χ. Does it approach zero for T ≪ Tc ? Problem 5.43. Veriﬁcation of (5.149) Verify the validity of the identity (5.149) by considering the diﬀerent possible values of si sj and using the identities 2 cosh x = ex + e−x and 2 sinh x = ex − e−x .
∗ (5.202) Problem 5.44. Lifetime of metastable state In Section 5.10.6 we discussed some of the features of metastable states in the Ising model. Suppose that we ﬂip the magnetic ﬁeld as in Problem 5.31 and deﬁne the lifetime of the metastable state as the number of Monte Carlo steps per spin from the time of the change of the ﬁeld to the occurrence of the nucleating droplet. Because the system can be treated by equilibrium considerations while it is in a metastable state, the probability of occurrence of the nucleating droplet during any time interval ∆t is independent of time. What is the form of the probability p(t)∆t that the lifetime of the metastable state is between t and t + ∆t (see Section 3.9)? Suggestions for Further Reading
Stephen G. Brush, “History of the LenzIsing model,” Rev. Mod. Phys. 39, 883–893 (1967). Cyril Domb, The Critical Point, Taylor & Francis (2996). This graduate level monograph discusses the history of our understanding of phase transitions. Much of the history is accessible to undergraduates. David L. Goodstein, States of Matter, Dover Publications (1985). This book discusses how to apply thermodynamics and statistical mechanics to gases, liquids, and solids, and reviews magnetism in solids. S. Kobe, “Ernst Ising 1900–1998,” Braz. J. Phys. 30 (4), 649–654 (2000). Available online from several web sites. CHAPTER 5. MAGNETIC SYSTEMS 291 D. P. Landau, ShanHo Tsai, and M. Exler, “A new approach to Monte Carlo simulations in statistical physics: WangLandau sampling,” Am. J. Phys. 72, 1294–1302 (2004). B. Liu and M. Gitterman, “The critical temperature of twodimensional and threedimensional Ising models,” Am. J. Phys. 71, 806–808 (2003). The authors give a simple argument for the value of the critical temperature of the Ising model which gives the exact result in two dimensions and an approximate result in three dimensions [see Problem 5.36(c)]. F. Mandl, Statistical Physics, second edition, John Wiley & Sons (1988). This text has a good discussion of the thermodynamics of magnetism. Daniel C. Mattis, The Theory of Magnetism Made Simple, World Scientiﬁc (2006). The Bohrvan Leeuwen theorem is discussed in Section 1.6. Peter PalﬀyMuhoraya, “The single particle potential in meanﬁeld theory,” Am. J. Phys. 70, 433–437 (2002). M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics, Clarendon Press (1999). H. Eugene Stanley, Introduction to Phase Transitions and Critical Phenomena, Oxford University Press (1971). The discussion of the high temperature expansion in Section 5.10.3 is based in part on this book. D. Stauﬀer, “Social applications of twodimensional Ising models,” Am. J. Phys. 76, 470–473 (2008). Jan Tobochnik and Harvey Gould, “Teaching statistical physics by thinking about models and algorithms,” Am. J. Phys. 76, 353–359 (2008). ...
View Full
Document
 Summer '10
 hcverma

Click to edit the document details