q-State Potts model - Theory and simulation

Tran, Hai
11069436

October 20, 2013

Abstract
This paper reports a studies into ferromagnetic (or likelihood pref-
erential) Potts model, its statistical mechanics meaning, and propose a
simulation tool to simulate the models with given sizes and number of
states (q) using Monte-Carlo method. Simulation performance is also
discussed.

1 Introduction
In statistical mechanics, a Potts model is a model of spins (a.k.a cells or sites)
interacting with each other on a lattice of d dimension. By having more than
two states of spin, which is the number of states a spin could have in Ising model,
Potts models are, in essence, the generalisation of the Ising model. In fact a
Potts model with 2 states is an Ising model with haft the critical temperature.
Ising model was invented by Wilhelm Lenz who gave it to his student
Ernst Ising as the subject of his doctoral dissertation. It originally was a one-
dimensional model to study the theory of ferromagnetism[1]. Today with the
assistance of computers, multi-dimensional Ising models could be constructed
to study models of individuals with binary states. Some applications of Ising
model could be named such as magnetism, lattice gas, and neuroscience.
Potts model on the other hands, is a general model to study order-disorder
transformations[2] without targeting any applications and therefore can be ap-
plied in many models. Like Ising model, Potts models are used to find the
critical points of a system (i.e. the temperature at which ferromagnetism be-
comes paramagnetism, or the velocity of the foam at which it starts flowing
uncontrollably), which makes two-dimensional Potts models more appealing as
the critical points can be found mathematically. With three or more dimen-
sions the critical points can be located only by numerical means [3], therefore
this project only studies two dimensional models.
This project delves into the simulations of the Potts model using Metropolis-
Hastings method. The simulation tool features an interactive interface with
statistical data diagrammed in real time as well as an option to store these data
for later examination.

1
2 Literature Reviews
One of the first contribution to the study of ferromagnetism is the one dimen-
sional Ising model [1], in which Ising solved the problem mathematically and
concluded that there is no phase transition in one dimensional lattices. He
further theorised the non-existence of phase transitions in high-dimensional lat-
tices, which then proven wrong in theory by Peierls[4] and in mathematical form
by Griffiths[5] and Dobrushin[6]. Ising model of three or more dimensions is yet
(and very unlikely) to be solved mathematically, thus later most of the studies
into the model are to approximate it using various algorithms.
The Potts models origins date back to 1943 with a model suggested by
Ashkin and Teller [7], whom later the vector Potts model with four states is re-
ferred to. Their study focuses on the Curie temperature at which ferromagnetic
materials lose their magnetic properties. Intrigued by this study, Cyril Domb
suggested to his PhD student Renfrey Potts a generalisation of Ising model
with q spins each pointing to one of the equally spaced direction specified by
the angles
2n
n = , n = 0, 1, ..., q 1. (1)
q
In this model the interaction between two neighbours depends only on the
angle between two vectors, hence the Hamiltonian:
X
H = J(ij ) (2)
<ij>

where ij is the angle between the two neighbouring vectors (spins) i and j.
Domb suggestion to the function J(ij ) is:

J(ij ) = 1 cos (3)

Potts instead used Kramers-Wanniers analysis[8] and solved the suggested model
with q = 2,3,4, and reported the critical point for all q following the model:

J(ij ) = 2 Kr (ni , nj ) (4)

The model Potts solved using Kramers-Wanniers analysis is called the stan-
dard Potts model, in contrast with the original model suggested by Domb,
which is called the planar Potts model. Two models are not related, and only
the original model is studied in this project.

3 Methodology
3.1 Energy of connexions
Physicists assume the magnetic properties are created by two or more particles
spinning to the same direction and contribute to the total magnetic field. Thus
particles of a ferromagnetic material prefer having neighbours with like spin
and particles with anti-magnetic properties prefer having neighbours with unlike
spins. To put the particles preference into mathematical equations, the energy
between two neighbouring particles are the Hamiltonian given by:

2
 hi,j = Ji,j (5)
In which is the inverse temperature given by:
1
= (6)
kT
where the Boltzmann constant k = 1.38 1023 joules/Kevin and T is the
temperature of the system. Also in Equation (5), the Kronecker delta-function
i,j is given by: 
1 if i = j
i,j = (7)
0 if i 6= j
With q = 2, equation 7 yields haft the energy of the equivalent Ising model,
which is:

