What is the most efficient way to simulate steady state configurations of the Ising model? I am just interested in having a large set of random steady state configurations of the 1D Ising model (with homogeneous coupling constants). A few ideas came to mind:

  1. Brute force sampling. Since the Ising model is exactly solvable in 1D and 2D, one has exact expressions for the probabilities of each state. However, random sampling over a set of $2^N$ will likely cause memory problems for small $N$ already.
  2. Monte Carlo dynamics. One could run the usual Monte Carlo algorithms (e.g. Glauber dynamics) on random initial states and wait until the system converges to thermal equilibrium. However, this seems inefficient when you are not interested in the dynamics and only want steady state configurations.
  3. Using the density of states. One could also first randomly sample the energy of the system, according to $P(E) \sim N(E) \exp(-\beta E)$, where $N(E)$ is the density of states, which is computable (at least numerically). Then one generates a random configuration with this energy, e.g. using a spin flip algorithm where one flips single spins to increase/decrease the energy until it matches the target energy. But I'm not sure if the configurations obtained this way statistically follow the Boltzmann distribution.

Note: in 1D there is also an exact expression for the Ising density of states, $g(E(k)) = 2 \binom{N-1}{k}$ with $E(k) = -N + 2k + 1$. See this other question: Ising model density of states.

Any ideas on what is the best way to approach this?


1 Answer 1


For the one-dimensional model, the most efficient way, by far, of simulating the Ising model is by using a Markov chain on $\{-1,1\}$, generating one spin at a time, conditionally on the values taken by the previous spins. Note also that in this way, you are sampling exactly from the Gibbs distribution, with no approximation (in contrast to a Monte Carlo approach).

For simplicity, let me consider the model with free boundary condition, that is, the model with Hamiltonian $$ \beta\mathcal{H} = - \beta\sum_{i=2}^N \sigma_{i-1}\sigma_i . $$ (You can also add a magnetic field, but I won't do it here to simplify the exposition). Then, $\sigma_1$ is equal to $+1$ or $-1$ with probability $\tfrac12$ by symmetry. Moreover, for any $k\geq 2$, $$ \mathrm{Prob}(\sigma_k=\sigma_{k-1} \,|\, \sigma_1, \dots, \sigma_{k-1}) = \mathrm{Prob}(\sigma_k=\sigma_{k-1}) = \frac{e^{\beta}}{e^{\beta} + e^{-\beta}} = \frac{1}{1+e^{-2\beta}}. $$ Let us call this probability $p$.

To summarize:

  • You sample $\sigma_1$: it is $+1$ with probability $\tfrac12$ and $-1$ with probability $\tfrac12$.
  • Given $\sigma_1$, you set $\sigma_2 = \sigma_1$ with probability $p$ and $\sigma_2 = -\sigma_1$ with probability $1-p$.
  • Given $\sigma_2$, you set $\sigma_3 = \sigma_2$ with probability $p$ and $\sigma_3 = -\sigma_2$ with probability $1-p$.
  • and so on...

This is very easy to implement and extremely fast (of course, compute $p=1/(1+e^{-2\beta})$ only once). Then most of the time is taken by the pseudorandom number generation. In this way, you can simulate chains of arbitrarily large length without any problem.

(See also this answer for another point of view of the relationship between one-dimensional models and Markov chains.)

Explanation of the formula for $p$.

The simplest way of seeing why the formula for $p$ given above holds is by using either the random-cluster or the high-temperature representations of the Ising model, if you're familiar with them (they are described, for instance, in Sections 3.7.3 and 3.10.6 in this book).

If you are not familiar with these representations, let me try to provide a direct argument.

Let $s_1,\dots,s_N \in \{-1,1\}$ and write $s=(s_1,\dots,s_{k-1},s_k,\dots,s_N)$ and $s'=(s_1,\dots,s_{k-1},-s_k,\dots,-s_N)$ (that is, the configuration $s'$ is obtained from the configuration $s$ by flipping the spins at $k, k+1, \dots N$).

Now, $$ \frac{{\rm Prob}(\sigma = s)}{{\rm Prob}(\sigma = s')} = \frac{\exp\bigl( -\beta \mathcal{H}(s) \bigr)}{\exp\bigl( -\beta\mathcal{H}(s') \bigr)} = \exp(2\beta\, s_{k-1}s_{k}). $$ In particular, $$ \frac{{\rm Prob}(\sigma_k=\sigma_{k-1})}{{\rm Prob}(\sigma_k = -\sigma_{k-1})} = \exp(2\beta). $$ But this implies that $$ {\rm Prob}(\sigma_k=\sigma_{k-1}) = e^{2\beta}\, {\rm Prob}(\sigma_k = -\sigma_{k-1}) = e^{2\beta} \bigl( 1 - {\rm Prob}(\sigma_k = \sigma_{k-1}) \bigr), $$ and therefore $$ (1+e^{2\beta})\, {\rm Prob}(\sigma_k=\sigma_{k-1}) = e^{2\beta}, $$ from which the formula for $p$ follows immediately.

  • $\begingroup$ Thanks, this is definitely what I was looking for! I was just wondering how to derive the conditional probability that you wrote down. It seems very intuitive, but somehow I don't see how to derive it. I'm wondering whether it is straightforward to see this or whether it follows from the Markov chain properties you refer to in your other answer (which I don't fully understand). $\endgroup$ Sep 21, 2020 at 14:19
  • 1
    $\begingroup$ I added a direct derivation. The idea is extremely simple: if you flip the spin $\sigma_k$ and all the spins at its right, then you change the energy by $\pm 2$, depending on whether the spin $\sigma_k$ agrees with $\sigma_{k-1}$ before or after the flip. In any case the ratio of the probabilities will be given by $e^{\pm 2\beta}$, from which the formula for $p$ follows easily. More details are given in the answer. $\endgroup$ Sep 21, 2020 at 15:11
  • $\begingroup$ I see the point, and also think this makes sense. However, in the derivation you assume that for the state $s'$, $s_k, \ldots s_N$ are all flipped, which is necessary to reduce the first formula to a single term. But I find it confusing, because it seems that we are imposing additional conditions on $s_{k+1}, \ldots s_{N}$ here. $\endgroup$ Sep 21, 2020 at 15:36
  • $\begingroup$ Actually, I'm having some doubts now. From this construction, $P(\sigma_{k+l} = \sigma_{k}) = \prod_{i=1}^{l} P(\sigma_{k+i} = \sigma_{k+i-1}) = p^{l}$. But this would imply that the two-point correlation function $\langle \sigma_{k} \sigma_{k+l} \rangle = 2 P(\sigma_{k+l} = \sigma_{k}) - 1$ would decay as a power-law, correct? But we know that instead it decays exponentially. $\endgroup$ Sep 21, 2020 at 15:59
  • 1
    $\begingroup$ No, we're not adding conditions on $s_{k+1},\dots,s_N$. When computing the probability that $\sigma_k=\sigma_{k-1}$, one has to sum over all admissible configurations $s$ (that is, those such that $s_{k}=s_{k-1}$). To each such configuration corresponds exactly one configuration in which $s_{k}=-s_{k-1}$, obtained by flipping the spins to the right of $s_{k-1}$. Moreover, for each of these pairs of configurations, the ratio of probabilities is the same $e^{2\beta}$. This implies that the ratio of the probabilities that $\sigma_k=\sigma_{k-1}$ and $\sigma_k=\sigma_{k-1}$ is also $e^{2\beta}$. $\endgroup$ Sep 21, 2020 at 16:00

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.