The Rutgers Scholar
An Electronic Bulletin of Undergraduate Research

Monte Carlo simulations of small model systems
with phase transitions

Stephanie Devlin* and E. Roger Cowley
Department of Physics, Rutgers University
Camden, New Jersey 08102

*Rutgers Undergraduate Research Fellow


We have studied the ability of a small, idealized system of magnets to predict accurately the properties of an infinitely large system at and around a phase transition. Specifically, by the Ising model, we have investigated the second order phase transition of a ferromagnetic system into a paramagnetic system.


A system of many interacting particles can exhibit a phase transition: that is, exhibit a change in a feature that characterizes the system. We have investigated the second order phase transition of a ferromagnetic system into a paramagnetic system.

A ferromagnet exhibits a permanent macroscopic magnetization whereas a paramagnet exhibits some magnetic properties, but not a macroscopic magnetization. At a certain temperature, the critical temperature or critical point, a ferromagnetic system changes from one displaying macroscopic magnetization into a seemingly unmagnetized system – a paramagnetic system.

This research is focused on the behavior of a small system in the critical temperature region. With suitable approximations it is the case that predictions made from such small sytems can be used to predict the properties of an infinite system in that same region.

Methods and background

This research needed a simple model with magnetic spins to show magnetic phase transitions. We chose one called the Ising Model (Figure 1).

Figure 1. Pictorial representation of the Ising model.

The Ising Model assumes that each spin is able to point in the +z direction or the –z direction and that the ith spin in the system has the value of ±si. There are interactive forces between spins. The interaction is strongest between nearest neighboring spins, and the Ising Model assumes that the interaction between nearest neighbors is the only interaction in the system. It neglects the forces associated with spins further away. The energy of the system, with no external field, is then:

where <nn> shows that the summing is over the nearest neighboring spins and J is the exchange constant. When the spins are parallel the energy is negative and the system attains greater stability. A system with every spin parallel to each other spin is a ferromagnetic system. The alignment of the spins creates a macroscopic magnetization. The probability of finding the model in any particular state for the Ising Model is proportional to the Boltzmann factor:

where E is the energy defined above, k is the Boltzmann constant, and T is the temperature.

Figure1 shows a small Ising Model lattice. Each arrow represents an individual spin and its associated direction. The central, 4 by 4 structure whose properties are of interest lie inside the lines. In so small a lattice, the objects at the edges are a source of concern because they have fewer nearest neighbors than do the objects in the center. Consideration of the lines and the spins outside of the lines make it possible to impose what are called periodic boundary conditions and thus obtain results more like those that would hold for a larger and more realistic system. In implementing periodic boundary conditions, we treat the left side of the lattice as if it were touching the right side -- wrapping the lattice around so opposite sides touch. The lines of spins on the outside of the lattice are purposely chosen to be identical to the spins on the opposite side of the 4 by 4 lattice to make the wrapping possible. Figure 1 also illustrates the energies associated with neighboring spins. The energy equals -J when the spins are parallel and equals +J when the spins are anti-parallel.

The Ising Model is an extremely simple model of a magnetic system, which, despite its simplifications, captures the essential physics of magnetism. The numerical approach used to analyze the Ising Model is called the Monte Carlo Method.

The system is thought of as if it were in a heat bath. As the system gains and loses energy from and to the heat bath, spins are flipped causing a move to a particular microstate. Figure 2, from Yeomans “Statistical Mechanics of Phase Transitions” shows how the Monte Carlo method works:

Figure 2. Monte Carlo flow chart.

In the Monte Carlo Method, first a lattice structure is set up. The initial values for the individual spins and for the external field H are defined. One element in the lattice is chosen, its spin changed from + to - or from - to plus (flipped), and then the energy of the system with the new flipped spin is calculated. The element is left in its new position if the energy of the system is lowered by the flip (if the change in energy is negative). If the change in energy is positive, then the flip may or may not be retained. To decide whether to retain the flip we calculate a random number between 0 and 1. Next we calcultate from the energy of the system the probability (from the Boltzmann factor) that the system will attain that state. If that probability is greater than the random number, the spin flip is retained; if the probability is less than the random number, the spin is not flipped. The comparison to the random number allows the magnetic system to enter states of higher energy than that of a ferromagnet. These steps complete one Monte Carlo time step. Another element in the lattice is then chosen, and the process repeated for that element. This procedure is repeated a large number of times, so that each spin is given many chances to flip. The variables - energy, magnetization, susceptibility, heat capacity and absolute magnetization - are calculated, as well as the averages of those variables. The system’s dependence on temperature is carried by the temperature's appearance in the probability term.

