Two-Dimensional Ising Model


This program provides a Monte Carlo simulation of the two-dimensional Ising model using the Metropolis and Wolff algorithms. The output includes the mean energy, magnetization, susceptibility, and the heat capacity. The instantaneous energy and magnetization are plotted as a function of time. See the one-dimensional Ising model for a discussion of the Metropolis algorithm.

    One of the limitations of the Metropolis algorithm is the occurrence of critical slowing down near the critical temperature. The existence of critical slowing down is related to the fact that the size of the correlated regions of spins becomes very large near the critical point, and hence the time required for a region to lose its coherence becomes very long if a local dynamics is used.

    If we are interested only in the static properties of the Ising model, the algorithm we use to sample the states is irrelevant as long as the transition probability satisfies what is known as detailed balance. The Wolff algorithm flips a cluster of spins rather than a single spin, and is an example of a global algorithm. The utility of the Wolff algorithm is that it allows us to sample states efficiently near the critical temperature.

    In the following we summarize the main features of the Wolff algorithm. More details and its justification can be found in the references.

    A naive definition of a cluster of spins is that it corresponds to a domain of parallel nearest neighbor spins. We can make this correspondence more general by introducing a bond between any two parallel nearest neighbor spins. We assume that such a bond exists with probability p. For reasons that are discussed in the references, the naive definition of a cluster, which is equivalent to asuming that p = 1 must be replaced by a bond probability which depends on the temperature and is given by

p = 1 - e-2J/kT.

    The idea is to first generate a cluster of parallel spins with this bond probability and then flip the cluster with probability one. The algorithm can be summarized as follows.

  1. Choose a spin (the seed spin) at random. The nearest neighbor spins that are parallel to the seed site are the perimeter spins.
  2. Choose a perimeter spin and generate a random number r. If r ≤ p, a bond exists between the two spins, and the perimeter spin becomes part of the cluster. Otherwise, we do not test this bond again.
  3. If a spin is added to the cluster, add its perimeter sites to the list of perimeter sites of the cluster.
  4. Repeat steps (2) and (3) until no perimeter spins remain.
  5. Flip all the spins in the cluster.
Note that bonds rather than sites are tested so that a spin might have more than one chance to join a cluster.


  1. Use the Metropolis algorithm and choose the linear dimension L = 32 and the temperature T = 2. Describe the nature of the spin configurations. Are there large correlated regions where the spins are parallel? Increase T to T = 4 and describe your observations.
  2. Choose L = 4 and compute the mean magnetization, energy, susceptibility, and the specific heat as a function of temperature for T between 2.0 and 2.5 in steps of ΔT = 0.05. Also obtain data for T < 2.0 and T > 2.5, but use larger values of ΔT. For each temperature press the Zero averages button after about 1000 MCS and then average over at least 5000 MCS. To speed up the simulation set Steps per display to 100. Plot your results as a function of the temperature. What qualitative evidence is there for a phase transition?
  3. Repeat Problem 2 for L = 8, 16, and 32. Do the changes that you observe for larger L provide more evidence for a phase transition? If so, why?
  4. Finite size scaling. Compute the mean value of the absolute value of the magnetization per spin |m|, the specific heat C, and the susceptibility χ at T = Tc ≅ 2.269 for L = 4, 8, 16, and 32. Use as many Monte Carlo steps per spin as possible. Try both the Metropolis and Wolff algorithms and choose the algorithm which gives results that appear to converge more quickly. Plot the logarithm of |m| and χ versus L and use the finite size scaling relations |m| ∝ L-β/ν and χ ∝ L-γ/ν to determine the critical exponents β and γ. Use the exact result ν = 1. Do a log-log plot for the specific heat. You should obtain some curvature in your plot instead of a straight line because the specific heat critical exponent is zero. Instead show that C is proportional to ln L.
  5. Metastability and nucleation. Use the Metropolis algorithm and choose L = 64, T = 1, and magnetic field H = 0.7. Run the simulation until the system reaches equilibrium. You will notice that most of the spins are aligned (up) with the magnetic field. Pause the simulation and let H = -0.7; we say that we have "flipped" the field. Continue the simulation after the changed field and watch the configuration of spins. Do the spins align themselves with the magnetic field immediately after the flip? What is the equilibrium state of the system? You will probably notice that spins do not immediately flip to align themselves with the magnetic field. Instead most of the spins remain up and the mean values of the magnetization and energy do not change for many Monte Carlo steps per spin. We say that the system is in a metastable state. Why do the spins not flip as soon as the field is flipped? If you wait a while, you will see isolated "droplets" of spins pointing in the stable (down) direction. If a droplet is too small, it will likely shrink and vanish. In contrast, if the droplet is bigger than a certain critical size, it will grow and the system will quickly reach its equilibrium state. If the droplet is a certain critical size, then it will grow with probability 50%. This droplet is called the nucleating droplet. The initial decay of the metastable state is called nucleation.
  6. *Do a simulation at T = 10 until the number of up and down spins is approximately the same. Then change the temperature to T = 0.1 thus quenching your system. Step the simulation one Monte Carlo step per spin at a time and record the energy for each time step. A measure of the linear dimension of the domains that grow after the quench is given by R = 2/(2 + E/N). Plot log R versus the log of the time after the quench, and determine the functional form of your plot. If your plot is a straight line for long times, then the growth of R depends on time as a power law; the slope of the line is the power law exponent. If the growth of R follows a power law in the time, what is the exponent? Repeat for different system sizes and quench temperatures and determine if you always obtain approximately the same form for the domain growth.


Java Classes

Updated 28 December 2009.