Simulation of an Ising 2D system using the Monte Carlo method and the Metropolis algorithm

This post provides a practical guide for creating a Monte Carlo program to simulate the square lattice Ising model. The example program is written in JavaScript and is intentionally kept simple. However, for these types of simulations, you might prefer to use programming languages such as C++, Rust, C, or Fortran.

You can find the code used for this simulation at this link.

Let’s start with some important concepts you may find useful before reading this post.

The Ising model is a simple and important idea in statistical physics, first created to explain how ferromagnets work. In a ferromagnet, each atom acts like a small “spin” that can point either up or down. These spins want to be next to other spins pointing the same way because of a quantum effect called the exchange interaction. When the temperature is very low, all the spins align in the same direction to reach the lowest possible energy. However, as the temperature increases, the atoms “move” more because of the heat energy. This thermal “motion” breaks their alignment, so the spins start to flip and become less ordered.

We need a way to measure how “ordered” the spins are. For this, we use a property called magnetization. At low temperatures, magnetization is high. At high temperatures, magnetization is low. There is a temperature where magnetization changes drastically. This is the critical temperature. The critical temperature (Curie temperature for ferromagnets) is the temperature at which a material loses ferromagnetic order and becomes paramagnetic.

Having said that, let’s describe the Ising model. The Ising model is a simple lattice model of magnetism where each site $i$ carries a spin $s_i = \pm 1$ (up or down), representing a tiny magnetic moment. Neighboring spins interact with energy $E_i = -J s_i s_j$, where $J$ is the exchange constant. If $J>0$, aligned neighbors (both +1 or both −1) are energetically favored, mimicking ferromagnetism.

It is worth mentioning that the Ising model can be applied in 1D, 2D, and 3D. For the sake of simplicity, we will focus on the 2D model, but the same principles apply to the 1D and 3D cases.

In 2D, each site at position $\left(i, j\right)$ on the square lattice has exactly four neighboring sites: $\left(i + 1, j\right)$ (below), $\left(i - 1, j\right)$ (above), $\left(i, j + 1\right)$ (to the right), and $\left(i, j - 1\right)$ (to the left). If one of these neighboring indices would fall outside the lattice, we use periodic boundary conditions: the lattice “wraps around” so that an edge site is neighbors with the site on the opposite edge. For example, in a lattice of size $L$ (where indices go from $0$ to $L-1$), the site at $i = L-1$ has $i = 0$ as its neighbor below, and similarly, the top and bottom rows are neighbors due to this wrapping. This creates a lattice surface shaped like a torus, where every site always has four neighbors, regardless of its position.

The local energy at a site (assuming all neighbor interactions have the same $J$) is calculated as follows:

$E_{(i,j)} = -J \left( s_{(i,j)} s_{(i+1, j)} + s_{(i,j)} s_{(i-1, j)} + s_{(i,j)} s_{(i, j+1)} + s_{(i,j)} s_{(i, j-1)} \right)$

$E_{(i,j)} = -J s_{(i,j)} \left( s_{(i+1, j)} + s_{(i-1, j)} + s_{(i, j+1)} + s_{(i, j-1)} \right)$

For simplicity, let’s refer to the site $\left(i, j\right)$ as site $i$. Please do not be confused by this. It is still a 2D lattice.

$E_i = -J s_i \displaystyle \sum_{\alpha} s_{\alpha}$

where $\sum_{\alpha}$ means the sum over the nearest neighbors.

The total energy $E$ is given by:

$E = \displaystyle -\dfrac{J}{2} \sum_i \sum_{\alpha} s_i s_{\alpha}$

Note that we divide by 2 because each interaction is counted twice: if site $\alpha$ is a neighbor of site $i$, then site $i$ is also a neighbor of site $\alpha$.

The magnetization $M$ can be calculated as:

$M = \displaystyle \sum_i s_i$

You’ve read quite a bit so far (I hope!). Let’s look at some diagrams to make these concepts clearer. In the plot below, you can see a

From now on, indices increase to the right and downward, and decrease to the left and upward. Specifically, index $i$ increases downward (and decreases upward), and index $j$ increases to the right (and decreases to the left).

Let’s focus on the site at position $(0,0)$ (top-left corner). The site at position $(0,0)$ has a spin of value . Its neighbors are located at the following positions (using periodic boundary conditions):

