Generating Ising model steady state configurations

Question

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?

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.