i,j = i j (8)
or 
1 if i=j
i,j = (9)
1 if i 6= j
since a spin in Ising model could only be either 1 or -1.

Figure 1: Possible states of a square lattice with two spins. The lowest energy
states are confined in dotted boundary, the highest are in solid.

To reflect the likeness preference, the factor J, or energy of connexion, is to

have negative value if the particles are ferromagnetic, and positive otherwise.
Since this project studies ferromagnetic Potts models, J is to be -1. However, to
simplify the calculation while keeping a comprehensible concept of low energy
preference, in this project the equation 7 is rewritten as:

0 if i = j
i,j = (10)
1 if i 6= j

3
 thus without introducing a factor J = 1, equation 10 still reflects low
energy preference, or ferromagnetism. Following this approach, the lowest and
highest energy states as depicted in Figure 1 are now H = 0 and H = 4
instead of H = 4J and H = 0J, and the medium energy states how have
energy H = 2 instead of H = 2J.

3.2 Probability function

Statistical mechanics, as well as quantum physics utilise the concept of likelihood
or probability theory, which describe the state of a system not by exact quantity
but probability of occurrence. Thus a connexion at a preferable low energy could
still transform to a higher energy state if the probability for such transformation
exists. These probabilities are distributed on a domain with the distribution
function:

eh(s)
P (s) = P h(s) (11)
se
in which the numerator is simply an exponential value of the Hamiltonian of a
particular state given in Equation (5) and is easy to calculate. The denominator,
usually called the partition function, is sum of all possible numerators, and is
hard to calculate for any system with meaningful lattice size.
For instance, in a system depicted in Figure 1, the probability for having
four white particles (more accurately four particles in white state) is:
e4J
p(sw,w,w,w ) = 12e2J +2e 4J +2e0J

In simulation the constant factor Jk is usually regards as -1 (ferromagnetic) or

1 (anti-ferromagnetic) as J and k cancel each other, thus the above probability
can be abstractedly calculated as:
4
e T
p(sw,w,w,w ) = 2 4 (12)
12 e T + 2 e T + 2
so the probability of occurrence of a state is now depend solely on the tem-
perature of the system. Table 1 shows how the probability changes regarding
the systems temperature.

Temperature Probability
0.01 0.5 or 50%
1.00 0.27 or 27%
2.00 0.149 or 14.9%
10.00 0.076 or 7.6%
1
100.00 0.0637 or approx. 16

Table 1: Probability for the system in Figure 1 to have all four particles in
white state.

As can be seen the probability of having all particles with the same state
is highest at low temperature, when the ferromagnetic particles prefer same
direction alignment, and the free energy in the system is not enough for most
of the particles to change state. Higher temperature brings more free energy to
the system, and when the temperature is high enough the probability for the

4
particles to stay aligned is approximately the same as when they stay misaligned.
It can also be observed that under any temperature, alignment is still preferred
over misalignment, and thus the probability for the system in Figure 1 to have
1
all particles in one state (black or white) is never equal or less than 16 .
Calculating the partition function is hard especially for system of meaningful
lattice size. For instance with d = 2 and q = 3, a system with 100100 spins
3
will have (100 100)2 = 1065536 states! Calculating the partition function with
this many states is unthinkable. Instead, Potts models consider the relative
probability that the system is in two states and decide whether or not to change
state.

p(s2 ) eh(s2 ) eh(s1 ) eh(s2 ) P

= P h(s) / P h(s) = h(s ) = eJ (s2 s1 ) (13)
p(s1 ) se se e 1

If this relative probability is negative, meaning that the state s2 has higher
probability than s1 , the system changes to s2 , otherwise it changes by the rel-
ative probability. In simulation terms that task is usually done by coining a
normally-distributed number between 0 and 1 and change state if the relative
probability is more than this random number.

3.3 Simulation method

This project features an application written in Java with rich graphical interface
to give users the highest flexibility to control the system. From the control
panel users can pause/resume the system, adjust the temperature, set the ideal
frame rate (FPS) on which the system will be running, change colours of the
states, as well as an option to record simulation data (including temperature,
average energy and magnetisation). Users can also define a new system with
different lattice size and number of states, and start a new system with a hot
start (all spins are random) or cold start (fill the system with the value of a
chosen state). Temperature of the system could be set static (unchanged) or
auto-increment/decrement towards the maximum/minimum temperature (0.0
and 10.0 respectively) after an adjustable number of sweep. For instance, if
chosen to be auto-increment when the temperature is 1.0 and number of sweeps
is 1000, the simulation will increase the temperature by (10.01.0)/1000 = 0.09
after each sweep, and at temperature = 10.0 the simulation will stop.
Each sweep is one iteration through the whole system. As iterating in any
specific order is not accurate, the simulation coin a pair of [x, y] coordinate in
which x [0..lattice width] and y [0..lattice height], and use this coordinate
to pick the corresponding spin at [x, y]. Once picked the relative probability
of the current state and a new, randomly chosen state, is calculated and the
system change according to the method described in Section 3.2. Ideally after
one sweep, all spins will be picked, in reality some will be picked more than
one, while some will not be picked at all. This is the typical application of
Metropolis-Hastings method, named after Metropolis[9] and Hastings[10] for
their invention and contribution to the model.
After the initialisation of a Potts model, total energy and magnetisation is
calculated as:
1X
E= hs (14)
2 s

5
 X
M= s (15)
s

Where the faction 12 in Equation 14 is to justify the total energy after each
connexion between two neighbouring spins is calculated twice.
These calculation is expensive if done once in each sweep. Instead, in this
simulation, they are adjusted each time a spin changes state:
+
E = s2 s1 (16)
+
M = s2 s1 (17)
These total values are divided by the lattice size to get the average energy and
average magnetisation, two main statistical data displayed on the applications
real time charts and also recorded to CSV format if requested. This CSV file is
ready to be used in any spreadsheet program for statistical purpose. Section 4
uses the data generated this way to discuss more about Potts models and their
results.

4 Simulation Result
4.1 Potts model results
The main focus of Potts models is to find the critical point and observe the phe-
nomenon occur during the phase transition between order/disorder states. With
d = 2 the critical point (usually denoted as Tc with T stands for temperature
by historical meaning) can be found mathematically[3] as follows:
1
Tc = (18)
ln(1 + q)
in which is the absolute difference between maximum and minimum state
values. In most Potts models this value is = 1 0 = 1, while in Ising model
it is = 1 (1) = 2. q is the number of states in the system.

q Tc
2 1.135
3 0.995
4 0.91
5 0.852
... ...
10 0.70
... ...

Table 2: Critical points in systems with different number of states.

To verify the model, data from different runs of Potts models with lattice
size [500 500] and q = 2, 4, 5, 10 is recorded and analysed, the results shown
in Figure 2 and 3 indeed verify the phase transitions occurrence at around
the critical temperature given in Tabel 2. It also shows the sharpness of the
transitions to be sharper when q increases. In fact with Q 4 the transition is

6
Figure 2: Average energy of various Potts models when temperatures auto
increment from 0.0 to 2.0.

Figure 3: Average energy of various Potts models when temperatures auto

increment from 0.0 to 2.0.

7
second order, or continuous, and with Q > 4 the model exhibits a first order,
abrupt transition from order to disorder.
Potts models with q > 2 dont have physical magnetisation and the word
itself carries a different implication depending on the particular application.
For the sole purpose of presenting, the data was recorded from cold-start Potts
models with value of all spin initialised to 0 at start. Unlike average energy which
could only vary from 0 to 2, magnetisation varies from 0 to q 1 and therefore
the data presented in Figure 3 has been normalised to the range [0.0,1.0]. As can
be observed the average magnetisation gradually increase toward the average
value (0.5), and the models transform from magnetised to non-magnetised at the
critical temperatures. This phenomenon occurs concurrently with the increase
of energy from low to high.
Alternatively if a cold-start with initial value different from 0, or a hot-
start with random initial value are used, this average magnetisation will be
misleading, as a completely ordered system can also have average magnetisation
at 0.5, simply by having all spins stabilised at the average value state. (E.g.
in a q=3 Potts model, if all spins have state value 1, the system is order and
magnetised but the average magnetisation (0.5) suggests the opposite).