• $(1, 0)$ (below)
• $(0, 1)$ (right)
• $(L-1, 0)$ (above, wraps around) $=(4, 0)$
• $(0, L-1)$ (left, wraps around) $=(0, 4)$

as shown here:

The spins at these neighbor sites are:

• $(1, 0)$:
• $(0, 1)$:
• $(4, 0)$:
• $(0, 4)$:

The local energy at site $(0,0)$ is $J$. You can calculate it easily by adding the number of neighboring spins that point in the same direction, subtracting those that point in the opposite direction, and multiplying the result by $-J$. The lattice has a total energy of $J$ and a magnetization of . You can compute the magnetization by counting the number of spins pointing up and subtracting the number of spins pointing down.

Metropolis algorithm

In a small lattice of $10 \times 10$ sites, there are a total of $2^{100}$ possible configurations, which is approximately $10^{30}$. Even if you had a very fast computer that could process $10^9$ configurations per second, it would still take $10^{21}$ seconds (about 1000 times the age of the universe) to analyze all possible states. Don’t worry about it! We have learned to solve this problem by performing importance sampling using the Metropolis algorithm.

We use the Metropolis algorithm to simulate the thermal fluctuations of the 2D Ising model at a given temperature $T$. The algorithm proceeds as follows:

1. Choose a random site $i$ in the lattice.

2. Flip the spin at site $i$ ($s_i \rightarrow -s_i$), and calculate the change in energy $\Delta E$ that this flip would cause. You may compute only the local energy $E_i$ at site $i$.

3. Decide whether to accept the flip:

   • If $\Delta E \leq 0$ (the energy decreases), always accept the flip.
   • If $\Delta E > 0$, accept the flip with probability $e^{-\Delta E / (k_B T)}$, where $k_B$ is the Boltzmann constant (here set to 1 for simplicity).

4. Repeat steps 1–3 $N$ times (one “Monte Carlo step”), where $N$ is the number of sites in the lattice.

5. Repeat step 4 (perform many Monte Carlo steps) until the system reaches equilibrium.

Below you can interact with a simulation at constant temperature. You can advance the simulation step by step, or press play to watch how the lattice evolves over time.

Temperature annealing

Now let’s see what happens when we vary the temperature. In the process called “annealing,” we slowly decrease the temperature and observe how the system evolves. This helps us understand how the magnetization and energy change as the system moves from a disordered (high-temperature) state to an ordered (low-temperature) state.

Two important physical quantities measured in Monte Carlo simulations of the Ising model are the magnetic susceptibility and the specific heat.

The magnetic susceptibility, usually denoted $\chi$, quantifies how strongly the magnetization of the system responds to an applied external magnetic field (even though we simulate at $H=0$ in our model, $\chi$ can be computed from equilibrium fluctuations via statistical mechanics). In practical simulations, it is computed via fluctuations of the magnetization:

$\chi = \frac{1}{k_{B} T} \left( \langle M^2 \rangle - {\langle M \rangle}^2 \right)$

where $T$ is the temperature, and angle brackets denote averaging over the simulation.

The specific heat $C_v$ describes how much the energy of the system fluctuates as we change the temperature:

$C_v = \frac{1}{k_{B} T^2} \left( \langle E^2 \rangle - \langle E \rangle^2 \right)$

where $E$ is the total energy of the system. Again, this is obtained from the variance of the energy measured during simulation.

Both quantities are examples of response functions: susceptibility measures the response of the magnetization to a field, while specific heat measures the response of energy to temperature. As you change the temperature, you will observe that both $\chi$ and $C_v$ typically exhibit sharp peaks near the critical temperature $T_c$, indicating a phase transition from a disordered to an ordered phase.

Finally, we have a simulation that explores a range of temperatures. You can use the controls to observe how the magnetization approaches $N$ (the total number of sites) at low temperature, and trends toward $0$ at high temperature. Note that we work in units where $k_B = 1$, so the dimensionless coupling is $J/(k_B T)$. You can also choose $J = -1$ to simulate antiferromagnetic behavior, but we will discuss this case later.

Note that the peaks happened around 2.26, which was demonstrated to be the theoretical critical temperature:

$\dfrac{k_B T_c}{J} = \dfrac{2}{\ln{\left( 1 + \sqrt{2} \right)}} \approx 2.26918531421$

You can find the code used for this simulation at this link.