The Ising Model Part 2:
Magnetism and Hysteresis

M. E. Mauel
AP1601

Introduction

Last week, we introduced the Ising model, wrote a Monte Carlo simulation, and demonstrated some of basic properties of an interacting two-dimensional lattice. Fundementally, the Ising model is a two-dimensional lattice of digital sites. The "spin" state of each site can be either "up" (+1) or "down" (-1). The sites interact only with its nearest neighbors. If the potential energy is lowered when neighboring spins align, the lattice is a simple model of a ferromagnet. If the potential energy is lowered when neighboring spins are antiparallel, then the lattice can represent antiferromagnetic material. The dynamics and statistical properties of the lattice are determined entirely by the exchange or interaction energy matrix.

This week, we extend the model to include the interaction energy between each lattice site and an applied magnetic field. We'll focus on paramagnetic (and ferromagnetic) materials. The energy will be lowest when the spins align with the magnetic field.

The total energy of our lattice is given by the sum:
	E = - (∑_(i,j))(Subscript[J, ij]) (Subscript[s, i]) (Subscript[s, j]) – B (∑_i)(Subscript[s, i]).
Recall that the interaction matrix, (Subscript[J, ij]), couples nearest neighbors in the simple Ising model and is zero for all {i, j} which are not nearest neighbors. The only possible values of the interaction energy are {+4, +2, 0, -2, -4} corresponding to the five cases {all aligned, three aligned, two aligned and two not aligned, three not aligned, and all not aligned}. The indicies, i and j, range across every site. (They are not x- and y- coordinate indicies!) The magnetic field adds a potential energy proportional to the direction of the spin and the strength of the magnetic field, B.

After describing our Ising program with an applied magnetic field, we'll investigate some of its properties and relate these to those of magnets.

A Mathematica Ising model you can use...

We use the Module[ ] structure to collect a series of Mathematica expressions into a convenient function that can be applied over and over again. This function performs m Monte Carlo steps. Each step consists on average of flip each site once. A long list of lattice configurations are returned as a result of evaluating the expression.
(Note: the meaning of the second argument, m, has changed from the function Ising1[...] discussed last week. Below, it refers to the total number of lattice configurations returned as a result of the simulation.)

Ising2[n_,m_,p_,T_,B_,startGrid_:0] :=
Module[{energyDiff,initConfig,flip,flipList},

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

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

initConfig=If[MatchQ[startGrid,_List],startGrid,(* else make a random starting grid *)Table[2 Floor[p+RandomReal[]]-1,{n},{n}]];

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 (n^2)];
flipList[[Range[1,m (n^2),(n^2)]]]
]

Magnetization

The magnetization function returns the average spin direction of all of the sites.
	M = ⟨ (Subscript[s, i])⟩ = (1/N) (∑_(i=1)^N)(Subscript[s, i])

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

Energy

The energy function returns the average energy of all of the sites. It includes the potential energy of spin alignment to the magnetic field, B.
	E = - (∑_(i=1)^N)(Subscript[s, i])((∑_(j=1)^4)(Subscript[s, j]) + B), where the sum over j represents the four neighbors.

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

Test 1: Adding magnetic field above the Curie temperature

Remember, Onsager showed that the critical temperature of the Ising model is T = 2.269. We call this the "Curie point" in recognition of the French physicist Pierre Curie, who in 1895 discovered the laws that relate ferromagnetic properties to temperature. Above the Curie point, a material is paramagnetic; whereas, below the Curie temperature, ferromagnetics can be permanently magnetized.

nSize = 20; (* 20 x 20 grid requires about 20 seconds to compute thermal equilibria *)

Positive Magnetic Field, B > 0

t1 = Ising2[nSize,100,0.5,4.5, 1.0]; // Timing

t1 is now set to a list of 100 two-dimensional lists...

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

mPos=ListPlot[magnetization/@t1,PlotRange->{-1.,1.},PlotStyle->{RGBColor[1,0,0]}]

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

Negative Magnetic Field, B < 0

t1 = Ising2[nSize,100,0.5,4.5, -1.0];

Show[mPos,ListPlot[magnetization/@t1,PlotRange->{-1.,1.},PlotStyle->{RGBColor[0,0,1]},DisplayFunction->Identity]]

Of course, the average energy is independent of the direction of the applied magnetic field.

Show[ePos,ListPlot[(energy[#1,-1.]&)/@t1,PlotRange->{0,-4},PlotStyle->{RGBColor[0,0,1]},DisplayFunction->Identity]]

Paramagnetism

magneticfields = Table[B, {B, -10, 10, 1}]

The following takes a very looonnnggggg time (about 1 minute on my laptop.) We perform a Monte Carlo simulation for each value of the magnetic field strength. We save the last 50 lattice configurations. We'll average the magnetization of these last 50 and plot the average magnetization as a function of applied field.

(iSims = ( Ising2[nSize,100,0.5, 4.5, #][[Range[51,100]]] )&
/@ magneticfields;) // Timing

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

ListPlot[Transpose[{magneticfields,avgMagnet}],PlotStyle->{PointSize[0.02]},PlotRange->{-1,1},PlotLabel->"Magnetization vs Applied B (T = 4.5)"]

This result will change with temperature. "Mean-field" theoretical approximations produce a prediction, magnetization  = Tanh[B / T]. You may be interested in testing this approximation as part of your project.

Test 2: Adding Field below the Curie temperature

When T < 2.269, the Ising lattice represents a ferromagnet. Ferromagnets exhibit hysteresis. The following program keeps a record of the magnetization as B slowing increases or decreases.

magneticfields = Table[B, {B, -1.5, 1.5, 0.1}]

Increasing magnetic field

iSims={};
lastLattice=Table[2 Floor[0.5+RandomReal[]]-1,{nSize},{nSize}];
Do[iSims=Append[iSims,Ising2[nSize,100,0.5,1.5,magneticfields[[i]],lastLattice][[Range[51,100]]]];lastLattice=Last[Last[iSims]],{i,1,Length[magneticfields]}];

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

mIncreasing=ListPlot[Transpose[{magneticfields,avgMagnet}],PlotStyle->{RGBColor[1,0,0]},PlotRange->{-1,1},Joined->True]

Decreasing magnetic field

magneticfields = Reverse[magneticfields]

iSims = {};
Do[iSims = Append[iSims , Ising2[nSize,100,0.5, 1.5, magneticfields[[i]],lastLattice][[Range[51,100]]] ];
lastLattice = Last[Last[iSims]],{i, 1, Length[magneticfields]}
];

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

Show[mIncreasing,ListPlot[Transpose[{magneticfields,avgMagnet}],PlotStyle->{RGBColor[0,0,1]},PlotRange->{-1,1},Joined->True,DisplayFunction->Identity],PlotLabel->"Magnetization vs Applied B (T = 1.5)"]

avgEnergy = (Plus @@@ Table[(energy[#, magneticfields[[i]]]& /@iSims[[i]]), {i, 1, Length[magneticfields]}])/50 // N;

ListPlot[Transpose[{magneticfields,avgEnergy}],PlotStyle->{PointSize[0.02`],RGBColor[0,0,1]},PlotRange->{-6,0},PlotLabel->"Energy vs B (as B Decreases)"]

Explanations and other measurements...

Phase-transistions are among the most difficult phenomona to derive from first principles.  At phase-transitions, matter becomes highly correlated. A long-range order exists extending across the entire lattice. These create significant fluctuations in thermodynamic (average) properties of matter at phase transitions.
The fluctuations of a system in thermodynamic equilibrium are related to the response of a system to applied forcings. For example, the specific heat, C, describes how the energy changes due to changes in temperature. This is also proportional to the average-squared energy fluctuations, ⟨((ΔE)^2)⟩. The average-squared fluctuations in the magnetic field, ⟨((ΔM)^2)⟩, are proportional to the magnetization.

Correlations

nSize = 40; (* a bigger lattice for improved statistics *)

correlation[lat_List, direction_] := 
(Plus @@ Flatten[lat RotateRight[lat,direction]])/Length[lat]^2 // N

t1 = Ising2[nSize,100,0.5,2.0, 0.0];

MatrixPlot[Last[t1],PlotLabel->"2.0 Degrees"]

correlation[Last[t1],#]& /@ {{1,0}, {2,0}, {3,0}, {4,0}}

Table[{k,0},{k,1,(nSize/2)}]

p1=ListPlot[Evaluate[(correlation[Last[t1],#1]&)/@Table[{k,0},{k,1,(nSize/2)}]],PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]

t2 = Ising2[nSize,100,0.5,3.0, 0.0];

MatrixPlot[Last[t2],PlotLabel->"3.0 Degrees"]

correlation[Last[t2],#]& /@ {{1,0}, {2,0}, {3,0}, {4,0}}

p2=ListPlot[Evaluate[(correlation[Last[t2],#1]&)/@Table[{k,0},{k,1,(nSize/2)}]],PlotStyle->{RGBColor[0,1,0],PointSize[0.02]}]

t3 = Ising2[nSize,100,0.5,5.0, 0.0];

MatrixPlot[Last[t3],PlotLabel->"5.0 Degrees"]

correlation[Last[t3],#]& /@ {{1,0}, {2,0}, {3,0}, {4,0}}

p3=ListPlot[Evaluate[(correlation[Last[t3],#1]&)/@Table[{k,0},{k,1,(nSize/2)}]],PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}]

Show[p1,p2,p3,PlotLabel->"Spin Correlation for T = 2.0, 3.0, 5.0 degrees",PlotRange->{0.,1.}]

Energy Fluctuations

The specific heat is related to the total energy of system by the relationship
	C = (∂_T) ⟨E⟩. 
It represents how much energy (i.e. "heat capacity") can be stored for a given temperature increment. According to the "fluctuation-dissipation theorem", the specific heat is also related to the thermal fluctuations of the stored energy as
	C = (((⟨ΔE⟩)^2)/(T^2)).

The fluctuations in the energy of the lattice are proportional to the specific heat, C = (((ΔE))^2)/ (T^2).

ListPlot[energy[#,0.0]& /@ t1]

energy[#,0.0]& /@ Take[t1,-50] (* only take the lattice state AFTER thermal equilibrium has been reached *)

nSize^2 Variance[energy[#,0.0]& /@ Take[t1,-50]]/(2.0)^2

nSize^2 Variance[energy[#,0.0]& /@ Take[t2,-50]]/(3.0)^2

nSize^2 Variance[energy[#,0.0]& /@ Take[t3,-50]]/(5.0)^2

Magnetization Fluctuations

For magnetic materials, the susceptibility is similarly related to the magnetization
	χ =  (∂_H) ⟨M⟩ =  (((⟨ΔM⟩)^2)/T).

The fluctuations in the magnetization of the lattice are proportional to the magnetic susceptibility, M = (((ΔM))^2)/ T.

ListPlot[magnetization /@ t1]

magnetization /@ Take[t1,-25]  (* only after equilibrium: last 25? *)

nSize^2 Variance[magnetization /@ Take[t1,-25]] /(2.0)

nSize^2 Variance[magnetization /@ Take[t2,-25]] /(3.0)

nSize^2 Variance[magnetization /@ Take[t3,-25]] /(5.0)

Summary

We have examined the effects of an applied magnetic field to our Ising model for a magnetic material. We've observed hystersis, and the abrupt transitions between up and down magnetization at low temperatures.

We have begun to look at statistical or average properties of the lattice. These are used to derive thermodynamic properties of the system.