8. 2D Ising model

For a system of interacting spins on a 2D lattice of sites, the Hamiltonian describing the energy of interaction is $$ {\cal H} = - \sum_{i,j} J_{ij} s_i s_j -\mu \sum_j h_j s_j $$ where the index $j$ runs over all spins on the lattice, while the index $i$ runs over those spins with which spin $s_j$ interacts, with the interaction strength $J_{ij}$. In the case of spin-1/2, each $s_j$ can take on values of +1 or -1. Typical example is a nearest-neighbour interaction, on a square lattice implying interactions with spins directly North, South, East and West of site $j$. The sign of the interaction $J_{ij}$ determines whether the interaction favors (lowers the energy for) spins aligned parallel to their neighbours ($J_{ij}>0$, ferromagnetic), or anti-parallel to their neighbours ($J_{ij}<0$, antiferromagnetic). In the absence of spin-spin interactions, $J_{ij}=0$.

The second sum describes the interaction of spins with an external field: $\mu$ is the unit of strength (a magnetic moment of a single spin) and $h_j$ represents the strength pf the field at site $j$. Positive $h_j$ value implies that the alignment along the external field is energetically favourable, $h_j < 0$ favours spins aligning against the field, and $h_j=0$ indicates the absence of the external field.

The probability of encountering a particular configuration $\{s\}$ of spins on the lattice is given by the Boltzmann factor: $$ P(\{s\}) = \frac{e^{-{\cal H(\{s\})}/k_B T}}{\sum_{\{s\}} e^{-{\cal H(\{s\})}/k_B T}} $$ with the sum over all possible configurations of the system in the denominator being the partition function $Z$. Often the notation $\beta = 1/k_b T$ (the inverse temperature) is used; with that, $$ P(\{s\}) = \frac{1}{Z} e^{-\beta{\cal H(\{s\})}} $$ In two or more dimensions, the Ising model exhibits interesting phase behaviour, including ferromagnetic phase transitions, and is widely studied.

The limitations of finite size of the model lattice of spins are normally overcome by imposing periodic boundary conditions.

Set up a 2D rectangular lattice of spins

Initial state will be a random orientation of $2N+1 \times 2N+1$ spins. Here we assume that all sites can only have the magnetic moment of $+\mu$ or $-\mu$. rand() produces a random number somewhere in the 0..1 range, and and a shift of -0.5 makes it equally likely for the sign to be up or down.

To visualize the state of the spin lattice, we can use colour to enhance the previous picture. While we are at it, we can make the function a little more robust by adding error checking:

System evolves in time according to the rules encoded in Evolve()

Just as a demonstration, go through all spins in the system and randomly "flip" some of them. Setting the cutoff at the midpoint of the 0..1 interval should flip about a half of them. Going over the entire spin system once corresponds to a single step in the time evolution of the system.

Monte-Carlo and $\pi$

The principle of Monte-Carlo sampling can be well illustrated by the following method of determining $\pi$, from the ratio of the area of a unit circle, $S_c = \pi r^2 = \pi$ and the area of the square that circumscribes it, $S_s = 2\times 2 = 4$. Imagine randomly dropping stones anywhere inside the square, but counting how many stones fall inside the circle ($N_{in}$), and how many outside of the circle ($N_{out}$). If the probability of dropping a stone anywhere inside the square is equal, then the fraction of the total number of stones that fall inside the circle will correspond to the fraction that the area of the square that is occupied by the circle, i.e. $$ \frac{N_{in}}{N_{in}+N_{out}} = \frac{S_c}{S_s} = \frac{\pi}{4} $$

Metropolis algorithm

Given a configuration $\{s\}$, the evolution of the system is always toward lower total energy. Therefore, for each spin one needs to:

Plan of action