Statistical Physics and 
Phase Transitions:
The Ising Model

M. E. Mauel
AP1601

Introduction

This notebook presents a short discussion of the collective interaction of two-dimensional systems generally related to Ising models. Ising models are used to study critical behavior of many complex, collective systems. (A good reference for this part of the course is the textbook Computational Physics by N. J. Giordano, Chapters 8 and 11.)

The Ising model was first introduced in 1925 as a model for ferromagnetism. It consists of a regular array in which each site can be occupied in one of two ways: up (+ 1) or down (- 1). When referring to magnetic materials, the direction of each site corresponds to the orientation of an atomic magnetic moment. We will use the symbol (Subscript[s, i]) = +- 1 to refer to the spin orientation. Each site interacts with its neighbors, and the model begins with a formula used to calculate the total potential energy, E, of the entire collection of sites based on the relative orientation of the sites. 
In general this energy will be equal to the sum of the individual interaction energies between all pairs of sites. Let's label two sites with the indices i and j. We will use the symbol (Subscript[J, ij]) to represent the interaction between the two sites i and j. (Subscript[J, ij]) is called the "exchange constant" and it will be a number. The contribution to the potential energy due to the interaction between these two sites will be (Subscript[E, ij]) = - (Subscript[J, ij])(Subscript[s, i])(Subscript[s, j]). (The negative sign is chosen for convention.) With (Subscript[J, ij]) > 0, the potential energy is lowered if the spins are parallel. This corresponds to ferromagnetism. If (Subscript[J, ij]) < 0, the potential energy is lowered if the spins are antiparallel.

The total potential energy of a square lattice can be written as

energy = - Sum[ j[ {i1,i2}, {j1,j2}]  s[[i1, i2]]  s[[j1, j2]], 
		{i1, 1, nSize}, {i2, 1, nSize}, 
		{j1, 1, nSize}, {j2, 1, nSize}]
or
	E = - (∑_ij)(Subscript[J, ij])(Subscript[s, i])(Subscript[s, j]) .

For an nSize x nSize grid, this is a large sum of nSize^4 terms! Fortunately, the Ising model for magnetic systems considers the interaction energy to be significant only between neighboring sites. For a square grid, there are four neighbors located to the right, left, up and down, and four more neighbors located on  diagonal sites.

Lars Onsager

The Norwegian-born American physical chemist Lars Onsager, b. Nov. 27, 1903, d. Oct. 5, 1976, won the 1968 Nobel Prize for chemistry for his theoretical studies on the thermodynamics of irreversible processes. Building on the work of others, Onsager laid the foundation for a comprehensive theory of the movement of ions in electrified solutions. He later became interested in low-temperature physics, and he predicted the existence of vortices in superfluid helium.

Lars Onsager was also the first to present an exact solution to the two-dimensional Ising model of ferromagnetism for a square lattice in the absence of an external magnetic field. Onsager's work with the Ising model constituted one of the few exact evaluations of a partition function for systems in which inter-particle interactions cannot be neglected. Onsager's complete solution is somewhat difficult to discuss at an introductory level, but we can reproduce many of his results (and the results of others) using numerical simulation of statistical processes. For example, we can determine (1) the critical Temperature below which the system is ferromagnetic (or antiferromagnetic), (2) calculate the specific heat as a function of temperature and examine the differences between a so-called first-order and second-order phase transition (ice melting is a first-order phase-transition; whereas, the Curie-point of a ferromagnet is a second-order phase-transition), (3) examine the hysteresis of magnets, and (5) show the relationship between collective or correlated activity of the sites and phase transitions. 
As a first step, we will reproduce Onsager's result that the critical temperature of a two-dimensional ferromagnetic Ising model is given by (Subscript[T, c]) = 2.269 J /k.

Calculating the Potential Energy of a Configuration

We need to calculate the potential energy of a configuration of "spin sites."

Let's start with a random lattice configuration. We use the symbol p to represent the fraction (on average) of sites which are up.

Random Numbers

ranLattice[nSize_Integer][p_]:=Table[2 Floor[p+RandomReal[]]-1,{nSize},{nSize}]

(ranArray = ranLattice[8][0.5]) // MatrixForm

Count[Flatten[ranArray], _?Positive]/Count[Flatten[ranArray], _]

A Random Lattice

nSize = 25;

initial=ranLattice[nSize][0.75];
MatrixPlot[initial]

(Notice that "white" represents up  and "black" represents down.)

In order to calculate the energy of this system, we need to sum over every site the contribution to the potential energy of that site with respect to each of its neighbors. 
For simplicity, we will look only at the four neighbors to the right, left, up, and down on the square lattice. We will further assume periodic boundary conditions. The top of the lattice is connected to the bottom and the left is connected to the right.

Mathematica makes it relatively easy to sum lists of numbers, and it has built-in functions for shifting the lists up, down, right, and left. For example,

Show[GraphicsRow[{MatrixPlot[initial,PlotLabel->"Original",DisplayFunction->Identity],MatrixPlot[RotateRight[initial,{1,0}],PlotLabel->"Shift Up One",DisplayFunction->Identity],MatrixPlot[RotateRight[initial,{0,1}],PlotLabel->"Shift Right One",DisplayFunction->Identity]}]]

We can add together the spin values of the four neighbors...

neighbors = Plus @@ (RotateRight[initial, #]& /@ {{1,0}, {-1, 0}, {0,1}, {0,-1}});

MatrixPlot[neighbors]

Assuming J  = 1, the energy contribution of every site is...

MatrixPlot[-initial neighbors]

Where "white" is positive energy and "black" is negative energy.

The total energy of the system is

- (Plus @@ Flatten[initial neighbors]) /nSize^2// N

If we make the spins "mostly down" (instead of "mostly up"), the energy is similar...

initial = ranLattice[nSize][0.25];
neighbors = Plus @@ (RotateRight[initial, #]& /@ {{1,0}, {-1, 0}, {0,1}, {0,-1}});- (Plus @@ Flatten[initial neighbors]) /nSize^2// N

But, if the spins are not ordered in any which way, ...

initial = ranLattice[nSize][0.5];
neighbors = Plus @@ (RotateRight[initial, #]& /@ {{1,0}, {-1, 0}, {0,1}, {0,-1}});- (Plus @@ Flatten[initial neighbors]) /nSize^2// N

The Ising system has the lowest potential energy when the spins are aligned either up or down!

energy[p_] := Module[{initial, neighbors},
initial = ranLattice[nSize][p];
neighbors = Plus @@ (RotateRight[initial, #]& /@ {{1,0}, {-1, 0}, {0,1}, {0,-1}});- (Plus @@ Flatten[initial neighbors]) /nSize^2// N]

ListPlot[Table[{p,energy[p]},{p,0.,1.,0.025}],PlotStyle->PointSize[0.02],Frame->True,PlotLabel->"Energy vs. Orientation",FrameLabel->{"Up fraction","E"}]

Monte Carlo Simulation (Rosenbluth, Teller, and Metropolis, 1953)

In this section, the famous "Rosenbluth, Teller, and Metropolis" algorithm (Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller, "Equation of State Calculations by Fast Computing Machines," J. Chem. Phys. 21 (1953) pp. 1087-1092) will be used to simulate the thermal fluctuations of a 2D magnetic lattice.

The random thermal agitation of the magnetic moments of a large collection interacting particles is "simulated" as previously used for interacting particles. Each particle can become "agitated" and "flip" the orientation of its spin. Our simulation will start at an initial state and continue through a series of time steps towards some sort of "thermal equilibrium" or a "relaxed state". 
How is this done?

We will follow a set of specific rules.

Calculating the Energy Difference

Show[Graphics[{Arrow[{{0,0},{0,1}}],Hue[0],Arrow[{{1,0},{1,1}}],Arrow[{{-1,0},{-1,1}}],Arrow[{{0,1},{0,2}}],Arrow[{{0,-1},{0,0}}]}],PlotLabel->"Low Energy = - 4"]

Show[Graphics[{Arrow[{{0,1},{0,0}}],Hue[0],Arrow[{{1,0},{1,1}}],Arrow[{{-1,0},{-1,1}}],Arrow[{{0,1},{0,2}}],Arrow[{{0,-1},{0,0}}]}],PlotLabel->"High Energy = + 4"]

What are all of the possible combinations?
If we list all possible cases, the first argument is the spin state of the site and the second argument is the total spin state of the neighbors.

energyDiff[ 1, 4] =  2 (4 );
energyDiff[ 1, 2] =  2 (2 );
energyDiff[ 1, 0] =  2 ( 0 );
energyDiff[ 1,-2] = 2 (- 2);
energyDiff[ 1,-4] = 2 (- 4);
energyDiff[-1,4] = -2 (4);
energyDiff[-1,2] = -2 (2);
energyDiff[-1,0] = -2 (0);
energyDiff[-1,-2] = -2 (-2);
energyDiff[-1,-4]= -2 (-4);

This function can be used to define a "flip" function...

flip[lattice_List,temp_] := Module[{lat,i1,i2,rt,lt,up,dn,n,neighbors},
lat=lattice;
n=Length[lat];

{i1,i2}=RandomInteger[{1,n},2];

If[i1==n,dn=1,dn=i1+1];

If[i2==n,rt=1,rt=i2+1];
If[i1==1,up=n,up=i1-1];
If[i2==1,lt=n,lt=i2-1];
neighbors=lat[[dn,i2]]+lat[[up,i2]]+lat[[i1,rt]]+lat[[i1,lt]];

If[energyDiff[lat[[i1,i2]],neighbors]<0||RandomReal[]<Exp[-((energyDiff[lat[[i1,i2]],neighbors])/temp)],lat[[i1,i2]]=-lat[[i1,i2]];lat,lat]
]

MatrixPlot[flip[initial,1.]]

MatrixPlot[Nest[flip[#1,1.]&,initial,10 (nSize^2)]]

MatrixPlot[Nest[flip[#1,4.]&,initial,10 (nSize^2)]]

Temperature and Phase-Transitions

Using the Module[ ] structure, we can write a convenient program...

Ising1[n_,m_,T_,p_] :=
Module[{initConfig,flip,flipList},

initConfig=ranLattice[n][p];

flip=(lat=#1;{i1,i2}=RandomInteger[{1,n},2];If[i1==n,dn=1,dn=i1+1];If[i2==n,rt=1,rt=i2+1];If[i1==1,up=n,up=i1-1];If[i2==1,lt=n,lt=i2-1];nnvalsum=lat[[dn,i2]]+lat[[up,i2]]+lat[[i1,rt]]+lat[[i1,lt]];If[energyDiff[lat[[i1,i2]],nnvalsum]<0||RandomReal[]<Exp[-((energyDiff[lat[[i1,i2]],nnvalsum])/T)],lat[[i1,i2]]=-lat[[i1,i2]];lat,lat])&;


flipList=NestList[flip,initConfig,m];
flipList[[Range[1,m,(n^2)]]]
]

Magnetization

magnetization[lat_List] := 
(Plus @@ Flatten[lat])/Length[lat]^2

Energy

energy[lat_List] := 
Module[{neighbors},
neighbors =Plus @@ (RotateRight[lat, #]& /@ {{1,0}, {-1, 0}, {0,1}, {0,-1}});
-(Plus @@ Flatten[lat neighbors])/Length[lat]^2]

Test 1: Cold Evolution (T = 1.0)

nSize = 25;

t1 = Ising1[nSize,100 nSize^2,1.0,0.5];

t1 // Shallow[#,{4, 1}]&

ListPlot[magnetization/@t1,PlotRange->{-1.`,1.`}]

ListPlot[(energy[#1]&)/@t1,PlotRange->{0,-4}]

Animation

ListAnimate[Table[MatrixPlot[t1[[n,All,All]], PlotRange->{-1,1}], {n, 1, 100}], AnimationRunning -> False]

Test 2: Hot Evolution (T = 4.0)

t2 = Ising1[nSize,100 nSize^2,4.0,0.9];

ListPlot[magnetization/@t2,PlotRange->{-1.`,1.`}]

ListPlot[(energy[#1]&)/@t2,PlotRange->{0,-4}]

Animation

ListAnimate[Table[MatrixPlot[t2[[n,All,All]], PlotRange->{-1,1}], {n, 1, 100}], AnimationRunning -> False]

The Critical Temperature

We calculate a critical temperature at which permanent magnetization is not possible.

nSize = 25;

temperatures = Table[t, {t, 1.5, 4.0, 0.25}];

iSims = ( Ising1[nSize,100 nSize^2,#,0.75][[Range[51,100]]] )&/@ temperatures;

iSims // Shallow[#, 3]&

avgMagnet = Plus @@@ Map[magnetization, iSims, {2}]/50 // N

avgEnergy = Plus @@@ Map[energy, iSims, {2}]/50 // N

ListPlot[Transpose[{temperatures,avgMagnet}],PlotStyle->{PointSize[0.02`]},PlotRange->{-1,1},PlotLabel->"Magnetization vs Temperature",AxesOrigin->{1.5`,0.`}]

Recall, Onsager gave the critical temperature at 2.269!

ListPlot[Transpose[{temperatures,avgEnergy}],PlotStyle->{PointSize[0.02`]},PlotRange->{-4,0},PlotLabel->"Energy vs Temperature",AxesOrigin->{1.5`,0.`}]

Summary and Other Numerical Experiments

We have written a "Monte Carlo" simulation of a two-dimensional lattice representing a ferromagnet. This so-called "Ising Model" has a critical temperature of 2.269 degrees. We represent temperature by the degree that positive changes in the lattice energy are allowed. Higher temperatures allow for greater energy fluctuations. At low temperature, the lattice is mostly aligned—either entirely "up" or entirely "down". As the temperature slowly increases, thermal fluctuations destroy the long-range order, or alignment, of the lattice. (For a real ferromagnet, the critical temperature is called the "Curie Temperature". Ferromagnets can be permanently magnetized only below the Curie temperature.)

The Ising model can be used for many additional experiments. We'll add the effects of an external magnetic field and explore hysteresis and first-order phase transitions. We'll display and try to understand the evolution of the collective structures or correlations within an Ising lattice and the relationship of the fluctuations to the heat capacity.
The Ising model can also be used to analyze lattices with very long-range interactions. It is related  to the use of neural networks to the problem of two-dimensional pattern recognition.