Ising

Introduction

Ising model was developed by Wilhelm Lenz and Ernst Ising as a simplified model to describe ferromagnetism in the perspective of statistical mechanics. It is one of the simplest models that experience phase transition. In short, it models magnetism as a process of electrons aligning their magnetic field with their immediate neighbors. This has similarity with flocking model and consensus problem in multi-robot system in that global properties are derived only from local interactions between the agents.

In this project, we will first introduce the mathematical foundation of the Ising model, then analyze the at equilibrium behavior of the 2D lattice Ising model, finally, we will show the simulation of 2D lattice Ising model using the Metropolis-Hasting algorithm.

Mathematical Formulation

Electron State

We are interested in the state of each electron called spin, denoted $x$. The spin can take on two values: 1) up ($1$), and 2) down ($-1$)

$$ x = \begin{cases} 1 & \text{(up)}\\ -1 & \text{(down)}\end{cases}$$

Electron Sea

Consider a sea of $N$ immobile electrons (such as in a block of iron), each electron $n$ has local interaction $J_{n,m}$ with every neighboring electron $m$. Additionally, each electron $n$ experience an external electric field $h_n$. This electron sea can be represented as a graph as follow:

Alignment Energy

The local interaction encodes the energy of misalignment between the pair of electrons, if the electrons are aligned (same spin), the energy is small; if the electrons are misaligned, the energy is high. In the Ising model, this is done by multiplying the spin of the pair

$$ E(x_n, x_m) = -J_{n,m}x_nx_m $$

If the electrons are aligned, the energy is $-J_{n,m}$; if the electrons are misaligned, the energy is $J_{n,m}$. Likewise, the energy associated with the external field is the alignment of the electron with its local external field:

$$E(x_n) = h_{n}x_n$$

The state of the entire system can be defined as a vector of the spin of every electron in the electron sea,

$$ X_i \triangleq (x_1, x_2, \cdots , x_n, \cdots, x_N) $$

The energy of any given system state is the sum of all alignments energies

$$ E(X_i) = -\sum_{(m,n)}J_{n,m}x_nx_m + \sum_{n} h_nx_n $$

State Probability

From statistical mechanics, we know that the probability of a system existing in a particular state $X_i$ follows the Boltzmann distribution

$$\mathbb{P}(X_i) = \frac{1}{Z}\exp(-\beta E(X_i))$$

where $\beta$ is the inverse temperature of the system, and $Z$ is the normalizing constant, aka the partition function,

$$Z\triangleq \sum_i \exp(-\beta E(X_i))$$

Magnetization

Finally, the macro quantity we are interested about the system is the magnetization $M$, which is the average spin of the entire electron sea

$$ M \triangleq \frac{1}{N}\sum_n x_n $$

If all the electrons has a up-spin $x=1$, the magnetization is $M=1$. If all the electrons has a down-spin, the magnetization is $M=-1$. If half of the electrons have a down-spin and other half have a up-spin, the magnetization is $M=0$.

Metropolis-Hasting Simulation

With the probability of the system being in each state defined, we can simulate the system using the Metropolis-Hasting (M-H) algorithm, which is in the class of algorithm called Markov chain monte carlo. The idea is to randomly perturb the state of the system and decide whether to keep the perturbation based on the probability of the state. By properly selecting this acceptance rate of perturbation, the M-H algorithm can ensure the frequency of the system being in each state during the simulation to be equal to the true probability defined by the Boltzmann distribution.

Detailed Balance

The M-H algorithm relies on a condition called detailed balance. In short, it is a condition where the in-flux and out-flux between two states are the same for every pair of states

$$ \mathbb{P}(X_i)\mathbb{P}(X_j | X_i) = \mathbb{P}(X_j)\mathbb{P}(X_i | X_j) $$

The goal is to design the transition probability $\mathbb{P}(X_j | X_i)$ such that the detailed balance condition hold for every pair of states for $\mathbb{P}(X_i) = \frac{1}{Z}e^{-\beta E(X_i)}$.

Acceptance Rate

The M-H algorithm is separates the transition probability into two steps: 1) propose $g(X_j|X_i)$ 2) accept $A(X_j|X_i)$. The total probability of transition is the product of the two.

$$ \mathbb{P}(X_j | X_i) = g(X_j|X_i)A(X_j|X_i) $$

The detailed balance can be rewritten as

$$ \begin{aligned} \mathbb{P}(X_i)g(X_j|X_i)A(X_j|X_i) &= \mathbb{P}(X_j)g(X_i|X_j)A(X_i|X_j) \\ \frac{A(X_j|X_i)}{A(X_i|X_j)} &= \frac{\mathbb{P}(X_j)g(X_i|X_j)}{\mathbb{P}(X_i)g(X_j|X_i)} \end{aligned} $$

Single Spin Flip Dynamic

The particular propose transition we use here is the single-spin-flip model, where we uniformly randomly select one of the electron, and flip its spin, then the proposed transition probability is just the chance of selecting a particular electron

$$ g(X_j|X_i) = g(X_i|X_j) = \frac{1}{N} $$

This simplifies the acceptance ratio to

$$ \begin{aligned} \frac{A(X_j|X_i)}{A(X_i|X_j)} & = \frac{\mathbb{P}(X_j)}{\mathbb{P}(X_i)} \\ & = \frac{\frac{1}{Z}e^{-\beta E(X_j)}}{\frac{1}{Z}e^{-\beta E(X_i)}}\\ & = e^{-\beta (E(X_j) - E(X_i))} \\ & =e^{-\beta \Delta E} \end{aligned} $$

Since this is a single flip, the change of energy can be evaluated only around the neighborhood of the selected electron.

There are still two variables, $A(X_j|X_i), A(X_i|X_j)$ with one equation. The standard method is to set $A(X_i|X_j)$ to $1$ if $\Delta E$ is positive; and $A(X_j|X_i)$ to $1$ if $\Delta E$ is negative. In other word, always transition if it leads to a decrease in energy, and sample from $e^{-\beta \Delta E}$ if the transition leads to an increase of energy. Mathematically,

$$ \begin{aligned} A(X_j|X_i) &= \begin{cases} 1 & \text{if } \Delta{E} \leq 0 \\ e^{-\beta \Delta E} & \text{if } \Delta{E} > 0 \end{cases}\\ &= \min\left\{1, e^{-\beta \Delta E}\right\} \end{aligned}$$

With this we have the full algorithm for simulating the Ising model

Full Algorithm

Uniformly randomly select electron

$$ n \sim \text{uniform}(\{1,2,\cdots, N\}) $$
Calculate energy different after flip

$$ \Delta E = -2Jx_n\sum_{m\in \mathcal{N}(n)} x_m $$
Calculate acceptance rate

$$ A = \min\left\{1, e^{-\beta \Delta E}\right\} $$
Sample acceptance rate, flip the spin if $a$ is smaller then $A$, reject the flip otherwise

$$ a \sim \text{uniform}([0, 1]) $$

$$ X = \begin{cases} X_j & \text{if } a\leq A \\X_i& \text{if } a> A \end{cases}$$

Simulation

We will simulate the system by using $J=1$ on a $100\times 100$ grid. The neighborhood of each electron is a $5\times 5$ square centered at the electron. Additionally, we slowly increase $\beta$ to observe the magnetization of the system.

We can see a sudden change of magnetization from $0$ to $1$ at around $\beta = 0.04$, this is called a phase change, where a system shows a sudden change of behavior after some critical point. Most common example being boiling and freezing of water. In this case, the magnetization of metal as it is been cooled down. Next section, we will derive the critical temperature for the magnetization of the Ising model using statistical mechanics.

Phase Transition

In this section, we want to study the magnetization at equilibrium for any given inverse temperature $\beta$ for a electron sea in a 2D square lattice. In other word, all electrons have 4 neighbors: up, down, left, right. And the distance to all the neighbors are equal.

No interaction model

Turns out, it is useful to first study the Ising model with no interaction ($J_{m,n} = 0$) between the electrons, which can be used for simplifying the full Ising model.

Assume a uniform external field $h_n = h$ for every electron site. The energy of the no interaction model is

$$ E(X_i) = h\sum_n x_n $$

We will rewrite this in terms of the total numbers of up-spin $k$, correspondingly, the total number of down-spin is $(N-k)$

$$ E(X_i) = h\left(\sum_{\text{up}}1 + \sum_{\text{down}}(-1)\right) = h(k - (N-k)) $$

The partition function is

$$ Z= \sum_{k} \textrm{C}^N_k e^{-\beta h(k - (N-k))} = \sum_{k}\textrm{C}^N_k e^{-\beta hk}e^{\beta h(N-k)} $$

where $\textrm{C}^N_k$ is the choose function since there are $N$ choose $k$ amount of states that has the same $k$ number of up-spin.

Substitute $a = e^{-\beta h}$ and $b = e^{\beta h}$

$$ Z = \sum_{k}\textrm{C}^N_k a^kb^{(N-k)} $$

this is the binomial expansion

$$ (a+b)^N = \sum_{k}\textrm{C}^N_k a^kb^{(N-k)} $$

Substitute $a$ and $b$ back, we get

$$ Z = (a+b)^N = ( e^{-\beta h}+e^{\beta h})^N = (2\cosh(\beta h))^N = 2^N(\cosh(\beta h))^N$$

The average energy of the system can be derived from the partition function

$$ \begin{aligned} E &= -\frac{\partial}{\partial \beta}\ln(Z) \\ &= -\frac{\partial}{\partial \beta} (N\ln 2 + N\ln (\cosh(\beta h)) )\\ &= -Nh\frac{\sinh\beta h}{\cosh \beta h} \\ &= -Nh\tanh(\beta h) \end{aligned} $$

The magnetization is related to the energy by

$$ E = NMh $$

Therefore

$$ M = -\tanh(\beta h) $$

when $\beta$ is $0$ (infinite temperature), the magnetization is $0$; the magnetization approaches $-1$ as $\beta$ goes up

2D-Lattice Model

Now, we introduce the local interaction $J_{m,n}$ back and remove the external field $h$ instead. Assume all the interactions have the same strength $J_{m,n} = J$. The energy equation is then

$$ E = -J \sum_{(m,n)} x_nx_m $$

Focus on the local energy associated with the spin of a single electron with its neighborhood $\mathcal{N}$

$$ E_n = -J x_n\sum_{m \in \mathcal{N}} x_m $$

the average of the neighborhood's state is approximately equal to the magnetization at equilibrium, $M$, this is called the mean field approximation

$$ E_n = -J x_n\sum_{m \in \mathcal{N}} x_m \approx -J x_n |\mathcal{N}| M $$

where $|\mathcal{N}$ is the number of neighbors. The bigger the neighborhood size is, the more accurate the approximation gets. Since $M$ and $ |\mathcal{N}|$ are constants, we can treat them as the external field for $x_n$

$$ E_n \approx (-J |\mathcal{N}| M) x_n. $$

From the no-interaction model, we know the expected spin, which is the magnetization, from only an external field is

$$ \mathbb{E}(x_n) = \tanh(\beta J |\mathcal{N}| M) $$

By symmetry, any electron is the same as electron $n$, then we can conclude that the expected magnetization of electron $n$ is just the global magnetization.

$$ M = \tanh(\beta J |\mathcal{N}| M) $$

To find the solution of $M$, we can plot the left hand side and right hand side and find the intersection

With $\beta J |\mathcal{N}| = 0.5$, the only solution for $M$ is $0$. Observe what happen when we start increasing $\beta$

When $\beta$ start increasing, the $\tanh$ function gets steeper and steeper. At some point, it gets so steep that it crosses the $y=M$ line, and three solutions exist. One close to $-1$, one close to $1$ and one that is $0$. That is, the system can be all electrons having a down-spin; all electrons having a up-spin; or exactly 50/50 up and down. The 50/50 solution is unstable because if one electron flip the spin due to random noise, the system will tip toward the spin direction that is the majority.

This cross-over point happens when $\tanh$ is tangent to $y=M$ at $0$

$$ \left.\left(\frac{\partial}{\partial M} \tanh(\beta J |\mathcal{N}| M)\right) \right|_{M=0} = \beta J |\mathcal{N}| $$

Set the slope at $M=0$ to be equal to the slope of $y=M$, which is $1$ to get the critical temperature $\beta_\text{c}$

$$ \begin{aligned} \beta_\text{c} J |\mathcal{N}| & = 1 \end{aligned} \\ \beta_\text{c} =\frac{1}{J |\mathcal{N}| } $$

Ising Model

Contents