\documentstyle[med_headings]{article}
\oddsidemargin=0in\textwidth=6in\topmargin=-0.80in\textheight=10in
\renewcommand{\textfraction}{0.01}
\pagestyle{empty}\begin{document}\title{}\author{David J.C. MacKay }\date{}
\noindent
\input{psfig.tex}
\input{/home/mackay/tex/newcommands1}
\begin{center}
{\large \bf
Computational Methods in Physics
\\ \medskip
}
\end{center}
%
\newcommand{\bb}{{\bf b}}
\section*{Ising Models}
An `Ising model' is an array of spins (\eg, atoms that can take
states $\pm 1$) that are
magnetically coupled to each other. In a ferromagnetic model it
is energetically favourable for neighbours to
be in the same state; in antiferromagnets, it is favourable for
them to be in opposite states.
Let the state $\bx$ of an Ising model with $N$ spins be a vector in which
each component $x_n$ takes values $-1$ or $+1$. If two spins $m$ and
$n$ are neighbours we write $(m,n) \in {\cal N}$. The energy of a state
$\bx$ is
\beq
E(\bx;J,H) = - \left[
% \frac{1}{2}
\sum_{(m,n) \in {\cal N}} J x_m x_n
+ \sum_{n} H x_n \right] ,
\eeq
where $J$ is the coupling between spins,
% $m$ and $n$,
and $H$ is
the applied field. If $J > 0$ then the
model is ferromagnetic, and if $J < 0$ it is antiferromagnetic.
(Note that each pair $(m,n)$ is counted just once in the first sum).
\begin{center}
\setlength{\unitlength}{2pt}
\begin{picture}(100,45)(0,-5)
\newsavebox{\verticalfour}
\savebox{\verticalfour}(0,0)[bl]{
\multiput(0,0)(0,10){4}{\circle{2}} % spins
\multiput(0,5)(0,10){4}{\line(0,-1){3}} % lines down
\multiput(0,-5)(0,10){4}{\line(0,1){3}} % lines up
\multiput(2,0)(0,10){4}{\line(1,0){6}} % lines right
}
\multiput(0,0)(10,0){12}{\usebox{\verticalfour}}
\put(-10,15){\makebox[0in]{$W$}}
\put(-10,20){\vector(0,1){10}}
\put(-10,10){\vector(0,-1){10}}
\put(50,38){\makebox[0in]{$C$}}
\put(42,40){\vector(-1,0){42}}
\put(58,40){\vector(1,0){52}}
\end{picture}
\end{center}
In this problem we will study a two-dimensional
planar Ising model of finite width $W$ and large length $C$
with periodic boundary conditions
using the `transfer matrix' method.
The aim is to compute the partition function
\beq
Z(\beta,\bJ,H) \equiv \sum_{\bx} \exp
\left[ - \beta E( \bx ; \bJ , H ) \right]
\eeq
where the summation is over all states $\bx$, and the inverse
temperature is $\beta = 1/T$.\footnote{Let $k_B = 1$.}
From the partition function we can obtain many
important physical properties --- for example,
the free energy is given by $F = -
\frac{1}{\beta} \log Z$. The number of states is $2^N$, where
$N=W C$, so direct
computation of the partition function by enumeration of all global states
$\bx$ is not possible for large
$N$.
In the transfer matrix method we rewrite the
partition function in the form
\beqan
Z& = & \sum_{s_1}\sum_{s_2}\cdots \sum_{s_C}
\left[ \prod_{c=1}^{C} M_{s_{c},s_{c+1}} \right] \\
& = & \Trace \left[ \bM^C \right]
\eeqan
where $s_1$ is an integer running over the $2^W$ microstates
of the spins in column 1, $s_2$ runs over the
microstates
of the spins in column 2, and $s_{C+1}$ and $s_1$ are identical, and
$\bM$ is a matrix made up of the Boltzmann factors
\newcommand{\lE}{{\cal E}}
\beq
M_{s s'} = \exp \left( -\b \lE(s,s') \right) ,
\eeq
where $\lE(s,s')$ is an appropriately defined energy
that counts the contribution to the energy made by
two adjacent columns whose microstates are $s$ and $s'$.
Once we have the matrix $\bM$, the partition function can be written
\beqan
Z& = & \Trace \left[ \bM^C \right] \\
& = & \sum_a \mu_a^C ,
\label{eq.Z.prods}
\eeqan
where $\{ \mu_a \}_{a=1}^{2^W}$ are the eigenvalues of $\bM$.
As the length of the strip $C$ increases, $Z$ becomes dominated by the
largest eigenvalue $\mu_{\max}$:
\beq
Z \rightarrow \mu_{\max}^C .
\eeq
So the free energy per spin in the limit of an infinite thin strip is
given by:
\beq
f = - kT \log Z / (WC) = - kT C \log \mu_{\max} / (WC )
= - kT \log \mu_{\max} / W .
\eeq
All we need to do is evaluate the matrix $\bM$ and find its
largest eigenvalue $\mu_{\max}$. There are NAG routines for
finding eigenvalues. One way to find $\mu_{\max}$ is iterative
multiplication of $\bM$ by a vector, starting, for example, from the
vector $(1,1,1,\ldots 1)$.
Evaluate the energies $\lE(s,s')$ for a long thin strip rectangular
Ising model with $J= \pm 1$ and $H=0$, taking care not to count
couplings twice. Evaluate the matrix $\bM$ and its largest eigenvalue,
and compute the free
energy per spin, $f(\beta,J,H)$ for widths ranging from $W = 2$ up to
$W=5$ and $W=6$ and perhaps $W=8$, as
a function of $\beta$ for $H=0$. Examine both the
cases $J = +1$ and $J = -1$ and compare their free energies,
especially at low temperature. Plot the graphs as a function
of $T$ also and think carefully about the results you expect.
Do the results for odd $W$ differ from
the results for even $W$? (Compare $W=5$ with $W=6$.)
Extract the mean energy per spin, the
entropy per spin, and the heat capacity per spin of the Ising models,
as a function of $\beta$,
from your computations. Also plot the variance of the energy
as a function of temperature. Compare the results for $J = +1$ and $J = -1$.
[Hints: $E = - \partial \log Z / \partial \beta$;
$F = - \frac{1}{\b} \log Z$;
$F = E - TS$
$\Rightarrow
S = \log Z + \b \partial \log Z / \partial \beta$;
$C = \partial E / \partial T
= \b^2 \partial^2 \log Z / \partial \beta^2
= \frac{ \partial^2 \log Z / \partial \beta^2 }{T^2 }$;
${\rm var}(E) = \partial^2 \log Z / \partial \beta^2$.
The interesting stuff should happen between $\beta = 0.1$ and $\beta = 10$.
You may find it helpful to work out exactly
the values of $F$, $E$, $S$ and ${\rm var}(E)$ for the simple case of
a two--level system, \ie, a single spin in an applied field.]
\subsection*{Part 2}
Repeat the above computations for a triangular Ising model of width
6 or 8, with $J = \pm 1$.
\begin{center}
\setlength{\unitlength}{2pt}
\begin{picture}(100,45)(0,-5)
\newsavebox{\verticalfourdiag}
\savebox{\verticalfourdiag}(0,0)[bl]{
\multiput(0,0)(0,10){4}{\circle{2}} % spins
\multiput(0,5)(0,10){4}{\line(0,-1){3}} % lines down
\multiput(0,-5)(0,10){4}{\line(0,1){3}} % lines up
\multiput(2,1)(0,10){4}{\line(2,1){6}} % lines rightup
\multiput(2,-1)(0,10){4}{\line(2,-1){6}} % lines rightdown
}
\multiput(0,0)(20,0){6}{\usebox{\verticalfourdiag}}
\multiput(10,5)(20,0){6}{\usebox{\verticalfourdiag}}
\put(-10,15){\makebox[0in]{$W$}}
\put(-10,20){\vector(0,1){10}}
\put(-10,10){\vector(0,-1){10}}
\end{picture}
\end{center}
Do you expect the triangular Ising model to show different physical
properties from the rectangular Ising model?
Look at the free energy as a function of $T = 1/\beta$
at low temperature.
Look at the variance of the energy as a function of $T$.
\subsection*{Optional extra}
Investigate the dependence of the free energy on the applied field $H$
and find the susceptibility.
\medskip
\hfill DJCM \today
\newpage
\section*{Notes for demonstrators}
\subsection*{Overall strategy}
It may be a good idea to make two programs, one which finds the
largest eigenvalue as a function of $\beta$, and another which
reads in the eigenvalues and computes the desired physical quantities.
The numbers should be written to files and the plotting should be
done separately using gnuplot.
\subsection*{Finding the free energy}
The central computation involves an $S \times S$ matrix where $S = 2^W$.
Label the possible states of a column of $W$ spins by integers $s$
between 0 and $2^W-1$.
Each element of the matrix $M_{ss'}$ is evaluated as follows. Count
up the number of positive and negative interactions between the spins
in the first column (not forgetting the wrap around connection from
the top spin to the bottom spin). Call the difference between these
integers $n_s$. The contribution to the energy made by these
interactions is $n_s J$. For example if $W=4$ and $s=9$ (1001 in
binary) then $n_s = 2-2=0$. Do the same for the second column,
obtaining $n_{s'}$. Then count the number of positive and negative
interactions between the spins in the first column and the second
column, $n_{ss'}$. For example between $s=9$ (binary 1001) and $s' = 1$
(binary 0001) there are 3 positive interactions and one negative
interaction. Then let $M_{ss'} = \exp [ \b J ( n_{ss'} + \frac{1}{2} n_{s}
+ \frac{1}{2} n_{s'})]$. The factors of $\frac{1}{2}$ are included so
that when we concatenate several matrices together, obtaining $\bM^N$,
we only count each within-column interaction once.
The largest eigenvalue of $\bM$, $\l_{\max}$,
is found. The log partition function per spin of the infinite thin strip is
$\log(\l_{\max}) / W$.
The free energy per
spin is $f = - \log(\l_{\max}) / ( W \b )$.
\subsection*{Computational ideas:}
When computing the eigenvalues, one is free to define the matrix
in a variety of ways, some of which are not symmetric. Is using a symmetric
matrix a good idea?
Only the largest eigenvalue is needed. There are several ways of getting
this quantity, for example, iterative multiplication of the matrix.
Because the matrix is all positive we know that the principal
eigenvector is all positive too (Frobenius--Perron theorem), so a
reasonable initial vector is $(1,1,\ldots,1)$. This iterative
procedure may be faster than explicit computation of all eigenvalues.
\subsection*{How to think about results}
The students should be encouraged to figure out what
the free energy, energy, entropy, heat capacity and variance of energy
are for the simple two--level system (the Schottky anomaly).
\subsection*{Physics insights:}
For large temperatures all Ising models should show the same
behaviour: the free energy is entropy-dominated, and the entropy per
spin is $\log(2)$. The mean energy per spin goes to zero.
The free energy per spin should tend to
$-\log(2)/\beta$. The variance of the total energy should tend to
a constant equal to the number of couplings.
In the limit of small temperature, the free energy per spin
can be thought of as a sum of two terms: the energy per spin of the
ground state, plus the entropy per spin of the ground state, if it is
degenerate. When $J=1$ and $b=0$, a ferromagnet has a unique
ground state with energy per spin -2.0 (four bonds, each shared between
two spins).
If $W$ is even then the antiferromagnet has an almost unique ground
state (degeneracy 2) with the same energy per spin, -2.0. (Note that at low
temperatures the unfrustrated antiferromagnetic is indistinguishable from
the ferromagnet.) But if $W$
is odd then the antiferromagnet is frustrated in the width direction;
this affects both the
energy per spin, which is not so negative; and also, in principle,
the entropy per spin, because the ground state of the frustrated
system may be significantly degenerate, with a non-zero entropy per
spin. In the case $W = 3$ the free energy per spin is -4/3 instead
of -2. The ground state only has degeneracy 2. For the rectangular
geometry I think that for any $W$ the ground state has
finite degeneracy. As $W$ increases this effect becomes negligible.
The degeneracy of the triangular system at low temperature is
different. The ground state is extensively
degenerate, at least for all even values of $W$.
(It is instructive to create ground states on the back of an
envelope.)
% \footnote{By constructing ground states for $W=3$
% on the back of an
%
% : f: -1.0807 lz: +2.3283 T +0.46416 log(beta) +0.76753:
% : f: -1.0482 lz: +3.7671 T +0.27826 log(beta) +1.2792:
% : f: -1.0289 lz: +6.1681 T +0.16681 log(beta) +1.7909:
% : f: -1.0173 lz: +10.1733 T +0.1 log(beta) +2.3026:
%
% using T = 0.464, obtain S = .0807 / 0.464 = 0.17
%
% envelope, I anticipated that the entropy per spin at low temperature
% might be about $\log(2)/3 \simeq 0.23$, because roughly every third spin
% seems undetermined by its neighbours.}
The zero-temperature degeneracy is
nicely revealed by a plot of the free energy versus temperature which
has gradient at $T=0$ equal to minus the entropy. If the ground state
is unique the gradient is zero; for a triangular antiferromagnet the
gradient is non-zero. See figure \ref{fig1} for an illustration with
width $W=4$. This graph has gradient corresponding to a zero temperature
entropy of $0.17$ per spin. With $W=8$ I found an entropy of 0.088 per spin.
I have not figured out whether the ground state entropy per spin vanishes
as $W$ increases---maybe it
goes as $\sqrt{N}$ rather than as $N$.
The variance of the energy, which is proportional to $\partial E/
\partial \beta$ is the most interesting curve. In the rectangular
geometries and in the triangular case with $J=+1$, the fluctuations
show a {\em peak\/} as the temperature passes through
the critical temperature. This gives a hint of the phase
transition shown by the infinite system. The triangular system with $J=-1$,
on the other hand, has no such peak in its fluctuations. It looks
much the same as the Schottky anomaly.
\begin{figure}[h]
\[
\psfig{figure=ising/figs/ferr4.ps,width=3in,angle=-90}
\psfig{figure=ising/figs/anti4.ps,width=3in,angle=-90}
\]
\caption[a]{Free Energy per Spin of Long Thin Strip Ising Models.
Note the non-zero gradient at $T=0$ in the case of
the triangular antiferromagnet.
}
\label{fig1}
\end{figure}
The mean energy as a function of temperature is plotted in figure \ref{fig2}.
It is evaluated using the identity $\left< E \right> = - \partial \log Z /
\partial \beta$.
\begin{figure}
\[
\psfig{figure=ising/figs/ebar.4.ps,width=3in,angle=-90}
\]
\caption[a]{Mean Energy versus temperature of Long Thin Strip Ising Models with width 4.
}
\label{fig2}
\end{figure}
Figure \ref{fig3} shows the estimated heat capacity (taking raw
derivatives of the mean energy) as a function of temperature
for the triangular antiferromagnet with widths 4 and 8. It is apparent
that the peak in the heat capacity is getting sharper as the width
increases. The nature of any phase transition is not obvious though.
% better range would b 0.3 to 1
%
\begin{figure}
\[
\psfig{figure=ising/figs/anti.H.4.8.C.ps,width=3in,angle=-90}
\psfig{figure=ising/figs/ferr.R.4.8.C.ps,width=3in,angle=-90}
\]
\caption[a]{Heat Capacities of (a) Triangular Antiferromagnets with different
widths; (b) Rectangular Ferromagnets.
}
\label{fig3}
\end{figure}
% \medskip
% \hfill DJCM \today
\end{document}
From the partition function we can obtain interesting thermodynamic properties
using the relations (which you should confirm):
As the temperature goes to zero, the Gibbs distribution becomes
concentrated in the ground state. If the ground state is degenerate (\ie,
there are multiple ground states with identical energy)
then the entropy as $T \to 0$ is non-zero.