4.2 Performance
On an average computer, the simulation could handle Potts model with lattice
size [300300] at the frame rate of 60fps. This number decreases sharply when
the lattice size increases beyond [300300]. At the maximum size of [500500]
the simulation runs at 1618fps.
This frame rate cannot be achieved without applying several techniques to
improve the performance. A naive version of this simulation can only produce
half the frame rate, here naive refers to a non-optimised version. Thus the
optimisations applied in this simulation is an important topic worth discussing
here, they are:

Using a faster (and more accurate) Random Number Generator (RNG)

algorithm.

Using a faster (though subtly less accurate) exponential function.

Using arrays of one dimension instead of two.

The faster RNG mentioned here is a XOR-Shift RNG proposed by George

Marsaglia[11], which is about 20 times faster than Javas built-in RNG. It is
also more accurate, as found by LEcuyer and Simard[12] to fail only 7 out of
160 tests of their BigCrush battery test bench. Javas RNG on the other hand,
is a mere Linear Congruential Generator failing 21 out of 160 tests1 . Not only
improving the overall accuracy of the system, using XOR-Shift RNG also helps
accelerating the simulation by 2025%.
Since there is no absolute accurate value of Eulers number e, all exponential
functions are approximation of the definition of ex :
1 The Javas RNG mentioned here is java.util.Random, since Java has other RNGs with

higher randomness but are much slower, i.e. java.security.SecureRandom is cryptographically

strong[13] but is 30 times slower than java.util.Random

8
 static int x = 123456789;
static int y = 362436069;
static int z = 521288629;
static int w = 88675123;
int t;
public int nextInt() {
t = x ^ (x << 11);
x = y; y = z; z = w;
return w = w ^ (w >> 19) ^ (t ^ (t >> 8));
}

Figure 4: XOR-Shift RNG algorithm in Java

x n
ex = lim (1 + ) (19)
n n
Obviously the accuracy and speed of an exponential function depends solely
on the number n, but n needs to be considerably large to produce statistically
accurate results. One solution is to build a look up table with highly accurate
predefined exponential value in a certain range. Nicol Schraudolph[14] proposed
an alternatively fast exponential function that utilises the exponent part of float-
ing point number representation in memory to produce much faster calculation
with a subtle inaccuracy. This simulation follows his approach in calculating
exponential value.

public static double exp(double val) {

final long tmp = (long) (1512775 * val + 1072632447);
return Double.longBitsToDouble(tmp << 32);
}

Figure 5: Schraudolphs algorithm in Java

( )
In this simulation, values passed to the exponential function ( s2T s1 )
are mostly in range (-20, 0). In this range the error of Schraudolphs algorithm
compared to Javas Math.exp() result is averagely 1.5% and 3.1% in the
worst cases. Highly accurate applications might not accept this error rate, in
which case this simulation could be modified to fall back to Javas Math.exp()
method. Nevertheless, Schraudolphs algorithm in the range (-20,0) is amazingly
300500 times faster and boost the overall performance by as much as 50%.
Javas JVM treats multi-dimensional arrays internally instead of storing
them as linearly with implicit conversion like the way C does. Thus it seems
trivial for this simulation to store arrays in one dimension and to explicitly
convert one dimensional indices to two dimensional coordinates. However, as
Metropolis algorithm requires random read and write to the arrays, using one
dimensional arrays reduce the workload by generating only one random num-
ber for the index instead of one for each dimension. Furthermore, all arrays
in this simulation has the same size, and thus their index perfectly match and
requires no one-to-two dimensional conversion. Experiment shows that using
one dimensional arrays further increases the performance by 10%.

9
 In total, when all of the above performance tweaks are applied, the simula-
tion enjoys a more than 100% increase in performance.

5 Applications of the Pott model

Recently Potts and Ising models have become subjects to a vast fields of appli-
cations. Google Scholar lists more than 247,000 publications on Ising model and
177,000 on Potts model up to the date of this report. Potts model in particular
can be found in many scientific models, including Computer Science (i.e. Neural
Networks), Physics (i.e. liquid-gas transitions, foam behaviours), Biology (i.e.
Protein Folds, Biological Membranes), Medication (i.e. cancer cells develop-
ment, beating heart cells), or Sociology (i.e. Flocking birds, Social Behaviours,)
etc.
Applications in these fields share the same model but might be mathemati-
cally and implicatively different from the original Potts model. For instance, a
study into foam drainage[15] using the Hamiltonian:
X X
H= J(1 i j ) + (an An )2 (20)
i,j n

