# Ising model and Metropolis Monte Carlo

For many problems of interest, it is much too difficult to compute an exact solution. For example, the exact partition function of the Ising model on a 3D lattice of 10x10x10 particles contains $2^{1000}\approx10^{300}$ terms in the sum. This is nearly the **fourth power** of the number of atoms in the visible universe. This partition function doesn't seem to have a simple exact formula, either, despite decades of effort by many researchers. The two options most easily accessible to us--brute force computation and elegant exact solution--are therefore unavailable to us. How are we to extract the physics of such a model, let alone a more complicated model that captures other neglected interactions?

The idea of Monte Carlo sampling is that we don't *need* to compute all of the terms in the Ising model partition function, or any other function for that matter. We just need to understand the "typical" behavior of that function.

As an example, consider the following. There are about 16 million registered voters in Australia. If we want to try to guess the outcome of the next election, do we need to query every single voter? No! It suffices to ask a **random sample** of voters what their preferences are. If that sample is not too small, and not too biased, then it should give accurate results.

If the sample *was* too small, then the statistical noise from taking so few voters would be too great to infer the outcome of a close election. Mathematically speaking, the variance of our estimate would be too large to extract useful results. Similarly, if our sample was too biased--for example, you only sampled voters in Dapto, or in Darlinghurst--you would not get an accurate picture of the federal election.

This idea is very simple and intuitive, but it is surprisingly powerful.

## Ising Model

The Ising model is the simplest model of a ferromagnet. The basic phenomenology of the Ising model is simple: there is a certain temperature $T_c$ below which the system will spontaneously magnetize. This is what we'll study with Monte Carlo.

We consider a system of interacting spins on a square lattice $\Lambda$. (The one-dimensional lattice case is too easy, and the three-dimensional case is surprisingly hard, so we'll stick to two dimensions.) The model is defined by choosing a Hamiltonian for the system. For every spin configuration, we'll assign a corresponding energy using the following rule. If two neighboring states are aligned, then we add $-J$ units of energy, and if they are anti-aligned add $+J$ units of energy. We can summarize this state of affairs by defining spin variables $s_i$ to every vertex in the lattice, where $s_i = \pm1$ are the possible values. Then the Hamiltonian is given by \[H = -J\sum_{(i,j)\in \Lambda} s_i s_j \,.\] Here the notation $(i,j)\in \Lambda$ simply means a sum over all neighbors (edges) on the lattice $\Lambda$.

The average magnetization density for the Ising model is given by \[\avg{M} = \frac{1}{N^2} \sum_{i\in \Lambda} \avg{s_i}\,,\] where the sum is over all spins in the lattice (a total of $N\times N = N^2$ spins), and the average in brackets is a thermal average with respect to the Gibbs state at temperature $T$. This quantity is implicitly a function of temperature and the coupling strength $J$ (as well as the number of spins and their geometry).

In the context of the Ising model, we can use Monte Carlo methods to compute various quantities of interest, such as the magnetization density or the energy density as a function of temperature. Because the 2D Ising exhibits much of the phenomonology of ferromagnetism but still admits an exact solution, it is a good test bed for Monte Carlo because we can compare the simulations to the exact answers.

## Metropolis sampling rule

Suppose we want to estimate the average magnetization density in thermal equilibrium at some given temperature. The problem is, we don't know the probabilities $P(s)$ for each of the states $s$. But we can *generate* a sample from this probability distribution using Monte Carlo methods.

Consider the following simple algorithm, called the "heat bath algorithm" for converging an initial spin configuration to the thermal equilibrium distribution.

#### Heat Bath Monte Carlo Sampling

- Pick a spin $j$ at random
- Compute $E(\pm1) = \pm \bigl(4 - 2\cdot [\# \text{ of +1 neighbors}]\bigr)$, the energy associated to spin $j$ in the $\pm1$ state given its current environment
- With probability $\P(\pm1) = \e^{-E(\pm1)/kT}/Z$, set spin $j$ to be $\pm 1$, where $Z = \P(+1) + \P(-1)$.
- Repeat

This simple algorithm uses two ideas: the mean-field approximation, and Monte Carlo sampling. It treats each spin one at a time in the mean-field approximation and randomly assigns it a spin direction consistent with the Boltzmann distribution.

This algorithm is an example of a **Markov chain**. A Markov chain is a probabilistic process where the next state depends probabilistically only on the current state of the system. There is a lot of elegant mathematics that has been developed to help us analyze Markov chains, but that is beyond the scope of this course. However, there are some important general features which hold for our Monte Carlo methods and which are crucial for their success. (Technically, we are only considering finite time-homogeneous Markov chains.)

A Markov chain is said to be **ergodic** if the probabilities of transitioning between states are such that any state can eventually transition to any other state (given enough time). A powerful theorem known as the Perron-Frobenius theorem asserts that for ergodic Markov chains, there is a *unique* probability distribution which is invariant under the transition rules. That is, there is a unique *equilibrium distribution*. Finally, the transition probabilities between states in the equilibrium distribution are said to obey the **detailed balance** condition if the transition probability from $s \to s'$, which we denote by $P_{s \to s'}$, acts on the equilibrium probability vector $\rho_s$ (a vector of probabilities for states $s$) and satisfies the equation
\[P_{s\to s'}\rho_s = P_{s'\to s}\rho_{s'} \,.\]
That is just a complicated way to say the following: in equilibrium, the conditional transition rates between two states are always *equal*. Note that this is *only* generally true in equilibrium, which is a good thing if you think about it, because we generally wish to transition preferentially from high energy states to low energy states. A Markov chain that obeys the detailed balance condition is sometimes (equivalently) said to be a **reversible Markov chain**.

The heat bath algorithm listed above satisfies the general rules outlined above: it is a Markov chain because the probabilities to transition to the next state depend only on the current state; it is ergodic because every spin has a chance to flip hence every state is reachable; the Boltzmann distribution is, by construction, an equilibrium distribution; it is therefore a unique equilibrium distribution. Thus, if we run our Monte Carlo algorithm long enough for all the transient effects to die down, then we will eventually reach the Boltzmann distribution that we wish to sample from. An important question here is *how long we need to wait* for those transients to die, and this idea can be quantified into the **mixing time** of the Markov chain (more on this later).

There is one more crucial idea for Monte Carlo sampling that will take us from the simple heat bath algorithm to the celebrated Metropolis algorithm. When we decide to flip a spin, we have a probability to take it from high to low energy, or vice versa. The idea of Metropolis sampling is: why should we ever reject a move that takes us to lower energy? After all, this is where we'd like to be. It makes sense that we sometimes accept a move to higher energy; that way we can try to get out of certain "traps" where the total energy is only *locally* a minimum. This idea is an example of **importance sampling**, a powerful general principle for Monte Carlo methods.

Let's state the specific Metropolis algorithm.

#### Metropolis Monte Carlo Sampling

- Pick a spin $j$ at random
- Compute $\Delta E$, the change in energy associated to spin $j$ if we flipped just that spin given its current environment
- If $\Delta E \lt 0$, flip spin $j$; if $\Delta E \gt 0$, flip spin $j$ with probability $\e^{-\Delta E / kT}$.
- Repeat

The only difference from the heat bath is that we *always* accept moves that lower the energy. It turns out that the Metropolis algorithm obeys all of the same nice properties listed above for the heat bath algorithm (it is a reversible ergodic Markov chain, hence the Gibbs distribution is the unique fixed point), but it is more efficient because it spends less time at higher energies.

## Ising model simulations in Matlab

Download the package `ising.zip` from the course website. Unzip the file and add the folder to your Matlab path. This folder contains the following files:

`ising.m`: the main file which contains a simulator for the 2D Ising model using Metropolis sampling and Wolff sampling, including parameters for temperature and coupling strength.`metropolis.m`: implements Metropolis sampling on a configuration of Ising spins for some given amount of time.`wolff.m`: implements Wolff sampling on a configuration of Ising spins for some given amount of time.`isingenergy.m`: compute the energy density of a spin configuration.`isingplot.m`: plot a configuration of spins.`correlations.m`: compute the 2D connected correlation function for a spin configuration.`radialavg.m`: average out the angular dependence from a 2D connected correlation function.

Now just play around with the code! Try to keep $N$ (the linear number of spins in the lattice) at a reasonable size; bigger is better up to a point, but you don't want to slow down too much. $N$ between 50 and 65 seems to work well on my laptop.

Here are some exercises to get you thinking about the physics of the Ising model.

- Can you locate the phase transition as a function of $kT/J$ to within, say, .25?
- How long does it take to reach a (roughly) stable configuration of spins as you approach the phase transition from below?
- Do you see any
*metastable*configurations? For example, states with low energy, but also low magnetization because they break up into magnetic*domains*? Why might it take a long time for such configurations to eventually reach a state of complete magnetization? - What happens if you change the sign of $J$? What do the ground states look like? What do the domains look like?
- Suppose you initialize the spins to be 45% up and 55% down. How often do you end up in the spin-down magnetized state when you start the simulation below the critical temperature? What thermodynamic quantity would control the probability of eventually landing in the up or down magnetized state?
- Plot the correlation function for a given spin configuration, defined as \[\frac{1}{N^2}\sum_{i,j}s_{i,j} s_{i+x,j+y} - \avg{M}^2\,,\] where $\avg{M}$ is the magnetization density. What do the correlations look like in a (mostly) magnetized state? What about in a metastable state? In the high-temperature phase?