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.

enter image description here

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

  • 3
    $\begingroup$ What exactly are the axes of your graph? $\endgroup$
    – Philip
    Jun 23, 2020 at 0:52
  • $\begingroup$ y label counts normalized x label energy per spin $\endgroup$
    – Simo7
    Jun 23, 2020 at 23:46

1 Answer 1


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):

enter image description here

(The computation using other boundary conditions are similar.)

  • $\begingroup$ it helps me so much thanks $\endgroup$
    – Simo7
    Jul 1, 2020 at 19:53
  • $\begingroup$ why the division by 2 in the binomial coefficient? $\endgroup$
    – Simo7
    Jul 1, 2020 at 20:52
  • $\begingroup$ 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$. $\endgroup$ Jul 2, 2020 at 7:42
  • $\begingroup$ What is $D(\sigma)$? $\endgroup$
    – AlphaOmega
    Feb 21, 2021 at 10:19
  • $\begingroup$ @AlphaOmega The number of pairs $\sigma_i, \sigma_{i+1}$ with $\sigma_i\neq\sigma_{i+1}$. $\endgroup$ Feb 21, 2021 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.