with the area of foam bubbles accounted in the total energy, found that
at the critical velocity, the foam starts to flow uncontrollably, and that larger
bubbles flow faster. Another study into the development of tumorous cells[16]
modifies the energy function to account for the volume of a cell and the amount
of nutrients, hence the following Hamiltonian:

XX X
H= J (ij ) (i0 j0 ) {1 ij i0 j0 } + (v VT )2 + Kp(i, j) (21)
ij ij

They found that tumours grow exponentially in the beginning, until the
critical level of nutrients is reached, when they start to grow slower and exhibit
a trend of migration toward nutrients.
The Potts model also contributes to Social Science, typically in Sociology
to study the behaviours of each individuals based on sociological variables such
as the preferences of individuals, sizes of the neighbourhoods or the number of
individuals in a society. Lately, one of such model has helped T.C.Schelling
won the 2005 Nobel prize on economics for his work: Dynamic models of
segregation[17].

6 Conclusion
This paper explains the formation of the Potts model in mathematical terms,
how it is simulated on computer and its applications in real life. This paper
also comes with a simulation tool with great flexibility to help visualising and
analysing the ferromagnetic Potts model in particular and the likelihood pref-
erential Potts model in general. The tools is optimised and can handle large
Potts models, thus its optimisation is also discussed in details.
Potts model is a great tools to study many probabilistic models of different
disciplines. It is very likely that Potts model is gaining more and more research

10
Figure 6: Number of publications on Ising model 1966-1998

11
interests in the future, as it has proven to be a successful research tool in recent
days.

References
[1] E. Ising. Beitrag zur Theorie des Ferromagnetismus. Zeitschrift fur Physik,
31:253258, February 1925.

[2] Renfrey B. Potts. Some generalized Order-Disorder transformation. In

Transformations, Proceedings of the Cambridge Philosophical Society, vol-
ume 48, pages 106109, 1952.
[3] F. Y. Wu. The potts model. Reviews of Modern Physics, 54(1):235268,
January 1982.

[4] R. Peierls. On isings model of ferromagnetism. Mathematical Proceedings

of the Cambridge Philosophical Society, 32:477481, 10 1936.
[5] Robert B. Griffiths. Peierls proof of spontaneous magnetization in a Two-
Dimensional ising ferromagnet. Physical Review Online Archive (Prola),
136(2A):A437A439, October 1964.

[6] L. Dobrushin, R. Existence of a phase transition in two-dimensional and

three-dimensional ising models. Theory Probab. Appl., 10:193213, 1965.
[7] J. Ashkin and E. Teller. Statistics of two-dimensional lattices with four
components. Phys. Rev., 64:178184, Sep 1943.

[8] H. A. Kramers and G. H. Wannier. Statistics of the Two-Dimensional

Ferromagnet. Part I. Physical Review, 60:252262, August 1941.
[9] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Au-
gusta H. Teller, and Edward Teller. Equation of state calculations by fast
computing machines. The Journal of Chemical Physics, 21(6):10871092,
1953.
[10] W. K. Hastings. Monte carlo sampling methods using markov chains and
their applications. Biometrika, 57(1):97109, April 1970.
[11] George Marsaglia. Xorshift RNGs. Journal of Statistical Software, 8(14):1
6, July 2003.

[12] Pierre LEcuyer and Richard Simard. TestU01: A c library for empirical
testing of random number generators. ACM Trans. Math. Softw., 33(4),
August 2007.
[13] Oracle Inc. Class SecureRandom.

[14] Nicol N. Schraudolph. A fast, compact approximation of the exponential

function. Neural Computation, 11(4):853862, 1999.
[15] Yi Jiang and James A. Glazier. Foam drainage: Extended large-q potts
model simulation, 1996.

12
[16] L. Sun, Y. F. Chang, and X. Cai. a Discrete Simulation of Tumor Growth
Concerning Nutrient Influence. International Journal of Modern Physics
B, 18:26512657, 2004.

[17] Thomas C. Schelling. Dynamic models of segregation. The Journal of

Mathematical Sociology, 1(2):143186, July 1971.

13