The Monte Carlo method introduced here with the example of the Ising model is very similar in spirit to the methods already encountered in the last lectures for the cases of random walks, the modelling of diffusion by them, or percolation, by using random numbers for the description of a system. However, here there is the added description of interactions through random numbers, triggering a dynamical evolution of the system driven by interactions. In order to use random numbers for their simulation also, it is a sine qua non to formulate the interactions in a probabilistic fashion. This then also accounts for the transition in language from using random methods for the description of a system to a proper Monte Carlo simulation.
Our goal is to simulate the system of spins - for definiteness now arranged in a two-dimensional square lattice - and their interactions. These interactions are facilitated by the system being in thermal equilibrium with a heat bath, allowing energy to shuffle back and forth between the spins and the heat bath. As the energy is exchanged, spins are flipped, causing the system to move from one micro-state to another, each of which corresponds to a particular arrangement of the spins. As already stressed before, macroscopic quantities such as the total magnetisation of the system then depend on the probabilities of the system being in the various micro-states contributing - these macroscopic observables can be calculated through a thermal (or time) average.
The Monte Carlo method uses a stochastic method to describe the exchange of energies with the heat bath, i.e. the flipping of spins. The procedure is as follows:
First of all, the system is initialised with a specific spin combination, like e.g. all spins pointing into the plus-direction, or all spin directions chosen one by one at random.
The interaction is then modelled by choosing a particular spin, say si, and calculating the energy Eflip required to flip it. This is done, in the case of the Ising model, through the energy equation above,
evaluating the energy before and after the flip and taking the difference. If Eflip is negative, i.e. if flipping the spin would lower the energy of the system, the spin is flipped. If it is positive the decision of whether to flip the spin or not is made on the basis of the Boltzmann factor, exp[-Eflip/(kBT)]. If this factor is larger than a random number, the spin is flipped, otherwise the system remains unaltered. Hence, in the case Eflip is positive, the system may or may not move to another micro-state, whereas in the case of Eflip being negative it certainly will move into a different micro-state. This describes one Monte Carlo step, and the algorithm outlined here is known as the Metropolis algorithm.
The procedure outlined above is repeated a huge number of times such that every spin has many chances to flip. We can interpret each time step in doing this as one interaction with the surrounding heat bath, the effect of this interaction clearly depends on the temperature, entering the Boltzmann factor. For T→0 this factor becomes unity for positive Eflip, meaning that spin flips are allowed only if they reduce the energy of the system. In contrast, for large T, the Boltzmann factor becomes less and less suppressing, allowing the spins to flip more easily, and the system thus to fluctuate more and more. In fact, in the limit T→∞, the system enters the paramagnetic phase.
In order to determine macroscopic observables, one would typically wait for a number of complete sweeps through the system, to allow for some equilibration. Then, the observables would emerge as averages over a sufficiently large number of sweeps through the lattice to damp down fluctuations.
These arguments show that the Monte Carlo algorithm will lead to the qualitatively correct behaviour. The natural question therefore is, whether it is also formally and quantitatively correct. To see this consider one spin flip, connecting two micro-states labelled with 1 and 2 and with energies E1 and E2, respectively. The rate of transitions W(a→ b) between the two states is given by the probabilities of the transitions. If, for instance, E1 > E2, then W(1→2)=1. Likewise, if the system is in state 2, the rate of the transition to state 1 is given by the Boltzmann factor and hence
This is because in this case Eflip = E1-E2 > 0. When our system is in equilibrium, the probability of finding it in a particular micro-state is independent of time (the heat bath is at constant temperature), therefore the number of transitions from state 1 to 2 must be equal to those in reverse direction. Denoting with P1 and P2 the probabilities of the system to be in state 1 or 2, respectively, therefore translates into
at equilibrium. Inserting the expressions for the rates thus yields
This is in perfect agreement with the thermodynamics of a system in thermal equilibrium: Boltzmann factors determine the probabilities for a system to be in a particular micro-state, with the energies of that micro-state entering. The Metropolis algorithm therefore lends itself to a description of the system in agreement with thermodynamics, since ultimately it will allow each micro-state to be populated with the thermodynamically correct probability. Care must only be taken to ensure that there is enough time for the system to explore all micro- states. In our case this is satisfied as long as enough MC sweeps are performed, since all micro-states of the system can be obtained with a finite, albeit potentially large number, of spin flips.
Before we move on to describe some first findings with this method, it is worth stressing that in principle all thermodynamics systems in thermal equilibrium can be studied with the Metropolis algorithm. The only things necessary are an enumeration of all micro-states in a useful way and a determination of transition probabilities between them in accordance with the laws of statistical physics. In that way, the Metropolis algorithm simulates the canonical ensemble (by building the phase space integral over the allowed states). An extension to the micro-canonical or the grand canonical ensemble are also feasible.