Cluster algorithms and the 2D Ising Model

Introduction

The Ising Model is one of the simplest concrete examples of a statistical system which can be analyzed by theoretical and computational techniques. In one dimension it is rather straightforward to derive an exact solution (and we shall do so in class). In two dimensions the Ising model resisted solution for many years until it was solved by Onsager in the 1940's. The solution of the 3D Ising model is considered, by many, to be one of the most important problems of statistical physics. One of the most important techniques used to study this and similar models is numerical simulation. Such methods correspond to finding smart ways of doing integrals over a very large number of variables. State-of-the-art investigations of such systems deal with degrees of freedom and may consume hundreds of hours of CPU time on powerful supercomputers.

The Ising model is a model of a ferromagnet. In two dimensions and above it shows behavior typical of such a magnet; namely for temperatures T below the Curie temperature the system is in a ferromagnetic phase and has a spontaneous magnetization - the system behaves like a bar magnet. For temperatures the magnetization M vanishes and the system is in a paramagnetic phase. Close to the critical (Curie) temperature various thermodynamic quantities exhibit singularities - that is their values become infinitely large - the system exhibits a phase transition. In this project we will be concerned with examining some of these observables close to the phase transition for the 2D Ising model. To do so we shall utilize a recently discovered algorithm - the cluster algorithm which is far more efficient than the simple metropolis algorithm discussed in class.

The Model

The Ising model consists of a set of spin variables located at the sites of an square lattice in two dimensions. The spins can therefore be labeled as , where i and j are the indices or coordinates for the two spatial directions or as , where is a generic site label. Each of these spin variables can be in one of two states: up ( ) or down ( ). Thus each site can be thought of as being equipped with a little bar magnet which can have only two orientations. It is also reminiscent of the two states of a spin-1/2 electron - this is no real coincidence as it is indeed the phenomenon of electron spin that is responsible for the common ferromagnetic materials.

The energy of the system is conventionally written

Here, the notation means that the sum is over nearest-neighbor pairs of spins; these interact with a strength J. Thus, the spin at site (i,j) interacts with the spins at and . Notice that we assume periodic boundary conditions, so that, for example, the lower neighbors of the spins with i=L are those with i=1 and the lefthand neighbors of those with with j=1 are those with j=L; the lattice has the topology of a torus. This choice of boundary conditions can be shown to minimize the effects of using small lattices. The term in H represents the interaction of the spins with an external magnetic field.

When J is positive, the energy is lower if a spin is in the same direction (state) as its neighbors. We will be interested in the thermodynamics of the system. In this case it is convenient to measure the couplings J, H in units of the temperature, so that heating the system corresponds to decreasing these couplings. All the thermodynamic quantities can be obtained from a knowledge of the partition function Z

 

Monte Carlo Simulation

The partition function sum ranges over the states of the system. Even for a very modest system with say L=16 this number is huge ( ) so carrying out the sum in eqn. 2 is impossible. Even if we could perform this sum, the vast majority of states would contribute essentially nothing to Z - states with energies far from the mean are exponentially damped. Monte Carlo simulation is a technique for approximating the sum in eqn. 2 by considering only the most important contributions - configurations close to equilibrium. These are generated by a stochastic algorithm (one which involves random numbers). The probability distribution of these configurations can be chosen in such a way that the calculation of observables is trivial - we merely average them down the Monte Carlo sequence. We have discussed one such algorithm in class - metropolis. Here we consider a new and very powerful scheme - a cluster algorithm - which allows us to study the 2D Ising model with rather moderate computing resources.

Cluster Method

This scheme proceeds by identifying (in a stochastic, coupling dependent way) a cluster of parallel spins on some initial configuration. These are then flipped (spin direction reversed) and a new configuration obtained. By iterating this procedure a number of times a sequence of (nearly) statistically independent configurations are obtained typical of those that characterize equilibrium.

The creation of a cluster proceeds as follows (assume H=0):

1. For each link (pair of nearest neighbor sites) on the lattice examine the spins at each end. If they are parallel the link will be denoted active with probability . If not the link will be considered inactive.
2. Pick a site at random. The cluster is defined by all the spins that can be reached by moving out along active links from that site
3. Flip all those cluster spins.
4. Repeat from step 1.

Observables

You will need to measure the mean magnetization and the mean energy E

Also of interest are the fluctuations in these quantities - the specific heat C and susceptibility

It will be convenient to write the measured M and E to a file for later analysis. This should be done for each configuration.

Errors

Consider some quantity Q. The Monte Carlo procedure will generate a sequence where N is the total number of measurements. Our estimate for Q is then simply the average

If successive values were all independent the Q's error - could be calculated from the formula

 

However, in practice, each configuration remembers the configuration that generated it and so the 's are not statistically independent. Only after a sufficiently large number of intervening updates can we safely assume that a configuration is not correlated with a previous one. This results in the naive formula eqn. 7 yielding typically too small an error estimate for Q.

This is tackled by so-called data binning. Divide the sequence into sets each containing P consecutive values and calculate the average in each set. It is assumed that P divides N and so this procedure yields new, partially-averaged variables .

For large enough P these new variables will be uncorrelated and formula
 7 can be used now with , . In practice we choose and compute a set of estimates for the error by letting . For large l should approach a constant - the true error.

Simulation

The exact solutions for the energy E and magnetization M for the infinite lattice are as follows.

The variable and the couplings and are given by

The function is called a complete elliptic integral of the first kind and is given by

1. Write a code to implement this cluster algorithm. Use as a basis the code you downloaded for the 307 lab. You just need to change the update function so that it no longer uses a metropolis algorithm to change the spins but a cluster algorithm.
2. (Optional)Plot your simulation data for L=32 against the exact solution (use Simpson's rule to evaluate the elliptic integral) for H=0 and in steps of 0.05. You will need at least 1000 cluster updates per data point.
3. Choose L=8, L=16 and L=32. Plot out the specific heat and magnetic susceptibility data over this range of coupling. What do you notice for ?

The nature of the singularities is encoded in so-called critical exponents. For example

These exponents can be measured using the theory of finite-size scaling. In essence this states that if simulations are performed at , these quantities scale with L as

The quantity is another exponent. It can be related to for many systems by the so-called hyperscaling relation (d is the dimension).

• By running at the transition coupling for different L estimate the critical exponents and . To extract an estimate of, for example, you will need to plot the log of against the log of L. This should be a straight line whose gradient will be the exponent ratio. It is not necessary to put errors into your analysis.