Results, discussion, and conclusions

Plot A displays the magnetization against the temperature. It shows how the combination of the Ising Model and the Monte Carlo Method does detect and calculate a second order phase transition in a small model system, 20 by 20 in this case.

Figure 3.Spin (20,20) Field = 0.0.

The 20 by 20 system starts out at low temperatures by displaying ferromagnetic properties, with the magnetization approximately equal to one. As the temperature is increased, the spins tend to flip, and eventually the system no longer exhibits the macroscopic magnetization it had originally. This is the phase transition.

Since the research centered about small system phase transitions, calculations were mostly on a 20 by 20 lattice and were concentrated around the critical temperature. A few other systems were studied for comparisons, however. The critical temperature for an infinite ferromagnetic to paramagnetic system studied through the Ising Model is approximately 2.269. The obstacle encountered with smaller sized systems is a dulling and broadening of the critical temperature region. Figure 4 and Figure 5 show, respectively, the temperature dependence of the heat capacity for an 80 by 80 lattice and for a 10 by 10 lattice.

Figure 4. Spin(80,80) Field = 0.0

Figure 5. Spin (10,10) Field = 0.0

For an infinite system, the heat capacity peak is finite at the critical temperature. One can see that for the 80 by 80 system, the peak is still a sharp one; but quickly the peak dulls and the critical temperature region broadens for the 10 by 10 system. Comparing these two graphs, one can see that the 10 by 10 system, which is represented by the triangular points, has a broader and more rounded peak at the critical temperature than does the 80 by 80 system. This effect of the decreasing lattice size on the critical temperature region makes calculation of the critical temperature more difficult.

There are two main regions to look at for a system with a phase transition – above the critical temperature and below the critical temperature. Figure 6 displays the field dependence of the magnetization for a 20 by 20 lattice at a temperature of 2.5, which is above the critical temperature.

Figure 6. Spin (20,20) Temperature = 2.5

The magnetization shows no discontinuity. The plot starts off linearly. The slope of that line is the susceptibility of the system. Due to the fact that the magnetization cannot be greater than one, the line gradually saturates to that limit.

Figure 7 shows magnetization versus magnetic field strength in the region below the critical temperature for the 20 by 20 and the 40 by 40 lattices. The circular points represent the 20 by 20 lattice. The triangular points represent the 40 by 40 lattice, and the data points where the two systems coincide are the squared points. This plot shows how the two systems’ data differ at low fields, and then coincide at larger fields – mimicking the infinite system.

Figure 7. Temperature = 2.2; Circles for Spin (20,20); Triangles for Spin (40,40).

Plot E shows the field dependence of the susceptibility.

Figure 8. Temperature = 2.1; Circles for Spin (20,20); Traingels for Spin (40,40).

The circular points represent the data for the 20 by 20 lattice and the triangles represent the 40 by 40 lattice data. Again in low fields the two systems show different values for the susceptibility, but as we increase the field, both systems’ data coincide within their deviations. This convergence illustrates that as the field is increased the small model systems can mimic the properties of an infinite system.

We conclude that there are two important regions to consider when dealing with a finite sample, as shown in Figure 9.

Figure 9. Findings.

Within Region I, the finite system and infinite systems behave similarly, but as the temperature and field decrease -- this is in Region II -- the finite system does not act in the same way as the infinite one. It is unfortunate, however, that the critical point and the phase transition lie within this region.

Further study of the magnetization and the extrapolation from the approximate infinite system data, the squared points, into the region of lower field would allow improved predictions of the infinite system within the region of the phase transition.


Giordano, Nicholas J. 1997. Computational Physics. Prentice-Hall, Inc., Upper Saddle River.

Yeomans, J. M. 1992. Statistical Mechanics of Phase Transitions. Oxford University Press, New York.

Copyright 2000 by Stephanie Devlin and E. Roger Cowley
Current URL: