A Student’s Guide to the Ising Model

Chapter 1

Figure 1.7
Boltzmann factor

Clear[f,En,kB,T]

Define the Boltzmann factor.

f=(E^(-En/(kB T)))

Plot the Boltzmann factor as a function of (Subscript[E, n])for constant T. Set (Subscript[k, B]) = 1.

Plot[f/.{kB->1,T->1},{En,0,4},AxesLabel->{En,f}]

Plot the Boltzmann factor as a function of T for constant (Subscript[E, n]). Set (Subscript[k, B]) = 1.

Plot[{f/.{kB->1,En->1},1},{T,0,15},PlotRange->{0,1},PlotStyle->{{},{Dashed}},AxesLabel->{T,f}]

Figure 1.8
Probabilities for two-state system

Clear[P1,P2,E1,kB,T]

Probabilities P1 and P2 for a two-state system, with energies E1 and E2.

P1=((E^(-E1/(kB T)))/((E^(-E1/(kB T)))+(E^(-E2/(kB T)))))

P2=((E^(-E2/(kB T)))/((E^(-E1/(kB T)))+(E^(-E2/(kB T)))))

Set E1 = 1 and E2 = 2, so that E2 > E1. Also, set kB = 1.

At low temperature the system is in state 1, the low-energy state, with a probability essentially equal to 1. At high temperature the system is equally likely to be found in one of the two states, state 1 or state 2.

Plot[{P1/.{E1->1,E2->2,kB->1},P2/.{E1->1,E2->2,kB->1},0.5},{T,0,5},PlotStyle->{{},{},{Dashed}},AxesLabel->{T,Probability}]

Chapter 2

Figure 2.2
Single-spin system with magnetic field

Clear[f,h]

Define the reduced free energy per site, f, for a single-spin system with magnetic field h, as given in Equation 2.4.

f=Log[2]+Log[Cosh[h]]

Plot the free energy as a function of h. Add limiting values (dashed lines) for h goes to zero and h goes to infinity.

Plot[{f,Log[2],h},{h,0,2.5},PlotStyle->{{},{Dashed},{Dashed}},AxesLabel->{h,"f"}]

Figure 2.6, Figure 2.9
Two-spin system with coupling K

Clear[f,K,c]

Define the reduced free energy per site, f, for a two-spin system with nearest-neighbor coupling K, as given in Equation 2.12.

f=(1/2)Log[4Cosh[K]]

Plot the free energy as a function of K. Also, plot the small-K limit of Log[2], and the large-K slope of 1/2, with dashed lines.

Plot[{f,Log[2],(1/2)K},{K,0,5.1},PlotStyle->{{},{Dashed},{Dashed}},AxesLabel->{K,"f"}]

Define the reduced specific heat per site, c, for a two-spin system with nearest-neighbor coupling K, as given in Equation 2.17.

c=(K^2)(∂_(K,K))f

Plot the specific heat as a function of K.

Plot[c,{K,0,5.1},AxesLabel->{K,"c"}]

Figure 2.11, Figure 2.13
Three-spin system with free boundary conditions

Clear[f,K,c]

Define the reduced free energy per site, f, for a three-spin system with free boundary conditions and nearest-neighbor coupling K, as given on Page 48.

f=(1/3)Log[8(((Cosh[K]))^2)]

Plot the free energy as a function of K.

Plot[{f,-(2/3)K,(2/3)K},{K,-2.1,2.1},PlotStyle->{{},{Dashed},{Dashed}},AxesLabel->{K,"f"},PlotRange->{0,2.1}]

Define the reduced specific heat per site, c, for a three-spin system with free boundary conditions and nearest-neighbor coupling K. This result is given on Page 49.

c=(K^2)(∂_(K,K))f

Plot the specific heat as a function of K.

Plot[c,{K,-4.1,4.1},AxesLabel->{K,"c"}]

Figure 2.17, Figure 2.19
Three-spin system with periodic boundary conditions

Clear[f,K,c]

Define the reduced free energy per site, f, for a three-spin system with periodic boundary conditions and nearest-neighbor coupling K, as given in Equation 2.20.

f=(1/3)Log[2((E^(3K))+3(E^(-K)))]

Plot the free energy as a function of K.

Plot[{f,-(1/3)K,K},{K,-2.1,2.1},PlotStyle->{{},{Dashed},{Dashed}},AxesLabel->{K,"f"},PlotRange->{0,2.1}]

Define the reduced specific heat per site, c, for a three-spin system with periodic boundary conditions and nearest-neighbor coupling K.

c=(K^2)(∂_(K,K))f

Plot the specific heat as a function of K.

Plot[c,{K,-4.1,4.1},AxesLabel->{K,"c"},PlotRange->{0,0.41}]

Figure 2.22
The 2 X 2 system

Clear[H22,K,s1,s2,s3,s4,Z22,f22,c22]

Define the reduced Hamiltonian, H, for the 2 X 2 system with periodic boundary conditions and nearest-neighbor coupling K and call it H22. Refer to Figure 2.21 for the numbering scheme used for this system.

H22=K(s1(s4+s2)+s2(s3+s1)+s3(s2+s4)+s4(s1+s3))

Sum Boltzmann factors over all configurations to find the partial partition function, Z22.

Z22=Sum[(E^H22),{s1,-1,1,2},{s2,-1,1,2},{s3,-1,1,2},{s4,-1,1,2}]

Define the reduced free energy per site, f, for the 2 X 2 system. Call it f22.

f22=(1/4)Log[Z22]

Plot the free energy as a function of K. Add limiting values (dashed lines) for large and small K.

Plot[{f22,2K,Log[2]},{K,0,1.1},PlotStyle->{{},{Dashed},{Dashed}},AxesLabel->{K,"f"},PlotRange->{0,2.1}]

Define the reduced specific heat per site, c, for the 2 X 2 system system with periodic boundary conditions and nearest-neighbor coupling K, and call it c22.

c22=(K^2)(∂_(K,K))f22

Plot the specific heat as a function of K.

Plot[c22,{K,0,1.41},AxesLabel->{K,"c22"},PlotRange->{0,0.42}]

Figure 2.24
The 3 X 3 system

Clear[H33,K,s1,s2,s3,s4,s5,s6,s7,s8,s9,Z33,f33,c33,c22]

Define the reduced Hamiltonian, H, for the 3 X 3 system with periodic boundary conditions and nearest-neighbor coupling K. Call it H33. Refer to Figure 2.23 for the numbering scheme used for this system.

H33=K(s1(s7+s2)+s2(s8+s3)+s3(s9+s1)+s4(s1+s5)+s5(s2+s6)+s6(s3+s4)+s7(s4+s8)+s8(s5+s9)+s9(s6+s7))

Sum Boltzmann factors over all configurations to find the partial partition function, Z33.

Z33=Sum[(E^H33),{s1,-1,1,2},{s2,-1,1,2},{s3,-1,1,2},{s4,-1,1,2},{s5,-1,1,2},{s6,-1,1,2},{s7,-1,1,2},{s8,-1,1,2},{s9,-1,1,2}]

Define the reduced free energy per site, f, for the 3 X 3 system. Call it f33.

f33=(1/9)Log[Z33]

Plot the free energy as a function of K. Add limiting results with dashed lines.

Plot[{f33,2K,Log[2]},{K,0,1.1},PlotStyle->{{},{Dashed},{Dashed}},AxesLabel->{K,"f"},PlotRange->{0,2.1}]

Define the reduced specific heat per site, c33, for the 3 X 3 system with periodic boundary conditions and nearest-neighbor coupling K.

c33=(K^2)(∂_(K,K))f33

Plot the specific heat as a function of K.

Plot[c33,{K,0,1.41},AxesLabel->{K,"c33"},PlotRange->{0,0.72}]

For the sake of comparison, define once again the reduced specific heat per site, c22, for the 2 X 2 system system with periodic boundary conditions and nearest-neighbor coupling K. This result is given in the code example for Figure 2.22.

c22=(-((((-16 (E^(-8 K))+16 (E^(8 K))))^2)/(4 (((12+2 (E^(-8 K))+2 (E^(8 K))))^2)))+((128 (E^(-8 K))+128 (E^(8 K)))/(4 (12+2 (E^(-8 K))+2 (E^(8 K)))))) (K^2)

Plot the specific heats for both the 2 X 2 and 3 X 3 systems.

Plot[{c22,c33},{K,0,1.41},AxesLabel->{K,"c22, c33"},PlotRange->{0,0.72}]

Chapter 3

Figure 3.4
Effective magnetic field

Clear[hp,K,h,hp1,pl1,hp11,hp12,pl2]

Define the effective magnetic field, hp, as given in Equation 3.10, for K and h. Set K = 1 for these results.

hp=(1/2)Log[(((E^(K+h))+(E^(-K-h)))/((E^(K-h))+(E^(-K+h))))]

hp1=hp/.K->1

Plot the effective magnetic field versus h. Also plot hp1 = h for comparison, and to used later in cob-webbing.

pl1=Plot[{hp1,h},{h,0,2.5},AxesLabel->{h,"hp1"},PlotRange->{0,1.05}]

Add cob-webbing, starting at h = 2, as in Figure 3.4. Do the first couple iterations by hand, and call them hp11, and hp12.

hp11=hp1/.h->2.0

hp12=hp1/.h->hp11

Put together the lines needed to do the first few steps of cob-webbing, and then show them in a plot, along with the effective magnetic field.

pl2=Graphics[{Dashed,Line[{{2,0},{2,hp11}}],Line[{{2,hp11},{hp11,hp11}}],Line[{{hp11,hp11},{hp11,hp12}}],Line[{{hp11,hp12},{hp12,hp12}}]}];

Show[pl1,pl2]

