entropic.tex
Competitive nucleation and the Ostwald rule
in a generalized Potts model with multiple metastable phases
David P. Sanders,∗ Hernán Larralde, and François Leyvraz
Instituto de Ciencias Fı́sicas, Universidad Nacional Autónoma de México, Apartado postal 48-3, 62551 Cuernavaca, Morelos, Mexico
(Dated: November 4, 2018)
We introduce a simple nearest-neighbor spin model with multiple metastable phases, the number and decay
pathways of which are explicitly controlled by the parameters of the system. With this model we can construct,
for example, a system which evolves through an arbitrarily long succession of metastable phases. We also
construct systems in which different phases may nucleate competitively from a single initial phase. For such a
system, we present a general method to extract from numerical simulations the individual nucleation rates of the
nucleating phases. The results show that the Ostwald rule, which predicts which phase will nucleate, must be
modified probabilistically when the new phases are almost equally stable. Finally, we show that the nucleation
rate of a phase depends, among other things, on the number of other phases accessible from it.
PACS numbers: 64.60.My, 64.60.Qb, 05.10.Ln, 05.50.+q
Metastability is a ubiquitous phenomenon in nature.
Broadly speaking, it occurs when a system is “trapped” in
a phase different from equilibrium. This non-equilibrium
phase, the metastable state, can last for extremely long times.
Thus, it is not surprising that metastable states play a cru-
cial role in many physical processes and are at the center of
much current research. For example, recently an intermedi-
ate metastable phase was shown to provide an easier pathway
for the growth of crystal nuclei from fluids (nucleation), with
implications for the crystallization of proteins [1, 2]. Proteins
themselves are known to get stuck in misfolded metastable
structures [3], preventing them from reaching their equilib-
rium configuration. The phenomenology observed in these
and in many other systems can be thought of as arising from a
complicated energy landscape, with several local “metastable
minima” where the trapping occurs [4]. The extreme situation
is that of glasses, in which the energy landscape can have ex-
temely many local minima hindering relaxation of the system
to a thermodynamically stable crystal [5].
The above systems present at least several metastable states.
These states and the transitions between them usually arise
from the microscopic interactions in a complicated way.
When this is the case, the study of phenomena such as compe-
tition between nucleating phases and specific nucleation path-
ways may be obscured. In view of this, in this work we
present a simple spin model with nearest-neighbour interac-
tions, where the number of metastable phases and the decay
pathways between them can be explicitly specified by varying
the model parameters. It thus serves as a test-bed for theoreti-
cal results relating to systems with multiple metastable phases
[6, 7, 8], just as the kinetic Ising model, a special case of our
model, has been central in the study of systems with a sin-
gle metastable phase [9]. As discussed below, the model also
describes the adsorption of multiple chemical species onto a
surface, an interesting physical problem in its own right.
After presenting the model, as an illustration of a possible
application, we construct a system with arbitrarily long suc-
cessions of metastable states. We then focus on competition
between phases nucleating from a single initial metastable
phase. An important question in this context is to understand
which phases nucleate under which conditions. The Ostwald
rule states that the nucleating phase is the one with the small-
est free energy barrier from the initial phase: see Ref. [10]
and references therein. Previous results have supported this
prediction [11].
We show that in general the Ostwald rule must be modified
probabilistically when the new phases are of similar stabil-
ity, using an argument based on individual nucleation rates of
each phase. We give a method by which these rates can be
measured in simulations or experiments, and show that there
is a parameter regime in which any of the new phases may
nucleate—only the nucleation probability of each phase can
be established, with the outcome in any given run being un-
predictable. We finally show that the nucleation probability of
a phase depends on the phases accessible from it.
Model details:- Our model is based on the Potts model, in
which each spin has one of q states [12] and each phase has a
majority of spins in one state; the Ising model corresponds to
q = 2. The relative stability of each phase is controlled by ex-
ternal fields, and the interplay of these fields with interactions
between different spin states allows us to obtain any desired
transition pathways between phases.
Viewing the fields as chemical potentials, we can recast the
model as a multi-component lattice gas which describes ad-
sorption on a lattice substrate (e.g., a crystal plane) of mul-
tiple chemical species with lateral interactions [13]. Much
experimental work has been done on the thermodynamics of
such systems, but little on the kinetics—see [14] and refer-
ences therein; nonetheless, our results should be testable in
that context. A more complicated system where the kinetics
has been characterized is a colloid–polymer system [15, 16],
where possible pathways were found from considerations of
the free energy landscape [17]. Our approach is complemen-
tary in that specific pathways result from microscopic interac-
tions.
We work on an L×L square lattice with N := L2 spins and
periodic boundary conditions, although the results are qualiti-
tavely unaffected by lattice type. Each lattice site i has a spin
http://arxiv.org/abs/0704.0472v3
σi taking values in {1, . . . ,q}, and the energy of a configura-
tion ß is given by the Hamiltonian
H(ß) :=− ∑
〈i, j〉
Jσi,σ j −
hαMα . (1)
Here, Mα := ∑i δσi,α is the magnetisation (=number of spins)
of the spin type α; δα ,γ = 1 if α = γ , and 0 otherwise. The
first term is a sum over nearest-neighbor pairs of spins of a
symmetric interaction energy Jα ,γ = Jγ,α , and the second de-
scribes the effect of external fields hα acting on spin type α .
We set the diagonal elements Jα ,α of the interaction ma-
trix to unity (Jα ,α := 1 for all α), so that in the absence of
non-diagonal interactions and fields, the model reduces to the
standard Potts model [12]. This has q symmetrical phases
coexisting below a critical temperature Tc = 1/(ln(1+
each phase has a majority of spins in one of the q spin states.
Including fields breaks the symmetry between phases. If
hα = hγ , then the α and γ phases coexist, with a first-order
phase transition between them, for T < Tc; this is at the ori-
gin of metastability. Weak non-diagonal interactions do not
qualititavely affect this coexistence.
To evolve the system we choose discrete-time Metropolis
dynamics [18]: at each time step, a spin and its new value
are chosen at random, the increment ∆H of the Hamiltonian
(1) for this change is calculated, and the update is accepted
with probability min{1,exp(−β ∆H)}, where β := 1/T is the
inverse temperature. This gives a Markov chain on the space
of all possible configurations.
This Markov chain has a unique equilibrium distribution,
concentrated on the phase(s) with the largest hα . The other
phases are metastable, that is, when started in such a phase α ,
the system stays there for some time, before a transition to a
more stable phase γ is nucleated by the appearance of a crit-
ical droplet of the γ phase. At sufficiently low temperatures,
the relative stability is determined by hγ > hα . The reverse
transition is exponentially unlikely.
In the standard Potts model, the equilibrium phase (almost)
always nucleates. To obtain non-trivial transition pathways,
nucleation of other phases must be promoted. This we achieve
using non-diagonal interactions between distinct spin types
α 6= γ: setting Jα ,γ > 0 favors nucleation of γ droplets inside
the α phase by lowering the surface tension between α and γ
regions, and hence decreasing the droplet free energy of for-
mation (nucleation barrier), whereas formation of γ droplets
in the α phase is suppressed if Jα ,γ < 0.
We can now construct models whose phases obey arbitrary
metastable transition graphs. These are directed graphs with
the restriction that no loops returning to a previously visited
phase are allowed. Each vertex corresponds to one phase, la-
beled by its dominant spin state, and each arrow to a desired
transition: α → γ means that phase γ can nucleate directly
from phase α . Fig. 1 shows example transition graphs.
To construct a model corresponding to a given transition
graph, we proceed as follows. The number of spin types, q, is
the number of vertices in the graph. To each spin type α we
2 3 4
FIG. 1: Metastable transition graphs: (a) kinetic Ising model; (b) suc-
cession of 3 phases; (c) single metastable phase decaying to two com-
peting phases; (d) as in (c), but such that one phase can decay further;
(e) three competing phases.
0 10000 20000 30000
time (Monte Carlo steps per site)
1 2 3 4 5
FIG. 2: (Color online) Time dependence of magnetisation per site
Mα/N of each phase α , in a single run of the model with transition
graph 1 → 2 → 3 → 4 → 5 (a succession of phases). Parameters are
L = 50, β = 1.25, hα = 0.1(α −1), K1 = 0.1, and K2 = 1.0. Labels
denote the dominant phase. Configuration snapshots depict post-
critical nuclei of each new phase embedded in the previous phase.
These grow to fill the system, producing the next phase in sequence.
assign a field hα , with hγ > hα if γ is below α in the graph.
The off-diagonal interactions are given by Jα ,γ := K1 > 0 (at-
tractive) if α → γ , and Jα ,γ :=−K2 < 0 (repulsive) otherwise.
K2 must be large enough to inhibit immediate formation of
non-adjacent phases with large fields.
As an illustration, we construct a model exhibiting a lin-
ear succession of metastable phases with transition graph
1 → 2 → ··· → q. We impose fields 0 = h1 < h2 < · · · < hq
and attractive interactions Jα ,α±1 := K1 > 0 between neigh-
bouring states, and set all other non-diagonal interactions to
−K2 < 0. With suitable, moderately robust, parameters, we
observe the desired behavior, shown for q = 5 in Fig. 2.
A three-phase succession was previously observed in a ki-
netic Blume–Capel model [19, 20], corresponding to a special
case of our model with q = 3 [21]. The physical reason for the
observed transitions is, however, much more transparent with
the Hamiltonian in the form (1), with its intuitive interpreta-
tion in terms of attractive and repulsive interactions.
Competitive nucleation:- We now turn to the decay of one
metastable phase into two competing phases (Fig. 1(c)). Sear
studied competitive heterogeneous nucleation (occurring on
impurities) in the 3-state standard Potts model [11]. In con-
trast, all behaviors discussed in this work are endogenous: ob-
0 0.1 0.2 0.3 0.4 0.5
(h3 −h2)/h2
0 50 100 150
β = 1.25
β = 1.40
L = 150
L = 200
h3 = 0.105
h3 = 0.11
FIG. 3: (Color online) Probability p2(h3) that phase 2 nucleates be-
fore phase 3 as a function of ∆h/h2, with h2 = 0.1, K1 = 0.1 and
K2 = 1, and β = 1.25, L = 50 unless otherwise noted. Up to n = 104
trials were used for each data point; statistical errors are of the or-
der of the symbol size. Inset: system-size dependence of p2 for two
values of h3 with fixed parameter values.
served transitions are not caused by external influences, but
rather arise spontaneously from within the system itself.
We fix 0 = h1 < h2 ≤ h3, J1,2 = J1,3 > 0 and J2,3 < 0. Let
∆h := h3 − h2 be the field difference between the new phases,
2 and 3. When ∆h = 0, these phases are symmetrical, each
nucleating half of the time, while for ∆h > 0, we expect the
1–3 free energy barrier to be lower than the 1–2 one, so that
according to the Ostwald rule, only phase 3 should nucleate.
To test this, we perform n simulations starting from phase 1
for each ∆h, in n2 of which phase 2 nucleates before phase 3.
The ratio n2/n is then an estimate of the probability p2(∆h)
that phase 2 nucleates first. For efficiency, we use a rejection-
free version of the Metropolis method [21, 22].
Fig. 3 plots p2 as a function of a non-dimensionalised ∆h.
For ∆h sufficiently close to 0, phase 2 can still nucleate first,
contrary to the simple Ostwald rule. The probability that it
does so rapidly decreases for larger ∆h, until a point beyond
which phase 3 effectively always nucleates.
To explain these results in a general context, we assume, as
in classical nucleation theory [23], that there are well-defined
nucleation rates λi(L) of phases i = 2,3, giving the number of
critical nuclei which form per unit time in a system of size L.
The nucleation rates per site are µi(L) := λi(L)/N.
A nucleation rate is the inverse of a mean nucleation time,
which can be measured in experiments or simulations by aver-
aging over many nucleation events in independent runs. In the
case of competitive nucleation, however, we can only measure
the mean time τ for the first phase to nucleate, after which
this phase invades the entire system. The rate of this first
nucleation is λ2 + λ3, since the total number of nucleation
events per unit time is the sum of those for each type, so that
τ = 1/(λ2 + λ3). For convenience, in simulations τ is taken
to be the time until the new phase occupies half the system.
Under the same assumptions, the probability that phase 2
0 0.1 0.2 0.3 0.4
(h3 −h2)/h2
λ2 direct
λ2 forward-flux
λ3 direct
λ3 forward-flux
-0.2 0 0.2 0.4 0.6
(h2 −h3)/h3
with phase 4
without phase 4
FIG. 4: (Color online) (a) Comparison of nucleation rates λi(h3) of
phases i = 2,3 from phase 1 in the model of Fig. 1(c) with h2 fixed,
calculated directly from simulations and using forward-flux sam-
pling. (b) Nucleation probability p2(h2) for h3 = 0.1 and h4 = 0.2
varying h2, with and without phase 4, in the model of Fig. 1(d).
Dashed lines show positions of equal nucleation probability and
equal field of the two phases. In both subfigures parameters are as in
Fig. 3, with L = 50 and β = 1.25.
nucleates first is p2 = Prob{T2 ≤ T3}, where Ti, the time for
phase i to nucleate, is an exponentially distributed random
variable with mean 1/λi. This gives p2 = λ2/(λ2+λ3) = λ2τ .
Individual nucleation rates of the two phases, which can-
not be obtained directly, can now be calculated as: λ2 =
p2/τ and λ3 = (1− p2)/τ . This generalizes to P compet-
ing phases, where the measurable quantities are the mean nu-
cleation time τ = 1/∑Pi=1 λi and the nucleation probabilities
pi = λi/∑ j λ j = µi/∑ j µ j of each phase i. The nucleation
rate of phase i is then λi = pi/τ .
Figure 4(a) shows λ2 and λ3 calculated in this way. To con-
firm the validity of such calculations, we use the “forward-flux
sampling” method [24, 25], which directly calculates the tran-
sition rate between two phases in a stochastic system. This
has previously been used to study nucleation rates in the Ising
model [26, 27]. In our case, the possibility of escape to an
additional new phase must be taken into account [21]. Fig-
ure 4(a) shows that the results indeed coincide with those of
the direct method, within statistical errors. We remark that we
are unaware of any analytical prediction giving the observed
variation of nucleation rates.
The above considerations used “in reverse” confirm that
the Ostwald rule must in general be modified when the new
phases have similar stabilities, as follows. Consider phases
which are equally stable for given parameter values. We ex-
pect nucleation barriers, and hence nucleation rates, to vary
continuously with the parameters, so that the nucleation prob-
abilities pi also vary continuously. Hence there is a region,
where the phases have similar stabilities, in which all nucle-
ation probabilities are non-zero—only the probability of each
phase nucleating is well-defined, with the outcome in any
given run being stochastic, as in Fig. 3. The definite prediction
given by the Ostwald rule is thus invalid in this region.
To see how our results depend on system size L, we note
FIG. 5: (Color online) Configuration snapshots of the multidroplet
regime in the model with transition graph Fig. 1(e), with L = 200,
β = 1.25, and h1 = 0, h2 = h3 = h4 = 0.15. Nucleation of droplets
of three equally stable phases from a single metastable phase is fol-
lowed by droplet growth and then domain coarsening. The “coating”
of the initial phase visible between domains in the final snapshot is
due to repulsive interactions between the new phases.
that in a broad region of L, the per-site nucleation rates µi,
and hence also the pi, are independent of L, as confirmed by
the plateaus in the inset of Fig. 3. This is valid when the nu-
cleation process is mediated by growth of a single droplet [9].
Note that this regime may be of relevance for macroscopically
large systems [9].
Above a certain system size, however, droplets of differ-
ent phases may nucleate before any dominates the system
(the ‘multidroplet’ regime) [9]. This results in coarsening,
as shown in Fig. 5 for three competing phases of equal sta-
bility. Even if a phase-α droplet nucleates first, droplets of
a more-stable phase may then nucleate and grow to dominate
the system before the α phase can do so, thus reducing pα , as
seen in Fig. 3.
Finally, a generic non-symmetric case can be obtained by
adding a new phase, 4, and a decay path, 3→ 4, as in Fig. 1(d).
Even when h2 = h3, phase 3 now nucleates more often, as
shown by the horizontally displaced nucleation probability
curve in Fig. 4(b): the presence of phase 4 reduces p2 by
roughly half. This is due to an entropic effect: there are more
3-dominated critical droplets than 2-dominated ones, since
spins of type 4 also appear in the former droplets, resulting
in nucleation of a binary mixture [23]; a similar effect is vis-
ible in Fig. 1. There is thus a lower free energy barrier to
form phase 3, and yet nucleation of both phases 2 and 3 is still
observed, again at odds with the simple Ostwald rule.
In summary, we have introduced a generalized Potts model
which can easily be tuned to have any given number of
metastable phases and arbitrary transitions between them. We
have shown generally that individual nucleation rates of com-
petitively nucleating phases can be calculated from exper-
imentally measurable quantities, and that the Ostwald rule
must be modified when the nucleating phases have compa-
rable stabilities. In future work [21], we will study the model
in detail and compare its properties with theoretical results [8]
on systems with multiple metastable phases.
DPS thanks A. Huerta for useful discussions and the Uni-
versidad Nacional Autónoma de México for financial sup-
port. The financial support of DGAPA-UNAM project PA-
PIIT IN112307 is also acknowledged.
∗ Electronic address: dsanders@fis.unam.mx
[1] J. F. Lutsko and G. Nicolis, Phys. Rev. Lett. 96, 046102 (2006).
[2] P. R. ten Wolde and D. Frenkel, Science 277, 1975 (1997).
[3] S. Takada and P. G. Wolynes, Phys. Rev. E 55, 4562 (1997).
[4] D. Wales, Energy Landscapes (Cambridge University Press,
Cambridge, 2003).
[5] G. Biroli and J. Kurchan, Phys. Rev. E 64, 016101 (2001).
[6] A. Bovier, M. Eckhoff, V. Gayrard, and M. Klein, Comm. Math.
Phys. 228, 219 (2002).
[7] B. Gaveau and L. S. Schulman, Phys. Rev. E 73, 036124 (2006).
[8] H. Larralde, F. Leyvraz, and D. P. Sanders, J. Stat. Mech. 2006,
P08013 (2006).
[9] P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides, Phys.
Rev. E 49, 5080 (1994).
[10] P. R. ten Wolde and D. Frenkel, Phys. Chem. Chem. Phys. 1,
2191 (1999).
[11] R. P. Sear, J. Phys. Cond. Matt. 17, 3997 (2005).
[12] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
[13] P. Rikvold, J. Collins, G. Hansen, and J. Gunton, Surf. Sci. 203
(1988).
[14] S. Manzi, W. Mas, R. Belardinelli, and V. Pereyra, Langmuir
20, 499 (2004).
[15] W. C. K. Poon, F. Renth, R. M. L. Evans, D. J. Fairhurst, M. E.
Cates, and P. N. Pusey, Phys. Rev. Lett. 83, 1239 (1999).
[16] F. Renth, W. C. K. Poon, and R. M. L. Evans, Phys. Rev. E 64,
031402 (2001).
[17] R. M. L. Evans, W. C. K. Poon, and F. Renth, Phys. Rev. E 64,
031403 (2001).
[18] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in
Statistical Physics (Oxford University Press, New York, 1999).
[19] E. N. M. Cirillo and E. Olivieri, J. Stat. Phys. 83, 473 (1996).
[20] T. Fiig, B. M. Gorman, P. A. Rikvold, and M. A. Novotny, Phys.
Rev. E 50, 1930 (1994).
[21] D. P. Sanders, H. Larralde, and F. Leyvraz, (in preparation).
[22] M. A. Novotny, Phys. Rev. Lett. 74, 1 (1995).
[23] P. G. Debenedetti, Metastable Liquids (Princeton University
Press, Princeton, 1996).
[24] R. J. Allen, P. B. Warren, and P. R. ten Wolde, Phys. Rev. Lett.
94, 018104 (2005).
[25] R. J. Allen, D. Frenkel, and P. R. ten Wolde, J. Chem. Phys.
124, 024102 (2006).
[26] R. P. Sear, J. Phys. Chem. B 110, 4985 (2006).
[27] A. J. Page and R. P. Sear, Phys. Rev. Lett. 97, 065701 (2006).
mailto:dsanders@fis.unam.mx

  We introduce a simple nearest-neighbor spin model with multiple metastable
phases, the number and decay pathways of which are explicitly controlled by the
parameters of the system. With this model we can construct, for example, a
system which evolves through an arbitrarily long succession of metastable
phases. We also construct systems in which different phases may nucleate
competitively from a single initial phase. For such a system, we present a
general method to extract from numerical simulations the individual nucleation
rates of the nucleating phases. The results show that the Ostwald rule, which
predicts which phase will nucleate, must be modified probabilistically when the
new phases are almost equally stable. Finally, we show that the nucleation rate
of a phase depends, among other things, on the number of other phases
accessible from it.

entropic.tex
Competitive nucleation and the Ostwald rule
in a generalized Potts model with multiple metastable phases
David P. Sanders,∗ Hernán Larralde, and François Leyvraz
Instituto de Ciencias Fı́sicas, Universidad Nacional Autónoma de México, Apartado postal 48-3, 62551 Cuernavaca, Morelos, Mexico
(Dated: November 4, 2018)
We introduce a simple nearest-neighbor spin model with multiple metastable phases, the number and decay
pathways of which are explicitly controlled by the parameters of the system. With this model we can construct,
for example, a system which evolves through an arbitrarily long succession of metastable phases. We also
construct systems in which different phases may nucleate competitively from a single initial phase. For such a
system, we present a general method to extract from numerical simulations the individual nucleation rates of the
nucleating phases. The results show that the Ostwald rule, which predicts which phase will nucleate, must be
modified probabilistically when the new phases are almost equally stable. Finally, we show that the nucleation
rate of a phase depends, among other things, on the number of other phases accessible from it.
PACS numbers: 64.60.My, 64.60.Qb, 05.10.Ln, 05.50.+q
Metastability is a ubiquitous phenomenon in nature.
Broadly speaking, it occurs when a system is “trapped” in
a phase different from equilibrium. This non-equilibrium
phase, the metastable state, can last for extremely long times.
Thus, it is not surprising that metastable states play a cru-
cial role in many physical processes and are at the center of
much current research. For example, recently an intermedi-
ate metastable phase was shown to provide an easier pathway
for the growth of crystal nuclei from fluids (nucleation), with
implications for the crystallization of proteins [1, 2]. Proteins
themselves are known to get stuck in misfolded metastable
structures [3], preventing them from reaching their equilib-
rium configuration. The phenomenology observed in these
and in many other systems can be thought of as arising from a
complicated energy landscape, with several local “metastable
minima” where the trapping occurs [4]. The extreme situation
is that of glasses, in which the energy landscape can have ex-
temely many local minima hindering relaxation of the system
to a thermodynamically stable crystal [5].
The above systems present at least several metastable states.
These states and the transitions between them usually arise
from the microscopic interactions in a complicated way.
When this is the case, the study of phenomena such as compe-
tition between nucleating phases and specific nucleation path-
ways may be obscured. In view of this, in this work we
present a simple spin model with nearest-neighbour interac-
tions, where the number of metastable phases and the decay
pathways between them can be explicitly specified by varying
the model parameters. It thus serves as a test-bed for theoreti-
cal results relating to systems with multiple metastable phases
[6, 7, 8], just as the kinetic Ising model, a special case of our
model, has been central in the study of systems with a sin-
gle metastable phase [9]. As discussed below, the model also
describes the adsorption of multiple chemical species onto a
surface, an interesting physical problem in its own right.
After presenting the model, as an illustration of a possible
application, we construct a system with arbitrarily long suc-
cessions of metastable states. We then focus on competition
between phases nucleating from a single initial metastable
phase. An important question in this context is to understand
which phases nucleate under which conditions. The Ostwald
rule states that the nucleating phase is the one with the small-
est free energy barrier from the initial phase: see Ref. [10]
and references therein. Previous results have supported this
prediction [11].
We show that in general the Ostwald rule must be modified
probabilistically when the new phases are of similar stabil-
ity, using an argument based on individual nucleation rates of
each phase. We give a method by which these rates can be
measured in simulations or experiments, and show that there
is a parameter regime in which any of the new phases may
nucleate—only the nucleation probability of each phase can
be established, with the outcome in any given run being un-
predictable. We finally show that the nucleation probability of
a phase depends on the phases accessible from it.
Model details:- Our model is based on the Potts model, in
which each spin has one of q states [12] and each phase has a
majority of spins in one state; the Ising model corresponds to
q = 2. The relative stability of each phase is controlled by ex-
ternal fields, and the interplay of these fields with interactions
between different spin states allows us to obtain any desired
transition pathways between phases.
Viewing the fields as chemical potentials, we can recast the
model as a multi-component lattice gas which describes ad-
sorption on a lattice substrate (e.g., a crystal plane) of mul-
tiple chemical species with lateral interactions [13]. Much
experimental work has been done on the thermodynamics of
such systems, but little on the kinetics—see [14] and refer-
ences therein; nonetheless, our results should be testable in
that context. A more complicated system where the kinetics
has been characterized is a colloid–polymer system [15, 16],
where possible pathways were found from considerations of
the free energy landscape [17]. Our approach is complemen-
tary in that specific pathways result from microscopic interac-
tions.
We work on an L×L square lattice with N := L2 spins and
periodic boundary conditions, although the results are qualiti-
tavely unaffected by lattice type. Each lattice site i has a spin
http://arxiv.org/abs/0704.0472v3
σi taking values in {1, . . . ,q}, and the energy of a configura-
tion ß is given by the Hamiltonian
H(ß) :=− ∑
〈i, j〉
Jσi,σ j −
hαMα . (1)
Here, Mα := ∑i δσi,α is the magnetisation (=number of spins)
of the spin type α; δα ,γ = 1 if α = γ , and 0 otherwise. The
first term is a sum over nearest-neighbor pairs of spins of a
symmetric interaction energy Jα ,γ = Jγ,α , and the second de-
scribes the effect of external fields hα acting on spin type α .
We set the diagonal elements Jα ,α of the interaction ma-
trix to unity (Jα ,α := 1 for all α), so that in the absence of
non-diagonal interactions and fields, the model reduces to the
standard Potts model [12]. This has q symmetrical phases
coexisting below a critical temperature Tc = 1/(ln(1+
each phase has a majority of spins in one of the q spin states.
Including fields breaks the symmetry between phases. If
hα = hγ , then the α and γ phases coexist, with a first-order
phase transition between them, for T < Tc; this is at the ori-
gin of metastability. Weak non-diagonal interactions do not
qualititavely affect this coexistence.
To evolve the system we choose discrete-time Metropolis
dynamics [18]: at each time step, a spin and its new value
are chosen at random, the increment ∆H of the Hamiltonian
(1) for this change is calculated, and the update is accepted
with probability min{1,exp(−β ∆H)}, where β := 1/T is the
inverse temperature. This gives a Markov chain on the space
of all possible configurations.
This Markov chain has a unique equilibrium distribution,
concentrated on the phase(s) with the largest hα . The other
phases are metastable, that is, when started in such a phase α ,
the system stays there for some time, before a transition to a
more stable phase γ is nucleated by the appearance of a crit-
ical droplet of the γ phase. At sufficiently low temperatures,
the relative stability is determined by hγ > hα . The reverse
transition is exponentially unlikely.
In the standard Potts model, the equilibrium phase (almost)
always nucleates. To obtain non-trivial transition pathways,
nucleation of other phases must be promoted. This we achieve
using non-diagonal interactions between distinct spin types
α 6= γ: setting Jα ,γ > 0 favors nucleation of γ droplets inside
the α phase by lowering the surface tension between α and γ
regions, and hence decreasing the droplet free energy of for-
mation (nucleation barrier), whereas formation of γ droplets
in the α phase is suppressed if Jα ,γ < 0.
We can now construct models whose phases obey arbitrary
metastable transition graphs. These are directed graphs with
the restriction that no loops returning to a previously visited
phase are allowed. Each vertex corresponds to one phase, la-
beled by its dominant spin state, and each arrow to a desired
transition: α → γ means that phase γ can nucleate directly
from phase α . Fig. 1 shows example transition graphs.
To construct a model corresponding to a given transition
graph, we proceed as follows. The number of spin types, q, is
the number of vertices in the graph. To each spin type α we
2 3 4
FIG. 1: Metastable transition graphs: (a) kinetic Ising model; (b) suc-
cession of 3 phases; (c) single metastable phase decaying to two com-
peting phases; (d) as in (c), but such that one phase can decay further;
(e) three competing phases.
0 10000 20000 30000
time (Monte Carlo steps per site)
1 2 3 4 5
FIG. 2: (Color online) Time dependence of magnetisation per site
Mα/N of each phase α , in a single run of the model with transition
graph 1 → 2 → 3 → 4 → 5 (a succession of phases). Parameters are
L = 50, β = 1.25, hα = 0.1(α −1), K1 = 0.1, and K2 = 1.0. Labels
denote the dominant phase. Configuration snapshots depict post-
critical nuclei of each new phase embedded in the previous phase.
These grow to fill the system, producing the next phase in sequence.
assign a field hα , with hγ > hα if γ is below α in the graph.
The off-diagonal interactions are given by Jα ,γ := K1 > 0 (at-
tractive) if α → γ , and Jα ,γ :=−K2 < 0 (repulsive) otherwise.
K2 must be large enough to inhibit immediate formation of
non-adjacent phases with large fields.
As an illustration, we construct a model exhibiting a lin-
ear succession of metastable phases with transition graph
1 → 2 → ··· → q. We impose fields 0 = h1 < h2 < · · · < hq
and attractive interactions Jα ,α±1 := K1 > 0 between neigh-
bouring states, and set all other non-diagonal interactions to
−K2 < 0. With suitable, moderately robust, parameters, we
observe the desired behavior, shown for q = 5 in Fig. 2.
A three-phase succession was previously observed in a ki-
netic Blume–Capel model [19, 20], corresponding to a special
case of our model with q = 3 [21]. The physical reason for the
observed transitions is, however, much more transparent with
the Hamiltonian in the form (1), with its intuitive interpreta-
tion in terms of attractive and repulsive interactions.
Competitive nucleation:- We now turn to the decay of one
metastable phase into two competing phases (Fig. 1(c)). Sear
studied competitive heterogeneous nucleation (occurring on
impurities) in the 3-state standard Potts model [11]. In con-
trast, all behaviors discussed in this work are endogenous: ob-
0 0.1 0.2 0.3 0.4 0.5
(h3 −h2)/h2
0 50 100 150
β = 1.25
β = 1.40
L = 150
L = 200
h3 = 0.105
h3 = 0.11
FIG. 3: (Color online) Probability p2(h3) that phase 2 nucleates be-
fore phase 3 as a function of ∆h/h2, with h2 = 0.1, K1 = 0.1 and
K2 = 1, and β = 1.25, L = 50 unless otherwise noted. Up to n = 104
trials were used for each data point; statistical errors are of the or-
der of the symbol size. Inset: system-size dependence of p2 for two
values of h3 with fixed parameter values.
served transitions are not caused by external influences, but
rather arise spontaneously from within the system itself.
We fix 0 = h1 < h2 ≤ h3, J1,2 = J1,3 > 0 and J2,3 < 0. Let
∆h := h3 − h2 be the field difference between the new phases,
2 and 3. When ∆h = 0, these phases are symmetrical, each
nucleating half of the time, while for ∆h > 0, we expect the
1–3 free energy barrier to be lower than the 1–2 one, so that
according to the Ostwald rule, only phase 3 should nucleate.
To test this, we perform n simulations starting from phase 1
for each ∆h, in n2 of which phase 2 nucleates before phase 3.
The ratio n2/n is then an estimate of the probability p2(∆h)
that phase 2 nucleates first. For efficiency, we use a rejection-
free version of the Metropolis method [21, 22].
Fig. 3 plots p2 as a function of a non-dimensionalised ∆h.
For ∆h sufficiently close to 0, phase 2 can still nucleate first,
contrary to the simple Ostwald rule. The probability that it
does so rapidly decreases for larger ∆h, until a point beyond
which phase 3 effectively always nucleates.
To explain these results in a general context, we assume, as
in classical nucleation theory [23], that there are well-defined
nucleation rates λi(L) of phases i = 2,3, giving the number of
critical nuclei which form per unit time in a system of size L.
The nucleation rates per site are µi(L) := λi(L)/N.
A nucleation rate is the inverse of a mean nucleation time,
which can be measured in experiments or simulations by aver-
aging over many nucleation events in independent runs. In the
case of competitive nucleation, however, we can only measure
the mean time τ for the first phase to nucleate, after which
this phase invades the entire system. The rate of this first
nucleation is λ2 + λ3, since the total number of nucleation
events per unit time is the sum of those for each type, so that
τ = 1/(λ2 + λ3). For convenience, in simulations τ is taken
to be the time until the new phase occupies half the system.
Under the same assumptions, the probability that phase 2
0 0.1 0.2 0.3 0.4
(h3 −h2)/h2
λ2 direct
λ2 forward-flux
λ3 direct
λ3 forward-flux
-0.2 0 0.2 0.4 0.6
(h2 −h3)/h3
with phase 4
without phase 4
FIG. 4: (Color online) (a) Comparison of nucleation rates λi(h3) of
phases i = 2,3 from phase 1 in the model of Fig. 1(c) with h2 fixed,
calculated directly from simulations and using forward-flux sam-
pling. (b) Nucleation probability p2(h2) for h3 = 0.1 and h4 = 0.2
varying h2, with and without phase 4, in the model of Fig. 1(d).
Dashed lines show positions of equal nucleation probability and
equal field of the two phases. In both subfigures parameters are as in
Fig. 3, with L = 50 and β = 1.25.
nucleates first is p2 = Prob{T2 ≤ T3}, where Ti, the time for
phase i to nucleate, is an exponentially distributed random
variable with mean 1/λi. This gives p2 = λ2/(λ2+λ3) = λ2τ .
Individual nucleation rates of the two phases, which can-
not be obtained directly, can now be calculated as: λ2 =
p2/τ and λ3 = (1− p2)/τ . This generalizes to P compet-
ing phases, where the measurable quantities are the mean nu-
cleation time τ = 1/∑Pi=1 λi and the nucleation probabilities
pi = λi/∑ j λ j = µi/∑ j µ j of each phase i. The nucleation
rate of phase i is then λi = pi/τ .
Figure 4(a) shows λ2 and λ3 calculated in this way. To con-
firm the validity of such calculations, we use the “forward-flux
sampling” method [24, 25], which directly calculates the tran-
sition rate between two phases in a stochastic system. This
has previously been used to study nucleation rates in the Ising
model [26, 27]. In our case, the possibility of escape to an
additional new phase must be taken into account [21]. Fig-
ure 4(a) shows that the results indeed coincide with those of
the direct method, within statistical errors. We remark that we
are unaware of any analytical prediction giving the observed
variation of nucleation rates.
The above considerations used “in reverse” confirm that
the Ostwald rule must in general be modified when the new
phases have similar stabilities, as follows. Consider phases
which are equally stable for given parameter values. We ex-
pect nucleation barriers, and hence nucleation rates, to vary
continuously with the parameters, so that the nucleation prob-
abilities pi also vary continuously. Hence there is a region,
where the phases have similar stabilities, in which all nucle-
ation probabilities are non-zero—only the probability of each
phase nucleating is well-defined, with the outcome in any
given run being stochastic, as in Fig. 3. The definite prediction
given by the Ostwald rule is thus invalid in this region.
To see how our results depend on system size L, we note
FIG. 5: (Color online) Configuration snapshots of the multidroplet
regime in the model with transition graph Fig. 1(e), with L = 200,
β = 1.25, and h1 = 0, h2 = h3 = h4 = 0.15. Nucleation of droplets
of three equally stable phases from a single metastable phase is fol-
lowed by droplet growth and then domain coarsening. The “coating”
of the initial phase visible between domains in the final snapshot is
due to repulsive interactions between the new phases.
that in a broad region of L, the per-site nucleation rates µi,
and hence also the pi, are independent of L, as confirmed by
the plateaus in the inset of Fig. 3. This is valid when the nu-
cleation process is mediated by growth of a single droplet [9].
Note that this regime may be of relevance for macroscopically
large systems [9].
Above a certain system size, however, droplets of differ-
ent phases may nucleate before any dominates the system
(the ‘multidroplet’ regime) [9]. This results in coarsening,
as shown in Fig. 5 for three competing phases of equal sta-
bility. Even if a phase-α droplet nucleates first, droplets of
a more-stable phase may then nucleate and grow to dominate
the system before the α phase can do so, thus reducing pα , as
seen in Fig. 3.
Finally, a generic non-symmetric case can be obtained by
adding a new phase, 4, and a decay path, 3→ 4, as in Fig. 1(d).
Even when h2 = h3, phase 3 now nucleates more often, as
shown by the horizontally displaced nucleation probability
curve in Fig. 4(b): the presence of phase 4 reduces p2 by
roughly half. This is due to an entropic effect: there are more
3-dominated critical droplets than 2-dominated ones, since
spins of type 4 also appear in the former droplets, resulting
in nucleation of a binary mixture [23]; a similar effect is vis-
ible in Fig. 1. There is thus a lower free energy barrier to
form phase 3, and yet nucleation of both phases 2 and 3 is still
observed, again at odds with the simple Ostwald rule.
In summary, we have introduced a generalized Potts model
which can easily be tuned to have any given number of
metastable phases and arbitrary transitions between them. We
have shown generally that individual nucleation rates of com-
petitively nucleating phases can be calculated from exper-
imentally measurable quantities, and that the Ostwald rule
must be modified when the nucleating phases have compa-
rable stabilities. In future work [21], we will study the model
in detail and compare its properties with theoretical results [8]
on systems with multiple metastable phases.
DPS thanks A. Huerta for useful discussions and the Uni-
versidad Nacional Autónoma de México for financial sup-
port. The financial support of DGAPA-UNAM project PA-
PIIT IN112307 is also acknowledged.
∗ Electronic address: dsanders@fis.unam.mx
[1] J. F. Lutsko and G. Nicolis, Phys. Rev. Lett. 96, 046102 (2006).
[2] P. R. ten Wolde and D. Frenkel, Science 277, 1975 (1997).
[3] S. Takada and P. G. Wolynes, Phys. Rev. E 55, 4562 (1997).
[4] D. Wales, Energy Landscapes (Cambridge University Press,
Cambridge, 2003).
[5] G. Biroli and J. Kurchan, Phys. Rev. E 64, 016101 (2001).
[6] A. Bovier, M. Eckhoff, V. Gayrard, and M. Klein, Comm. Math.
Phys. 228, 219 (2002).
[7] B. Gaveau and L. S. Schulman, Phys. Rev. E 73, 036124 (2006).
[8] H. Larralde, F. Leyvraz, and D. P. Sanders, J. Stat. Mech. 2006,
P08013 (2006).
[9] P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides, Phys.
Rev. E 49, 5080 (1994).
[10] P. R. ten Wolde and D. Frenkel, Phys. Chem. Chem. Phys. 1,
2191 (1999).
[11] R. P. Sear, J. Phys. Cond. Matt. 17, 3997 (2005).
[12] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
[13] P. Rikvold, J. Collins, G. Hansen, and J. Gunton, Surf. Sci. 203
(1988).
[14] S. Manzi, W. Mas, R. Belardinelli, and V. Pereyra, Langmuir
20, 499 (2004).
[15] W. C. K. Poon, F. Renth, R. M. L. Evans, D. J. Fairhurst, M. E.
Cates, and P. N. Pusey, Phys. Rev. Lett. 83, 1239 (1999).
[16] F. Renth, W. C. K. Poon, and R. M. L. Evans, Phys. Rev. E 64,
031402 (2001).
[17] R. M. L. Evans, W. C. K. Poon, and F. Renth, Phys. Rev. E 64,
031403 (2001).
[18] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in
Statistical Physics (Oxford University Press, New York, 1999).
[19] E. N. M. Cirillo and E. Olivieri, J. Stat. Phys. 83, 473 (1996).
[20] T. Fiig, B. M. Gorman, P. A. Rikvold, and M. A. Novotny, Phys.
Rev. E 50, 1930 (1994).
[21] D. P. Sanders, H. Larralde, and F. Leyvraz, (in preparation).
[22] M. A. Novotny, Phys. Rev. Lett. 74, 1 (1995).
[23] P. G. Debenedetti, Metastable Liquids (Princeton University
Press, Princeton, 1996).
[24] R. J. Allen, P. B. Warren, and P. R. ten Wolde, Phys. Rev. Lett.
94, 018104 (2005).
[25] R. J. Allen, D. Frenkel, and P. R. ten Wolde, J. Chem. Phys.
124, 024102 (2006).
[26] R. P. Sear, J. Phys. Chem. B 110, 4985 (2006).
[27] A. J. Page and R. P. Sear, Phys. Rev. Lett. 97, 065701 (2006).
mailto:dsanders@fis.unam.mx
