Simulation of the 2D Ising model

One of the most interesting phenomena in nature is ferromagnetism. A FM material exhibits a non-zero spontaneous magnetization in the absence of an applied magnetic field. This occurs below a well-defined critical temperature $T_c$ known as the Curie temperature. For $T > T_c$ the magnetization vanishes. Hence $T_c$ separates two phases, a disordered one for $T>T_c$, and a ferromagnetic one for $T<T_c$.

Although the Ising model is too simple, it already contains much of the physics of the FM phase transition. In order to explore the properties of this model, we need to calculate some physical quantities of interest, including the mean energy $\langle E \rangle$, the mean magnetization $\langle M \rangle$, the heat capacity $C$, and the magnetic susceptibility $\chi$.

The Ising model

Consider a lattice with $N$ sites, where each site $i$ can assume two possible states $s_i=+1,-1$, or spin “up” and spin “down”. A particular configuration or microstate of the lattice is specified by the set of variables $\{s_1,s_2,...s_N\}$ for all lattice sites.

Now we need to know the dependence of the energy $E$ of a given microstate, according to the configuration of spins. The total energy in the presence of a uniform magnetic field is given by the “Ising model”: $$E=-J\sum_{\langle ij \rangle}s_is_j-h\sum_{i=1}^Ns_i,$$ where the first summation is over all nearest neighbor pairs and the second summation is over all the spins of the lattice. The “exchange constant” $J$ is a measure of the strength of the interaction between nearest neighbor spins. If $J>0$, the states with the spins aligned $\uparrow \uparrow$ and $\downarrow \downarrow$ are energetically favored, while for $J<0$ the configurations with the spins antiparallel $\uparrow \downarrow$ and $\downarrow \uparrow$ are the ones that are preferred. In the first case, we expect that the state with lower energy is “ferromagnetic”, while in the second case, we expect it to be “antiferromagnetic”. If we subject the system to a uniform magnetic field $h$ directed upward, the spins $\uparrow$ and $\downarrow$ possess and additional energy $-h$ and $+h$ respectively. Note that we chose the units of $h$ such that the magnetic moment per spin is unity.

Instead of obeying Newton’s laws, the dynamics of the Ising model corresponds to “spin flip” processes: a spin is chosen randomly, and the trial change corresponds to a flip of the spin $\uparrow \rightarrow \downarrow$ or $\downarrow \rightarrow \uparrow$.

Physical quantities

The net magnetic moment or "magnetization" $M$ is given by $$M=\sum_{i=1}^N s_i.$$ Usually we are interested in the average $\langle M \rangle$ and the fluctuations $\langle M^2 \rangle - \langle M \rangle ^2$ as a function of the temperature of the system and the applied magnetic field. Notice that in the absence of a magnetic field, the value of $M$ should average zero. This is because the system is as likely to be found with all the spins pointing up or down. Therefore, one typically measures $m=\sqrt{M^2}$ or $m=|M|$, which is always a positive quantity.

The heat capacity

One way to measure the heat capacity at constant external field id from the definition: $$C=\frac{\partial \langle E \rangle}{\partial T}.$$ Another way is to use the statistical fluctuations for the total energy in the canonical ensemble: $$C=\frac{1}{(k_BT)^2}\left( \langle E^2 \rangle - \langle E \rangle ^2 \right).$$

The magnetic susceptibility

The magnetic susceptibility $\chi$ is an example of a “response function ”, since it measures the ability of a spin to “respond” or flip with a change in the external magnetic field. The zero isothermal magnetic susceptibility is defined by the thermodynamic derivative $$\chi=\lim _{H \rightarrow 0} \frac{\partial \langle M \rangle}{\partial H}.$$ The zero field susceptibility can be related to the magnetization fluctuations in the system: $$\chi=\frac{1}{k_BT} \left( \langle M^2 \rangle - \langle M \rangle ^2 \right),$$ where $\langle M^2 \rangle$ and $\langle M \rangle ^2$ are zero field values.

Metropolis algorithm

Moves

One typically picks a random spin and it is flipped or not by calculating the energy difference between the considered spin and its 4 nearest neighbors using the formula:

The energy change is $\Delta U = 2U$. The spin then directly flips if $\Delta U \le 0$. Otherwise, it only flips if a randomly chosen number between 0 or 1 is smaller than the Boltzmann factor $\exp(-kT\Delta U)$.

Notice that the sum in brackets can run from -4 to +4, so we can easily store the values for all the possible configurations in lookup tables.

Boundary conditions

Since we are interested in the properties of an infinite system, we have to consider the boundary conditions. The simplest case corresponds to “free boundary condition”. In a 1D chain, the spins at sites $1$ and $N$ are open ends and have one interaction each. In general a better choice is periodic boundary conditions, where sites 1 and $N$ interact with each other closing a loop. In this situation, the chain has the topology of a ring, and all the spins have the same number of interactions. We also say that there is translational invariance, since the origin can be picked arbitrarily.

As discussed earlier, the use of PBC minimizes the finite size effects. However, a disadvantage is that they reduce the minimum separation between two spins to half the length of the system. In a 2D system we need to account for the periodicity on both the $x$ and $y$ directions, and our system will have the topology of a torus.

Initial conditions and equilibration

We can pick a random initial configuration. However, as we shall see, in some simulations the equilibration process can account for a substantial fraction of the total computer time. The most practical choice of an initial condition is an “equilibrium” configuration of a previous run which is at a temperature close to the desired temperature.

Tricks

It is convenient that we store all the transition probabilities in lookup tables, so we do not have to calculate them at each step. Another trick consists of storing all the positions of the spins and their neighbors to avoid calculating many random numbers. If you need to perform several runs for different values of the temperature $T$, you can do it at the same time using the same random numbers.

Exercise 11.2: Equilibration of the 2D Ising model

1. Run your simulation with $L=8$ and $T=2$ and choose the initial spins to be all up. Plot the variation of the energy and the magnetization with time. How much time is necessary for the system to reach equilibrium?

2. Visually inspect several “equilibrium” configurations. Is the system ordered or disordered?

3. Run the program with $T=1.5$ and chose the same initial configuration with all the spins up. How long does it take for the system to reach equilibrium?

4. Visually inspect several equilibrium configurations with $T=1.5$. Are they more or less ordered than those in part 2?

5. Is the acceptance ratio and increasing or decreasing function of $T$? Does the Metropolis algorithm become more or less efficient at low temperatures?

Challenge 11.2

Modify the code to calculate the correlation function

Plot the correlations along the $x$ direction across the phase transition and describe its behavior.