Figure 3.9
Effective interaction

Clear[K,K1,K2,K3,K4,y]

Define the effective interaction between end spins separated by 1, 2, 3, or 4 bonds, as shown in Figure 3.9.

K1=K

K2=(1/2)Log[Cosh[2K]]

K3=(1/2)Log[(((E^(3K))+3(E^(-K)))/((E^(-3K))+3(E^K)))]

K4=(1/2)Log[Cosh[2K2]]

Plot the various effective couplings as a function of K.

pl1=Plot[{K1,K2,K3,K4},{K,0,1.21},AxesLabel->{"K","y"},PlotRange->{0,1.22}]

Chapter 4

Figure 4.4
Specific heats in 1-D, Free boundary conditions

Clear[f,K,c,a]

First, write out the reduced free energy per site, f, for an Ising chain of N sites, as given in Equation 4.3.

f=Log[2]+((N-1)/N)Log[Cosh[K]]

Next, define the reduced specific heat per site, c, and plot as a function of K for various values of N.

c=(K^2)(∂_(K,K))f

Plot[{c/.N->5,c/.N->10,c/.N->10000},{K,0,4.1},AxesLabel->{K,"c"},PlotRange->{0,0.46}]

Now, let’s animate the specific heat for N ranging from 2 to 100. The dashed curve is the specific heat for large N.

Animate[Plot[{(c /. N->a),(c /. N->10000)},{K,0,4.1},PlotRange->{0,0.46},PlotStyle->{Black,Dashed}],{a,2,100},AnimationRunning->False]

Figure 4.14
Specific heats in 1-D, Periodic boundary conditions

Clear[f,K,c,a]

First, write out the reduced free energy per site, f, for an Ising ring of N sites, as given in Equation 4.35.

f=Log[2]+Log[Cosh[K]]+(1/N)Log[1+(((Tanh[K]))^N)]

Next, define the reduced specific heat per site, c, and plot as a function of K for various values of N.

c=(K^2)(∂_(K,K))f

Plot[{c/.N->5,c/.N->15,c/.N->45,c/.N->135,c/.N->10000},{K,0,3.6},AxesLabel->{K,"c"},PlotRange->{0,0.555},PlotStyle->{{},{},{},{},{Dashed}}]

Now, let’s animate the specific heat for N ranging from 2 to 100. The dashed curve is the specific heat for large N.

Animate[Plot[{(c /. N->a),(c /. N->10000)},{K,0,4.1},PlotRange->{0,0.55},PlotStyle->{Black,Dashed}],{a,2,200},AnimationRunning->False]

Figure 4.15
Magnetization per site for infinite lattice

Clear[m,K,h,a]

Write out the magnetization for the infinite lattice, m, as a function of K and h.

m=(((E^K)Sinh[h])/(Sqrt[(E^(-2K))+(E^(2K))(((Sinh[h]))^2)]))

Plot m versus h for various values of the coupling K.

Plot[{m/.K->0,m/.K->1,m/.K->2,m/.K->3,1,-1},{h,-2.1,2.1},AxesLabel->{h,"m"},PlotRange->{-1.1,1.1},PlotStyle->{{},{},{},{},{Dashed},{Dashed}}]

Now, let’s animate the magnetization for K ranging from 0 to 3.

Animate[Plot[{(m /. K->a),1,-1},{h,-2.1,2.1},PlotRange->{-1.1,1.1},PlotStyle->{Black,Dashed,Dashed}],{a,0,3},AnimationRunning->False]

Figure 4.19
Nibble calculation fixed point

Clear[K1p,K]

Write an expression for (Subscript[K, 1])’ in terms of K and (Subscript[K, 1]). Call it K1p.

K1p=K+(1/2)Log[(((E^(-2K+(Subscript[K, 1])))+(E^(2K+(Subscript[K, 1])))+2(E^(-(Subscript[K, 1]))))/((E^(-2K-(Subscript[K, 1])))+(E^(2K-(Subscript[K, 1])))+2(E^(Subscript[K, 1]))))]

Plot K1p versus K1 for K = 1.0.

Plot[{K1p/.K->1,(Subscript[K, 1])},{(Subscript[K, 1]),0,2.1},AxesLabel->{(Subscript[K, 1]),"K1p"},PlotRange->{0,2.1}]

Chapter 5

Figure 5.2
Exact free energy for 2-D Ising model

Clear[u,K,g,th,f,plfreeenergy,Kc,pldashed]

Define an important variable in this calculation.

u=(1/(((Sinh[2K]))^2))

Define the integrand and call it g. The integration variable is th.

g=Log[((Cosh[2K])^2)+(1/u)(Sqrt[1+(u^2)-2u Cos[2th]])]

Define the free energy f in terms of the integrand g as follows:

f = (1/2)Log[2] + (1/(2π))(∫_0^π)gdth

We’ll put this expression into the Plot command. We’ll also plot the K = 0 limit, Log[2], and the K –> ∞ limit, 2K.

