# Boltzmann distribution in Ising model

I've written in Matlab a code for a Ising model in 1 dimension with 40 spins at $$k_{B}T=1$$.

I record the energy of every step in a Metropolis Monte Carlo algorithm, and then I made an histogram like this.

I want to show the theoretical Boltzmann distribution. What is the exact formula to get this shape? $$Ax \exp(-\beta x)$$?

• What exactly are the axes of your graph? Jun 23, 2020 at 0:52
• y label counts normalized x label energy per spin Jun 23, 2020 at 23:46

I have to make a number of assumptions, as you did not state all the necessary information. So, I am going to assume that you are using periodic boundary conditions, that is, your Hamiltonian is $$\mathcal{H}(\sigma) = -\sum_{i=1}^N \sigma_i\sigma_{i+1},$$ where I have denoted by $$N$$ the number of spins (that is, $$N=40$$ in your case) and used the convention that $$\sigma_{N+1}=\sigma_1$$.

Write $$D(\sigma)$$ the number of $$i\in \{1,\dots,N\}$$ such that $$\sigma_i\neq \sigma_{i+1}$$, still using the convention that $$\sigma_{N+1} = \sigma_1$$. Note that $$D(\sigma)$$ is necessarily an even number (because of the periodic boundary conditions).

Then, the total energy can be rewritten as $$\mathcal{H}(\sigma) = D(\sigma) - (N-D(\sigma)) = 2D(\sigma) - N \,.$$ The probability that $$D(\sigma)=\delta N$$ (with $$\delta \in \{0,\frac{2}{N},\frac{4}{N},\dots,\frac{2\lfloor N/2 \rfloor}{N}\}$$) is $$\mathrm{Prob}(D = \delta N) = \binom{N}{\delta N}\frac{\exp(-2\beta\delta N)}{\frac{1}{2}\bigl(1-\exp(-2\beta)\bigr)^N + \frac{1}{2}\bigl(1+\exp(-2\beta)\bigr)^N},$$ since there are $$\binom{N}{\delta N}$$ ways of choosing the $$\delta N$$ pairs of disagreeing neighbors.

This can be easily reformulated in terms of the energy per spin $$\frac{1}{N}\mathcal{H}(\sigma) = \frac{2}{N}D(\sigma) - 1.$$ Note that the possible values of the latter are $$-1$$, $$-1+\frac{4}{N}$$, $$-1+\frac{8}{N}$$, ... , $$-1+\frac{4\lfloor N/2\rfloor}{N}$$.

The probability of observing an energy per spin equal to $$\epsilon$$ is then given by $$\mathrm{Prob}(\mathcal{H}=\epsilon N) = \binom{N}{\frac{1+\epsilon}{2}N} \frac{\exp\bigl(-\beta (1+\epsilon) N\bigr)}{\frac{1}{2}\bigl(1-\exp(-2\beta)\bigr)^N + \frac{1}{2}\bigl(1+\exp(-2\beta)\bigr)^N} .$$

Here is a plot of the distribution for your parameters $$N=40$$ and $$\beta=1$$ (only values of $$\epsilon$$ smaller than $$-0.2$$ are indicated as higher values have too small probability at this temperature):

(The computation using other boundary conditions are similar.)

• it helps me so much thanks Jul 1, 2020 at 19:53
• why the division by 2 in the binomial coefficient? Jul 1, 2020 at 20:52
• Do you mean in $\binom{N}{(1+\epsilon)N/2}$? I have just rewritten the original binomial coefficient $\binom{N}{\delta N}$ using the identity $\epsilon=2\delta - 1$ (which follows from the identity $\mathcal{H} = 2D - N$ by dividing by $N$ and using $D=\delta N$ and $\mathcal{H}=\epsilon N$), which gives $\delta = (1+\epsilon)/2$. Jul 2, 2020 at 7:42
• What is $D(\sigma)$? Feb 21, 2021 at 10:19
• @AlphaOmega The number of pairs $\sigma_i, \sigma_{i+1}$ with $\sigma_i\neq\sigma_{i+1}$. Feb 21, 2021 at 14:15