Classical Monte Carlo Lab

Samples instead of a sum

The thermal average of an observable is

where the configurations \(\mathbf S_t\) are drawn with \(P\propto e^{-\beta H}\). Uniform random sampling does not work: the important configurations are an exponentially small fraction of phase space.

The error decreases as \(M^{-1/2}\) independently of the system size, which makes the method powerful for high-dimensional sums. However, \(Z\) itself (and hence the free energy and entropy) is not measured directly.

Detailed balance

If the transition probability \(W\) for \(\mathbf S\to\mathbf S'\) satisfies detailed balance and the chain is ergodic, its stationary distribution is the Boltzmann distribution:

Metropolis: a change is proposed and accepted with probability \(\min(1,e^{-\beta\Delta E})\). Heat bath: the chosen spin is redrawn directly from the Boltzmann distribution in the local field of its neighbours.

One sweep is \(N\) update attempts. In Metropolis the sites are chosen at random. Strictly sequential Metropolis can get stuck in the wrong distribution in 1D because moves with \(\Delta E=0\) are always accepted; this version uses random site selection.

Forgetting the start

The chain starts from a random or ordered state and reaches equilibrium only after a few \(\tau_{\rm exp}\). Samples from this period are discarded (thermalization). When a parameter changes, the measurement is reset automatically and as many samples as given in the "Thermalization" box are discarded.

A practical test: two runs with the same parameters from a random (hot) and an ordered (cold) start must converge to the same average. At first-order transitions and in frustrated systems this test can fail: the chain stays in a metastable state.

Effective number of samples

Consecutive Markov samples are not independent. The variance of the mean is

\(M\) samples are equivalent to only \(M/2\tau_{\rm int}\) independent ones. \(\tau_{\rm int}\) depends on the observable: \(|m|\) is usually slower than \(e\).

\(\tau_{\rm int}\) in the readout is measured in sweeps. The autocorrelation view in Figure 1 computes the same quantity directly from \(\rho(t)\); agreement between the two estimates is a sign of a reliable measurement.

How the errors are computed

The series is repeatedly averaged in pairs; once the block length exceeds \(\tau_{\rm int}\) the block means become independent and the error reaches a plateau. For nonlinear functions of averages such as \(\chi\), \(C\) and \(U\) the error is obtained from a 32-block jackknife:

If the measurement is shorter than a few hundred \(\tau_{\rm int}\), both \(\tau_{\rm int}\) and the error come out too small. This is common for Metropolis at the critical point; trust the agreement of two independent runs rather than the error bar.

A fair comparison of algorithms

Comparing algorithms by \(\tau_{\rm int}\) per sweep is misleading: a Wolff sweep is more expensive than a Metropolis sweep. The right measure is the number of independent samples per second:

At the critical point \(\tau\sim L^{z}\): for the 2D Ising model \(z\simeq2.17\) for Metropolis and \(z\lesssim0.3\) for Swendsen–Wang and Wolff. For large \(L\) the difference therefore reaches hundreds. Away from the critical point local algorithms can be more efficient.

Swendsen–Wang

The Ising/Potts partition function can be rewritten with bond variables. Every satisfied bond (equal spins for FM) is activated with probability

and each cluster connected by active bonds is flipped independently (assigned a random state for Potts). The acceptance probability is 1.

The clusters carry the physical correlations: at \(T_c\) they are at the percolation threshold. Large-scale changes are therefore made in a single step.

Single cluster

A single FK cluster is grown from a random seed site and flipped. Because large clusters are chosen more often it is somewhat more efficient than SW. In the disordered phase the mean cluster size equals the susceptibility: \(\langle|C|\rangle=\chi'/\beta\).

A hidden trap: the rule "generate clusters until the number of flipped sites exceeds \(N\), then measure" correlates the measurement with the size of the last cluster and produces a measurable bias. Here the number of clusters per sweep is set during thermalization and then frozen.

Embedding and overrelaxation

For XY and Heisenberg a random \(\hat r\) is chosen; the reflection \(\mathbf S\to\mathbf S-2(\mathbf S\cdot\hat r)\hat r\) defines an embedded Ising model that flips the sign of \(\hat r\cdot\mathbf S\):

With anisotropy only reflections that leave \(H\) invariant (\(\hat r\in xy\) or \(\hat r=\hat z\)) can be used. They map the \(xy\) plane onto itself, so the chain is not ergodic; a Metropolis sweep is added to every sweep. Overrelaxation reflects the spin about the local field, \(\mathbf S\to2(\mathbf S\cdot\hat{\mathbf H})\hat{\mathbf H}-\mathbf S\); it conserves energy and is cheap, but on its own it is microcanonical.

A finite system has no transition

At finite \(L\) every observable is analytic; singularities are rounded once \(\xi\sim L\). At the critical point

2D Ising: \(\beta=1/8,\ \gamma=7/4,\ \nu=1\), \(C\) diverges logarithmically. 3D Ising: \(\nu\simeq0.630\), \(\gamma\simeq1.237\). The positions of the \(\chi\) and \(C\) peaks shift with \(L\); reading \(T_c\) off a single \(L\) gives a systematic error.

The crossing point

For an order parameter with \(n\) components

In the ordered phase \(U\to1\), in the disordered phase (Gaussian fluctuations) \(U\to0\). Being dimensionless under scaling, the curves for different \(L\) cross at a single point at \(T_c\), up to leading corrections.

For Ising this is \(3/2\) times the standard definition \(1-\langle m^4\rangle/3\langle m^2\rangle^2\); the crossing value \(U^*\) is universal but depends on boundary conditions and aspect ratio (2D square Ising: \(U^*\simeq0.916\)). No crossing (1D, 2D Heisenberg) means no long-range order. At a first-order transition \(U\) develops a negative dip.

A curve from one run

Samples collected at \(\beta_0\) also give the average at a nearby \(\beta\) (Ferrenberg–Swendsen):

The method relies on the tails of the \(P(E)\) histogram; the reliable range is of order \(|\Delta T|\lesssim T^2/\sqrt{NC}\) and shrinks as \(L\) grows. The dashed curve in Figure 2 is drawn only where the effective number of independent samples \(N_{\rm eff}=(\sum w)^2/\sum w^2/2\tau\ge30\).

No continuous symmetry breaking in 2D

In 2D XY and Heisenberg models with short-range interactions \(\langle\mathbf m\rangle=0\) at \(T>0\). Spin-wave fluctuations diverge logarithmically, like \(\int d^2k/k^2\).

On a finite lattice \(|m|\) is nonzero and decreases slowly with \(L\). In 2D Heisenberg \(\xi\propto e^{2\pi/T}\) easily exceeds the lattice at low \(T\); the state that looks "ordered" on screen is a finite-size effect. The Binder curves do not cross.

Unbinding of vortex pairs

In the 2D XY model correlations decay algebraically at low \(T\) (quasi-long-range order). At \(T_{\rm BKT}\simeq0.893\) bound vortex–antivortex pairs unbind and the helicity modulus makes a universal jump:

At finite \(L\), \(\Upsilon\) decreases continuously and crosses the \(2T/\pi\) line slightly above \(T_{\rm BKT}\); the corrections decay very slowly, like \(1/\ln L\). The Binder crossing is not a good criterion for BKT: the curves merge below \(T_{\rm BKT}\).

The universality class depends on \(n\)

\(T_c\) on the simple cubic lattice: Ising \(4.5115\), XY \(2.2018\), Heisenberg \(1.4430\). In the Heisenberg model the effective anisotropy \(A=(1-\lambda)-K\) changes the number of order-parameter components:

A small anisotropy produces a crossover near \(T_c\): \(n=3\) behaviour at large distances from \(T_c\), \(n=1\) or \(2\) close to it. In 2D an easy axis restores an Ising transition and an easy plane a BKT transition. The page computes \(U\) with the \(n\) appropriate for the chosen anisotropy.

From continuous to first order

On the 2D square lattice self-duality gives \(T_c=1/\ln(1+\sqrt q)\). For \(q\le4\) the transition is continuous (\(q=2\) Ising, \(q=3\): \(\nu=5/6\), \(q=4\): \(\nu=2/3\)); for \(q>4\) it is first order: ordered and disordered phases coexist at \(T_c\), \(P(e)\) is double peaked and there is latent heat.

At \(q=5\) the transition is very weakly first order; the correlation length at \(T_c\) is \(\sim2500\) lattice spacings. At the sizes accessible in a browser it cannot be distinguished from a continuous transition. Choose \(q\ge8\) to see the double peak.

AFM Ising on the triangular lattice

The three bonds of a triangle cannot all be satisfied. In the ground state every triangle has exactly one unsatisfied bond: \(e_0=-|J|\), but the number of such states grows exponentially. The Wannier entropy is

There is no transition at \(T>0\); correlations become algebraic as \(T\to0\). The page shows the three-sublattice \(\sqrt3\times\sqrt3\) order parameter and, as "defects", the fraction of triangles whose three spins are equal. The lattice side must be a multiple of 3.

Slow cooling and defects

In simulated annealing \(T\) is lowered over time. A system cooled slowly enough approaches the ground state; a system cooled quickly freezes in domain walls, vortices or frustration defects. The remaining defect density decreases with the cooling rate.

Cluster algorithms fail in frustrated systems: at low \(T\) the FK clusters cover the whole lattice and the defects cannot be removed. Compare Wolff with Metropolis for the triangular AFM. This is an example of an algorithm that is correct but effectively frozen.