plfreeenergy=Plot[{(((Log[2])/2)+((1/(2π)))NIntegrate[g,{th,0,π}]),2K,Log[2]},{K,0,1},PlotStyle->{Black,Dashed,Dashed},AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["f",24,Plain,Black,FontFamily->"Times New Roman"]}]

Add a vertical line to indicate the location of the critical point.

First, define the location of the critical point, (Subscript[K, c]).

Kc=0.5Log[1+(Sqrt[2])]

Now, generate a vertical dashed line at the critical value of K.

pldashed=Graphics[{Dashed,Line[{{Kc,0},{Kc,3}}]}];

Now, combine the graphics.

Show[plfreeenergy,pldashed,AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["f",24,Plain,Black,FontFamily->"Times New Roman"]}]

Figure 5.3
Exact internal energy for 2-D Ising model

Clear[u,K,g,th,dg,e,f,plenergy,Kc,pldashed]

Define an important variable in this calculation.

u=(1/(((Sinh[2K]))^2))

Define the integrand, and call it g. The integration variable is th.

g=Log[((Cosh[2K])^2)+(1/u)(Sqrt[1+(u^2)-2u Cos[2th]])]

Define the derivative of the integrand with respect to K, and call it dg.

dg=(∂_K)g

Recall that the internal energy is K times the derivative of the free energy with respect to K, as in Equation 2.14.

Therefore, we define the internal energy, e, in terms of the derivative of the integrand as follows: 

e = (K/(2π))(∫_0^π)dgdth

Now, plot e, as well as its large-K limit with a dashed line.

plenergy=Plot[{K((1/(2π)))NIntegrate[dg,{th,0,π}],2K},{K,0,1},AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["e",24,Plain,Black,FontFamily->"Times New Roman"]},PlotStyle->{Black,{Gray,Dashed}}]

Add a vertical line to indicate the location of the critical point.

First, define the location of the critical point, (Subscript[K, c]).

Kc=0.5Log[1+(Sqrt[2])]

Now, generate a vertical dashed line at the critical value of K.

pldashed=Graphics[{Dashed,Line[{{Kc,0},{Kc,3}}]}];

Combine the graphics.

Show[plenergy,pldashed,AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["f",24,Plain,Black,FontFamily->"Times New Roman"]}]

Figure 5.4
Exact specific heat for 2-D Ising model

Clear[u,K,g,th,d2g,c,plspecificheat,Kc,pldashed]

Define an important variable in this calculation.

u=(1/(((Sinh[2K]))^2))

Define the integrand and call it g. The integration variable is th.

g=Log[((Cosh[2K])^2)+(1/u)(Sqrt[1+(u^2)-2u Cos[2th]])]

Define the second derivative of the integrand with respect to K, and call it d2g.

d2g=(∂_(K,K))g;

Recall that the reduced specific heat, c, is K^2 times the second derivative of the free energy with respect to K, as in Equation 2.16.

Therefore, we define the specific heat as follows:

c = ((K^2)/(2π))(∫_0^π)d2gdth

We now plot this function.

plspecificheat=Plot[{(K^2)((1/(2π)))NIntegrate[d2g,{th,0,π}]},{K,0,1},PlotRange->{0,2.3},AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["c",24,Plain,Black,FontFamily->"Times New Roman"]},PlotStyle->{Black}]

Add a vertical line to indicate the location of the critical point.

First, define the location of the critical point, (Subscript[K, c]).

Kc=0.5Log[1+(Sqrt[2])]

Now, generate a vertical dashed line at the critical value of K.

pldashed=Graphics[{Dashed,Line[{{Kc,0},{Kc,3}}]}];

Combine the graphics.

Show[plspecificheat,pldashed,AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["c",24,Plain,Black,FontFamily->"Times New Roman"]}]

Figure 5.6
Exact magnetization for 2-D Ising model

Clear[Kc,Tc,m,x,pl1,list,pl2]

Define the critical coupling, (Subscript[K, c]).

Kc=0.5Log[1+(Sqrt[2])]

Define the critical temperature, (Subscript[T, c]).

Tc=1/Kc

Define the magnetization, m, using Equation 5.7. Also, let T/(Subscript[T, c]) = x for simplicity.

m=(((1-(1/(((Sinh[(2/(x Tc))]))^4))))^(1/8))

Plot m versus x. Avoid x = 1 where the singular behavior of m needs to be dealt with carefully.

pl1=Plot[m,{x,0,0.99999999},PlotRange->{0,1.1}]

Add a final segment to include x = 1. First, define the connecting line segment.

list={{1,0},{0.9999999,m/.x->0.9999999}}

Plot this tag end by itself.

pl2=ListPlot[{list},Joined->True,PlotRange->{{0,1.1},{0,1.1}}]

Now, combine the plots.

