PY 502, Computational Physics (Fall 2013)Monte Carlo simulations in classical statistical physicsAnders W. Sandvik, Department of Physics, Boston University1 IntroductionMonte Carlo simulation is a very important class of stochastic methods for calculating thermalproperties of many-particle systems—arguably these are the most important numerical techniquesin statistical physics. Monte Carlo simulation methods are related to the elementary Monte Carlointegration methods that we discussed earlier, but are based on more efficient non-uniform samplingschemes. By using importance sampling, the configurations (particle positions, spin directions, etc.)of a finite but large many-body system (up to millions of degrees of freedom) can be generatedaccording to the Bolzmann distribution, so that thermal expectation values are obtained as simplearithmetic averages of functions “measured” on the configurations.As a simple illustration of the advantages of non-uniform Monte Carlo sampling, consider a one-dimensional integral similar to a thermal expectation value in statistical physics (the discussionhere can be directly generalized to multi-dimensional integrals);hAi =ZL−LP (x)A(x)dx,ZL−LP (x)dx = 1, (1)where P (x) is an arbitrary probability distribution. By randomly sampling M points x1, . . . , xMin the range [−L, L], the expectation value is estimated ashAi ≈2LMMXi=1P (xi)A(xi). (2)As we discussed before, if P (x) is sharply peaked in a small region, the statistical fluctuations of thisestimate will be large as only a small fraction of the generated points will fall within the dominantregion. If we instead sample the points according to some probability distribution W (x), i.e., theprobability of picking a point in an infinitesimal range [x, x + dx] is W (x)dx (we now assume thatthis can be done for arbitrary W (x), and leave for later discussion how this is accomplished inpractice), then the estimate for the expectation value ishAi ≈1MMXi=1P (x)W (x)A(xi). (3)This has less statistical fluctuations than the estimate (2) of the uniform sampling if W (x) is peakedin the same region as P (x) and if the function A(x) is well-behaved, in the sense of being reasonablysmooth and not very small where P (x) is large and vice versa. More precisely, the fluctuation inthe values sampled using the distribution W (x) is given byσ2W[A] =ZLLP (x)W (x)A(x) − hAi2W (x)dx, (4)1which in principle can be minimized by choosing a particular W (x). In general (for the multi-dimensional integrals or sums encountered in statistical physics) it is not possible in practice tofind the optimal W (x) that minimizes the fluctuations, but if P (x) has much larger variations thanA(x) a very good solution is to use W (x) = P (x). The expectation value is then just the simplearithmetic average of A(x) over the sampled configurationshAi ≈1MMXi=1A(xi), (5)and the expected fluctuation of the measured values isσ2P[A] =ZLL[A(x) − hAi]2P (x)dx. (6)It should be noted that the distribution of the values Aiis typically not Gaussian, and hence tocalculate the statistical errors of hAi estimated as (5) the values should first be binned, in the sameway as we discussed previously in the chapter on Monte Carlo integration.In statistical physics, P is a sharply peaked exponential function e−E/kBTof the energy and A is typ-ically a linear or low-order polynomial function of the system degrees of freedom. The fluctuationsin P are thus very large relative to those of A and the sampling using P as the probability distri-bution is then close to optimal. This is what is normally meant by the term importance sampling.Using importance sampling instead of uniform random sampling is crucial when a small fractionof the configuration space dominates the partition function, which is always the case with theBolzmann probability in statistical mechanics models at temperatures of interest. How to achievethe correct distribution in practice is the main theme of this chapter; we will discuss importancesampling schemes for both lattice and continuous-space models.One of the primary utilities of Monte Carlo simulation is in studies of phase transitions and criticalphenomena. This will be the focus of applications discussed here. Although there are analogoussimulations methods available also for quantum systems (called quantum Monte Carlo methods),we will here consider only Monte Carlo simulations of classical many-body models. In addition todescribing simulation algorithms, we will also discuss how simulation data is analyzed in order tolocate phase transitions and extract critical exponents.In the following sections we will first briefly review the the expressions for thermal expectation valuesin systems of particles identified by coordinates in continuous space. We then consider models withdiscrete degrees of freedom on a lattice, focusing on spin models, the Ising model in particular.We will discuss the general form of the detailed balance condition that can be used to sampleconfigurations according to any desired probability distribution. Monte Carlo simulation algorithmsare for instructional purposes often developed in the context of the Ising model, and we will followthis path here as well (it should also be noted that Ising models are of continued importancein research). We will develop a standard program for simulations using the Metropolis algorithm,which is based on evolving (updating) configurations by flipping individual spins. We will study theperformance of this method using autocorrelation functions, which characterize the way in whichgenerated configurations gradually become statistically independent of the past configurations.This leads us to the problem of critical slowing down, which makes accurate studies close to phasetransitions difficult. Critical slowing down can be often be greatly reduced, in some cases completely2eliminated, using cluster algorithms, where clusters of spins are flipped collectively (using cluster-building rules that satisfy detailed balance). We will develop a cluster Monte Carlo program and useit to study the ferromagnetic phase transition in the two-dimensional Ising model. For this purpose,we will also discuss finite size scaling methods used to extract critical points and exponents. Finally,we will return to problems involving particles in continuous space; we will develop a program forsimulating a mono-atomic gas and its phase transition into the liquid state.2 Statistical