\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.