Show[pl1,pl2,AxesStyle->Directive[Black,20],AxesLabel->{Style["T/Tc",24,Plain,Black,FontFamily->"Times New Roman"],Style["m",24,Plain,Black,FontFamily->"Times New Roman"]}]

Figure 5.8
Specific heat compared with series expansions

Clear[clt1,K,clt2,clt3,plclt,cht1,cht2,cht3,plcht]

Clear[u,K,g,th,d2g,c,plspecificheat,dv,Kc,pldashed]

Define and plot the low-temperature series approximations to the specific heat. Let clt1 be the approximation with one term, clt2 be the approximation with two terms, and clt3 be the approximation with 3 terms.

clt1=(K^2)64((((E^(-2K))))^4)

clt2=(K^2)64((((E^(-2K))))^4)+(K^2)288((((E^(-2K))))^6)

clt3=(K^2)64((((E^(-2K))))^4)+(K^2)288((((E^(-2K))))^6)+(K^2)1152((((E^(-2K))))^8)

plclt=Plot[{clt1,clt2,clt3},{K,0.45,1.1},PlotStyle->{Dashed}]

Define and plot the high-temperature series approximations to the specific heat. Use a similar notation in this case.

cht1=(K^2)(∂_(K,K))(Log[2]+2Log[Cosh[K]])

cht2=(K^2)(∂_(K,K))(Log[2]+2Log[Cosh[K]]+(((Tanh[K]))^4))

cht3=(K^2)(∂_(K,K))(Log[2]+2Log[Cosh[K]]+(((Tanh[K]))^4)+2(((Tanh[K]))^6))

plcht=Plot[{cht1,cht2,cht3},{K,0,0.43},PlotStyle->{Dashed}]

Now, repeat the calculation of the exact specific heat.

Define an important variable in this calculation.

u=(1/(((Sinh[2K]))^2))

Integrand the integrand and call it g. The integration variable is th.

g=Log[((Cosh[2K])^2)+(1/u)(Sqrt[1+(u^2)-2u Cos[2th]])]

Define the second derivative of the integrand with respect to K, and call it d2g.

d2g=(∂_(K,K))g;

Recall that the reduced specific heat, c, is K^2 times the second derivative of the free energy with respect to K, as in Equation 2.16.

Therefore, we define the specific heat as follows:

c = ((K^2)/(2π))(∫_0^π)d2gdth

We now plot this function.

plspecificheat=Plot[{(K^2)((1/(2π)))NIntegrate[d2g,{th,0,π}]},{K,0,1},PlotRange->{0,2.3},AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["c",24,Plain,Black,FontFamily->"Times New Roman"]},PlotStyle->{Black}]

Add a vertical line to indicate the location of the critical point.

First, define the location of the critical point, (Subscript[K, c]).

Kc=0.5Log[1+(Sqrt[2])]

pldashed=Graphics[{Dashed,Line[{{Kc,0},{Kc,3}}]}];

Combine the graphics.

Show[plspecificheat,pldashed,AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["c",24,Plain,Black,FontFamily->"Times New Roman"]}]

Finally, add the low- and high-temperature series expansions to the plot of the exact specific heat.

Show[plspecificheat,pldashed,plclt,plcht,AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["c",24,Plain,Black,FontFamily->"Times New Roman"]}]

Chapter 6

Figure 6.2
Self-consistent relation for single-spin mean-field

Clear[m,mmf,x,a]

Define the mean-field magnetization, mmf, in terms of m and T/(Subscript[T, c]) = x.

Clear[m]

mmf=Tanh[x^(-1) m]

Plot the mean-field magnetization for three different temperatures.

Plot[{mmf/.x->0.5,mmf/.x->1,mmf/.x->1.5,m,1,-1},{m,-1.1,1.1},PlotStyle->{{},{},{},{Dashed},{Dashed},{Dashed}}]

Animate the self-consistent relation for T/(Subscript[T, c]) varying from 2 to 0.1.

Animate[Plot[{(mmf /. x->a),1,-1,m},{m,-1.05,1.05},PlotRange->{{-1.05,1.05},{-1.05,1.05}},PlotStyle->{Black,Dashed,Dashed,Dashed}],{a,2,0.1},AnimationRunning->False]

Table 6.1
Fixed point values m* as a function of T/Tc

Clear[x,mn1,mn2,mn,mn0]

Define the magnetization, (Subscript[m, n + 1]) in terms of (Subscript[m, n]) and T/(Subscript[T, c]) = x.

mn1=Tanh[(x^(-1))mn]

Now, iterate until the fixed point is reached. Set x to whatever value of T/Tc you would like. At the moment, the calculation is set up to find the fixed point for x = 0.6, which we see from Table 6.1 is equal to 0.9073.... We let mn0 = 0.5 be the starting point for the iterations to the fixed point.

mn0=0.5;
Do[
mn2=mn1/.{mn->mn0,x->0.6};
Print["m = ",mn2];
If[Abs[mn2-mn0] < 0.000001,Break[]];
mn0=mn2,
{i,1,100}]

Figure 6.6
Mean-field specific heat

Clear[c,m,x,clist]

Define the specific heat, as given in Equation 6.12, in terms of m and T/(Subscript[T, c]) = x.

c=((m^2)/(((x^2)/((1-(m^2))))-x))

Here is a table of fixed-point values of m obtained with the self-consistent relation for various values of T/(Subscript[T, c]) = x.

T/Tc		m
0		1
0.1		1
0.2		0.9999
0.3		0.9974
0.4		0.9856
0.5		0.9575
0.6		0.9073
0.7		0.8286
0.8		0.7104
0.9		0.5253
1.0		0

We now use this table to make a list of values of c for various temperatures.

clist={{0,0},{0.1,c /. {m->0.9999999,x->0.1}},{0.2,c /. {m->0.9999,x->0.2}},{0.3,c /. {m->0.9974,x->0.3}},{0.4,c /. {m->0.9856,x->0.4}},{0.5,c /. {m->0.9575,x->0.5}},{0.6,c /. {m->0.9073,x->0.6}},{0.7,c /. {m->0.8286,x->0.7}},{0.8,c /. {m->0.7104,x->0.8}},{0.9,c /. {m->0.5253,x->0.9}},{1,1.5}}

Plot the specific heat for the temperatures in the above table.

ListPlot[clist,Joined->False]

Figure 6.7, Figure 6.8
Landau-Ginsburg free energy

Clear[f,m,a,b,h,x,a]

Define the Landau-Ginsburg free energy, f.

f=a (m^2)+b (m^4)-h m

Plot the free energy f for three different values of a.

Plot[{f/.{a->1,b->1,h->0},f/.{a->0,b->1,h->0},f/.{a->-1,b->1,h->0}},{m,-1.5,1.5},PlotRange->{-0.5,1}]

Animate the free energy for a ranging from 1 to –1, with b = 1 and h = 0.

Animate[Plot[{(f/.{a->x,b->1,h->0})},{m,-2,2},PlotRange->{{-1.5,1.5},{-0.5,1.0}},PlotStyle->{Black}],{x,1,-1},AnimationRunning->False]

Now, consider non-zero h. Plot the free energy for three different values of h.

Plot[{f/.{a->-1,b->1,h->0.25},f/.{a->-1,b->1,h->0},f/.{a->-1,b->1,h->-0.25}},{m,-1.5,1.5},PlotRange->{-0.5,1}]

Animate the free energy for h ranging from 0.25 to –0.25, with a = –1 and b = 1.

Animate[Plot[{(f/.{a->-1,b->1,h->x})},{m,-2,2},PlotRange->{{-1.5,1.5},{-0.5,1.0}},PlotStyle->{Black}],{x,0.25,-0.25},AnimationRunning->False]

Chapter 7

Figure 7.3 and Table 7.1
Renormalized coupling for 1-D decimation

Clear[Kp,K,K1,K2,n]

Define the renormalized coupling, K’ and call it Kp.

Kp=(1/2)Log[Cosh[2K]]

Plot Kp as a function of K.

Plot[{Kp,K},{K,0,1.3},PlotStyle->{{Black},{Dashed}},PlotRange->{{0,1.3},{0,1.3}}]

Now, do successive iterations of Kp. You can set K1 to any initial value you like. In this case we have K1 = 1.2, which corresponds to Table 7.1 and the cob-webbing in Figure 7.3.

K1=1.2;
Print["n = 0","   K = ",K1];
Do[
K2=Kp /. K->K1;
Print["n = ",n,"   K = ",K2];
If[Abs[K1-K2] < 0.000001,Break[]];
K1=K2,
{n,1,100}]

Table 7.2
Renormalized coupling and free energy for 1-D decimation

Clear[Kp,K,K1,K2,n,f,fexact]

Define the renormalized coupling, K’ and call it Kp.

Kp=(1/2)Log[Cosh[2K]]

Now, do successive iterations of Kp and sum up the contributions to the free energy. You can set K1 to any initial value you like. In this case we have K1 = 1.2, which corresponds to Table 7.2.

K1=1.2;
f=0;
Do[
K2=Kp /. K->K1;
f = f+((1/(2^(n+1))))(Log[2]+K2);
Print["n = ",n,"        K = ",K1,"                      f = ",f];
If[Abs[((1/(2^(n+1))))(Log[2]+K2)] < 0.000001,Break[]];
K1=K2,
{n,0,50}]

For comparison, here’s the exact free energy for this system.

fexact=Log[2]+Log[Cosh[K]]/.K->1.2

Table 7.3
Renormalized couplings and free energy for 1-D decimation

Clear[Kp,K,K1,K2,n,f,fexact,zpp,zpm,zmp,zmm,hp,K0p,h1,h2]

Define the partial partition functions.

zpp=(E^(2K+2h))+(E^(-2K))

zpm=(E^h)+(E^(-h))

zmp=zpm

zmm=(E^(2K-2h))+(E^(-2K))

Define the renormalized couplings and the renormalized zero level.

Kp=(1/4)Log[((zpp zmm)/(zpm zmp))]

hp=(1/2)Log[(zpp/zmm)]

K0p=(1/4)Log[zpp zpm zmp zmm]

Now, do successive iterations of Kp and hp, and sum up the contributions to the free energy. You can set K1 and h1 equal to any initial values you like. In this case we have K1 = 1.2 and h = 1.2, which corresponds to Table 7.3.

K1=1.2;
h1=1.2;
f=0;
Do[
K2=Kp /.{ K->K1,h->h1};
h2=hp/.{ K->K1,h->h1};
f = f+((1/(2^(n+1))))(K0p/.{ K->K1,h->h1});
Print["n = ",n,"   K = ",K1,"         h = ",h1,"             f = ",f];
If[Abs[((1/(2^(n+1))))(Log[2]+K2)] < 0.000001,Break[]];
K1=K2;
h1=h2,
{n,0,50}]

For comparison, here’s the exact free energy for this system.

fexact=Log[(E^K)Cosh[h]+(Sqrt[(E^(-2K))+(E^(2K))(((Sinh[h]))^2)])]/.{K->1.2, h->1.2}

Figure 7.8
Renormalized coupling for MK transformation

Clear[Kp,K]

Define the renormalized coupling, K’ and call it Kp.

Kp=(1/2)Log[Cosh[4K]]

Plot Kp as a function of K. The fixed point is at approximately K = 0.304....

Plot[{Kp,K},{K,0,0.61},PlotStyle->{{Black},{Dashed}},PlotRange->{0,0.61}]

Table 7.4
Renormalized couplings for MK transformation

Clear[K,Kp,K1,K2]

Define the renormalized coupling, Kp.

Kp=(1/2)Log[Cosh[4K]]

Now, do successive iterations of Kp. Below we show the initial conditions K1 = 0.25 and K1 = 0.35, corresponding to Table 7.4.

K1=0.25;
Do[
K2=Kp /. K->K1;
Print["n = ",n,"   K = ",K1];
K1=K2,
{n,1,10}]

K1=0.35;
Do[
K2=Kp /. K->K1;
Print["n = ",n,"   K = ",K1];
If[Abs[K1-K2] < 0.000001,Break[]];
K1=K2,
{n,1,10}]

Figure 7.9 and Table 7.5
Renormalized coupling for MK transformation and the reversed flow

Clear[K,Kp,K1,K2,a,Krf]

Define the renormalized coupling, K’ and call it Kp.

Kp=(1/2)Log[Cosh[4K]]

Plot Kp as a function of K. The fixed point is at approximately K = 0.304....

Plot[{Kp,K},{K,0,0.61},PlotStyle->{{Black},{Dashed}},PlotRange->{0,0.61}]

Define the reversed flow.

Krf=K-a(Kp-K)

Plot Kp and Krf as a function of K. The fixed point is at approximately K = 0.304....

Plot[{Kp,K,Krf/.a->1},{K,0,0.61},PlotStyle->{{Black},{Dashed}},PlotRange->{0,0.61}]

Determine the fixed point with the reversed flow transformation for two different initial conditions, K1 = 0.1 and K1 = 0.5, corresponding to Table 7.5.

K1=0.1;
Do[
K2=Krf/.{a->1,K->K1};
Print["n = ",n,"   Krf = ",K1];
If[Abs[K1-K2] < 0.000001,Break[]];
K1=K2,
{n,1,20}]

K1=0.5;
Do[
K2=Krf/.{a->1,K->K1};
Print["n = ",n,"   Krf = ",K1];
If[Abs[K1-K2] < 0.000001,Break[]];
K1=K2,
{n,1,20}]

Table 7.6
Renormalized coupling and free energy for MK transformation

Clear[Kp,K,K1,K2,n,f]

Define the renormalized coupling, K’ and call it Kp.

Kp=(1/2)Log[Cosh[4K]]

Now, do successive iterations of Kp and sum up the contributions to the free energy. You can set K1 to any initial value you like. In this case we have K1 = 1.0, which corresponds to Table 7.6.

K1=1.0;
f=Log[2.0];
Print["n = ",n,"   K = ",K1,"       f = ",f];
Do[
K2=Kp /. K->K1;
f = f+((1/2))((1/(4^(n-1))))(K2);
Print["n = ",n,"   K = ",K2,"       f = ",f];
If[Abs[((1/2))((1/(4^(n-1))))(K2)] < 0.000001,Break[]];
K1=K2,
{n,1,100}]

Figure 7.12
Bond-strengthening factor

Clear[y,K,Ktil]

Define an intermediate expression.

y=(1/2)(1+Cosh[8K])

Now, define K tilde, and call it Ktil.

Ktil=(1/4)Log[y+(Sqrt[(y^2)-1])]

Plot the bond-strengthening factor Ktil/K as a function of K.

Plot[{Ktil/K,2,(Sqrt[2])},{K,0,2.1},PlotStyle->{{Black},{Dashed},{Dashed}},PlotRange->{1.35,2.05}]

Figure 7.18
Cell-cluster transformation

Clear[Hcc,K,P1,P2,Zppcc,Zpmcc,sp1,sp2,Kpcc,K0pcc,flist,Kpcc0,fe1,fe2,j,pl1,pl2,u,g,f,dv,th,plfreeenergy]

Define the Hamiltonian for the cell-cluster (cc) transformation, as presented in Figure 7.17.

Hcc=K(s1 (s8 + s2)+s2 (s7 + s1)+s3(s2+s4)+s4(s1+s3)+s5(s4+s6)+s6 (s3+s5)+s7(s6+s8)+s8(s5+s7))

Define the projection for cell 1, call it P1, and the projection for cell 2, call it P2.

P1=((1/2)(1+Sign[(s1+s2+s3+s4)] sp1))

P2=((1/2)(1+Sign[(s5+s6+s7+s8)] sp2))

Define the partial partition function the cell-cluster (cc) transformation with both cell spins equal to +1 (pp). Call it Zppcc.

Zppcc=Sum[P1 P2 (E^Hcc),{s1,-1,1,2},{s2,-1,1,2},{s3,-1,1,2},{s4,-1,1,2},{s5,-1,1,2},{s6,-1,1,2},{s7,-1,1,2},{s8,-1,1,2}] /. {sp1->1,sp2->1}

Define the partial partition function for one cell spin equal to +1, and one cell spin equal to –1. Call it Zpmcc.

Zpmcc=Sum[P1 P2 (E^Hcc),{s1,-1,1,2},{s2,-1,1,2},{s3,-1,1,2},{s4,-1,1,2},{s5,-1,1,2},{s6,-1,1,2},{s7,-1,1,2},{s8,-1,1,2}] /. {sp1->1,sp2->-1}

Define the renormalized coupling, Kpcc.

Kpcc=(1/4)Log[(Zppcc/Zpmcc)]

Plot the renormalized coupling, Kpcc,  as a function of the original coupling, K. The fixed point is at K = 0.391....

Plot[{Kpcc,K},{K,0,0.8},PlotRange->{{0,0.8},{0,0.8}},PlotStyle->{{Black},{Dashed}}]

Define the contribution to the free energy.

K0pcc=(1/2)Log [Zpmcc]

Now calculate the free energy. In this case, we calculate it for 100 points from K = 0 to K = 1.0.

flist={{0,Log[2.]}};
Do[
Kpcc0=(0.01)j;
fe1=0;
Do[K1=Kpcc /. K->Kpcc0;
fe2=fe1+(1/(4^i)) K0pcc /. K->Kpcc0;
If[Abs[fe2-fe1]<0.000001,Break[]];
Kpcc0=K1;
fe1=fe2,
{i,150}];
flist = Append[flist,{(0.01)j,fe2}],
{j,100}]
pl1=ListPlot[flist,Joined->True]

Plot limiting values.

pl2=Plot[{Log[2],2K},{K,0,1.5},PlotStyle->{Dashed}]

Combine the free energy with its limiting values.

Show[pl1,pl2]

To compare with the exact free energy, we repeat the calculation given for Figure 5.2.

To start, define an important variable in the calculation.

u=(1/(((Sinh[2K]))^2))

Define the integrand and call it g. The integration variable is th.

g=Log[((Cosh[2K])^2)+(1/u)(Sqrt[1+(u^2)-2u Cos[2th]])]

Define the free energy f in terms of the integrand g as follows:

f = (1/2)Log[2] + (1/(2π))(∫_0^π)gdth

We’ll put this expression into the Plot command. We’ll also plot the K = 0 limit, Log[2], and the K –> ∞ limit, 2K.

plfreeenergy=Plot[{(((Log[2])/2)+((1/(2π)))NIntegrate[g,{th,0,π}]),2K,Log[2]},{K,0,1},PlotStyle->{Black,Dashed,Dashed},AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["f",24,Plain,Black,FontFamily->"Times New Roman"]}]

Add a vertical line to indicate the location of the critical point.

First, define the location of the critical point, (Subscript[K, c]).

Kc=0.5Log[1+(Sqrt[2])]

pldashed=Graphics[{Dashed,Line[{{Kc,0},{Kc,3}}]}];

Now, combine the graphics.

Show[plfreeenergy,pldashed,AxesStyle->Directive[Black,20],AxesLabel->{Style["K",24,Plain,Black,FontFamily->"Times New Roman"],Style["f",24,Plain,Black,FontFamily->"Times New Roman"]}]

Combine with the cell-cluster free energy from pl1. The two results are almost identical, with the cell-cluster free energy slightly lower than the exact result near the critical point.

Show[plfreeenergy,pl1,pldashed]