Approaching equilibrium and the distribution of clusters
Approaching equilibrium and the distribution of clusters
Hui Wang,1 Kipton Barros,2 Harvey Gould,1 and W. Klein2
1Department of Physics, Clark University, Worcester, MA 01610
2Department of Physics and the Center for Computational Science,
Boston University, Boston, MA 02215
Abstract
We investigate the approach to stable and metastable equilibrium in Ising models using a cluster
representation. The distribution of nucleation times is determined using the Metropolis algorithm
and the corresponding φ4 model using Langevin dynamics. We find that the nucleation rate is
suppressed at early times even after global variables such as the magnetization and energy have
apparently reached their time independent values. The mean number of clusters whose size is
comparable to the size of the nucleating droplet becomes time independent at about the same
time that the nucleation rate reaches its constant value. We also find subtle structural differences
between the nucleating droplets formed before and after apparent metastable equilibrium has been
established.
http://arxiv.org/abs/0704.0938v4
I. INTRODUCTION
Understanding nucleation is important in fields as diverse as materials science, biological
physics, and meteorology [1, 2, 3, 4, 5, 6, 7, 8, 9]. Fundamental progress was made when
Gibbs assumed that the nucleating droplet can be considered to be a fluctuation about
metastable equilibrium, and hence the probability of a nucleating droplet is independent of
time [10]. Langer [11] has shown that the probability of a nucleating droplet can be related
to the analytic continuation of the stable state free energy in the limit that the metastable
state lifetime approaches infinity. Hence the assumption by Gibbs is valid in this limit. It has
also been shown that the Gibbs assumption is correct in systems for which the interaction
range R → ∞ [12, 13].
For metastable states with finite lifetimes equilibrium is never reached because a large
enough fluctuation would initiate the transformation to the stable state. However, if the
probability of such a fluctuation is sufficiently small, it is possible that systems investigated
by simulations and experiments can be well approximated as being in equilibrium. Hence,
for metastable lifetimes that are very long, we expect the Gibbs assumption to be a good
approximation.
In practice, nucleation is not usually observed when the lifetime of the metastable state
is very long. Processes such as alloy formation, decay of the false vacuum, and protein
crystallization generally occur during a continuous quench of a control parameter such as
the temperature. It is natural to ask if the nucleation process that is observed occurs when
the system can be reasonably approximated by one in metastable equilibrium. If so, the
nucleation rate will be independent of time.
It is usually assumed that metastable equilibrium is a good approximation when the
mean value of the order parameter and various global quantities are no longer changing
with time. As an example, we consider the nearest-neighbor Ising model on a square lattice
and equilibrate the system at temperature T = 4Tc/9 in a magnetic field h = 0.44. The
relatively small value of the linear dimension L = 200 was chosen in order to avoid nucleation
occurring too quickly. At time t = 0 the sign of the magnetic field is reversed. In Fig. 1 we
plot the evolution of the magnetizationm(t) and the energy e(t) per spin using the Metropolis
algorithm. The solid lines are the fits to an exponential function with the relaxation time
τg ≈ 1.5. In the following we will measure the time in terms of Monte Carlo steps per spin.
(a) m(t). (b) e(t).
FIG. 1: The evolution of the magnetization m(t) and the energy e(t) per spin of the nearest-
neighbor Ising model on a square lattice with linear dimension L = 200 using the Metropolis
algorithm. The system was prepared at temperature T = 4Tc/9 in the external magnetic field
h = 0.44. At time t = 0 the sign of the magnetic field is reversed. The solid lines are fits to
an exponential function with relaxation time τg = 1.5 and 1.2 respectively. (Time is measured in
Monte Carlo steps per spin.) The data is averaged over 5000 runs.
A major goal of our work is to address the question, “Can the system be treated as being
in metastable equilibrium for t >∼ τg?”
If the nucleation rate is independent of time, the probability of a nucleating droplet
occurring at time t after the change of magnetic field is an exponentially decreasing function
of time. To understand this dependence we divide the time into intervals ∆t and write the
probability that the system nucleates in a time interval ∆t as λ∆t, where the nucleation
rate λ is a constant. The probability that nucleation occurs in the time interval (N + 1) is
given by
PN = (1− λ∆t)
Nλ∆t. (1)
If we assume that λ∆t is small and write N = t/∆t, we can write
P (t)∆t = (1− λ∆t)t/∆tλ∆t → e−λtλ∆t, (2)
where P (t)∆t is the probability that the system nucleates at a time between t and t+∆t after
the change of the magnetic field. In the following we ask if the nucleation rate and the mean
values of the order parameter and other thermodynamic quantities become independent
of time at approximately the same time after a quench or is the approach to metastable
equilibrium more complicated?
In Sec. II we determine the probability distribution of the nucleation times and find that
the nucleation rate becomes a constant only after a time τnequil that is much longer than the
relaxation time τg of m(t) and e(t). In Sec. III we study the microscopic behavior of the
system and determine the relaxation time τs for ns, the mean number of clusters of size s,
to approach its equilibrium value [14]. Our main result is that τs is an increasing function of
s, and the time required for ns to reach its equilibrium value is the same order of magnitude
as τnequil for values of s comparable to the nucleating droplet. That is, the time for the
number of clusters that are the size of the nucleating droplet to reach its equilibrium value
is considerably longer than the time for the mean value of the order parameter to become
independent of time within the accuracy that we can determine.
In Secs. IV and V we show that there are subtle differences between the structure of the
nucleating droplets which occur before and after metastable equilibrium appears to have
been achieved. This difference suggests the possibility of finding even greater differences in
the nucleating droplets in systems of physical and technological importance. We summarize
and discuss our results in Sec. VI. In the Appendix we study the evolution of the clusters
after a quench to the critical temperature of the Ising model and again find that that the
clusters equilibrate in size order, with the smaller clusters equilibrating first. Hence in
principle, an infinite system will never equilibrate. How close to equilibrium a system needs
to be and on what spatial scale so that it can be treated by equilibrium methods depends
on the physical process of interest.
II. DISTRIBUTION OF NUCLEATION TIMES
We simulate the Ising model on a square lattice with interaction range R with the Hamil-
tonian
H = −J
<i,j>
sisj − h
si, (3)
where h is the external field. The notation <i, j> in the first sum means that the distance
between spins i and j is within the interaction range R. We studied both nearest-neighbor
(R = 1) and long-range interactions (R ≥ 20). The interaction strength J is scaled as
J = 4/q, where q = 2R(R + 1) is the number of interaction neighbors per spin. The
external field h and the temperature are measured in terms of J . All of our simulations
are at temperature T = 4Tc/9, where Tc is the critical temperature. For R = 1 the critical
temperature is Tc ≈ 2.269. For R >∼ 20 the mean field result Tc = 4 is a good approximation
to the exact value of the critical temperature [21]. As discussed in Sec. I the system is
equilibrated in a magnetic field h. The time t = 0 corresponds to the time immediately after
the magnetic field is reversed.
The clusters in the Ising model are defined rigorously by a mapping of the Ising critical
point onto the percolation transition of a properly chosen percolation model [9, 22, 23]. Two
parallel spins that are within the interaction range R are connected only if there is a bond
between them. The bonds are assigned with the probability pb = 1 − e
−2βJ for R = 1 and
pb = 1− e
−2βJ(1−ρ) near the spinodal, where ρ is the density of the stable spins, and β is the
inverse temperature. Spins that are connected by bonds form a cluster.
Because the intervention method [15] of identifying the nucleating droplet is time con-
suming (see Sec. IV), we use a simpler criterion in this section to estimate the nucleation
time. We monitor the size of the largest cluster (averaged over 20 bond realizations) and
estimate the nucleation time as the time when the largest cluster first reaches a threshold
size s∗. The threshold size s∗ is chosen so that the largest cluster begins to grow rapidly
once its size is greater than or equal to s∗. Because s∗ is larger than the actual size of the nu-
cleating droplet, the nucleation time that we estimate by this criterion will be 1 to 2 Monte
Carlo steps per spin later than the nucleation time determined by the intervention method.
Although the distribution function P (t) is shifted to slightly later times, the nucleation rate
is found to be insensitive to the choice of the threshold.
Figure 2 shows P (t) for R = 1 and h = 0.44, where P (t)∆t is the probability that
nucleation has occurred between time t and t+∆t. The results for P (t) were averaged over
5000 runs. The mean size of the nucleating droplet is estimated to be approximately 25
spins for this value of h. Note that P (t) is an increasing function of t for early times, reaches
a maximum at t = τnequil ≈ 60, and fits to the expected exponential form for t >∼ τnequil.
The fact that P (t) falls below the expected exponential for t < τnequil indicates that the
nucleation rate is reduced from its equilibrium value and that the system is not in metastable
equilibrium. Similar nonequilibrium effects have been observed in Ising-like [16, 17] and
continuous systems [18]. We conclude that the time for the nucleation rate to become
independent of the time after the change of magnetic field is much longer than the relaxation
time τg ≃ 1.5 of the magnetization and energy. We will refer to nucleation that occurs before
metastable equilibrium has been reached as transient nucleation.
(a) P (t). (b) lnP (t).
FIG. 2: The distribution of nucleation times P (t) averaged over 5000 runs for the same system as
in Fig. 1. The threshold size was chosen to be s∗ = 30. (The mean size of the nucleating droplet
is ≈ 25 spins.) (a) P (t) begins to decay exponentially at τnequil ≈ 60. The nucleation rate after
equilibrium has been established is determined from the log-linear plot in (b) and is λ ≈ 9× 10−4
(see Eq. (2)).
In order to see if the same qualitative behavior holds near the pseudospinodal, we simu-
lated the long-range Ising model with R = 20 and h = 1.258. In the mean-field limit R → ∞
the spinodal field is at hs = 1.2704 (for T = 4Tc/9). A plot of m(t) for this system is shown
in Fig. 3(a) and is seen to have the same qualitative behavior as in Fig. 2 for R = 1; the
relaxation time τg ≈ 4.5. In Fig. 3(b) the distribution of nucleation times is shown, and we
see that P (t) does not decay exponentially until t >∼ τnequil = 40. According to Ref. 19, τnequil
should become comparable to τg in the limit R → ∞ because the free energy is described
only by the magnetization in the mean-field limit. We find that the difference between τnequil
and τg is smaller for R = 20 than for R = 1, consistent with Ref. 19.
III. RELAXATION OF CLUSTERS TO METASTABLE EQUILIBRIUM
Given that there is a significant time delay between the relaxation of the magnetization
and the energy and the equilibration of the system as measured by the nucleation rate,
it is interesting to monitor the time-dependence of the cluster-size distribution after the
reverse of the magnetic field. After the change the system gradually relaxes to metastable
equilibrium by forming clusters of spins in the stable direction. How long is required for
the number of clusters of size s to reach equilibrium? In particular, we are interested in the
(a) m(t). (b) ln(P (t)).
FIG. 3: (a) The evolution of m(t) for the long-range Ising model on a square lattice with R = 20,
h = 1.258, and L = 500. The solid line is an exponential fit with the relaxation time τg ≈ 4.5.
The data is averaged over 2000 runs. (b) The distribution of nucleation times P (t) for the same
system and number of runs. P (t) decays exponentially for t >∼ τnequil ≈ 40. The nucleation rate
once equilibrium has been established is λ = 6.4 × 10−2. The mean size of the nucleating droplet
is ≈ 300 spins.
time required for clusters that are comparable in size to the nucleating droplet.
We first consider R = 1 and monitor the number of clusters ns of size s at time t. To
obtain good statistics we chose L = 200 and averaged over 5000 runs. Figure 4 shows the
FIG. 4: The evolution of the number of clusters of size s = 6 averaged over 5000 runs for R = 1
and the same conditions as in Fig. 1. The fit is to the exponential form in Eq. (4) with τs ≈ 8.1
and ns,∞ = 0.0175.
(a) R = 1. (b) R = 20.
FIG. 5: (a) The equilibration time τs as a function of the cluster size s for R = 1 and h = 0.44 the
same conditions as in Fig. 1. The s-dependence of τs is approximately linear. The extrapolated
value of τs corresponding to the mean size of the nucleating droplet (≈ 25 spins) is τextrap ≈ 34,
which is the same order of magnitude as time τnequil ≈ 60 for the system to reach metastable
equilibrium. (b) Log-log plot of the equilibration time τs versus s for R = 20 and h = 1.258 and
the same conditions as in Fig. 3(b). We find that τs ∼ s
x with the exponent x ≈ 0.56. The
extrapolated value of τs corresponding to the mean size of the nucleating droplet (≈ 300 spins)
is τextrap ≈ 30, which is comparable to the time τnequil ≈ 40 for the system to reach metastable
equilibrium.
evolution of n6(t), which can be fitted to the exponential form:
ns(t) = ns,∞[1− e
−t/τs ]. (4)
We find that τs ≈ 8.1 for s = 6. By doing similar fits for a range of s, we find that the time
τs for the mean number of clusters of size s to become time independent increases linearly
with s over the range of s that we can simulate (see Fig. 5). The extrapolated value of τs
corresponding to the mean size of the nucleating droplet (≈ 25 spins by direct simulation)
is τextrap ≈ 34. That is, it takes a time of τextrap ≈ 34 for the mean number of clusters whose
size is the order of the nucleating droplets to become time independent. The time τextrap is
much longer than the relaxation time τg ≈ 1.5 of the macroscopic quantities m(t) and e(t)
and is comparable to the time τnequil ≈ 60 for the nucleation rate to become independent of
time.
Because the number of clusters in the nucleating droplet is relatively small for R = 1
except very close to coexistence (small h), we also consider a long-range Ising model with
R = 20 and h = 1.258 (as in Fig. 3). The relaxation time τs of the clusters near the
pseudospinodal fits to a power law τs ∼ s
x with x ≈ 0.56 (see Fig. 5(b)). We know of no
theoretical explanation for the qualitatively different dependence f the relaxation time τs
on s near coexistence (τs ≃ s) and near the spinodal (τs ≃ s
1/2). If we extrapolate τs to
s = 300, the approximate size of the nucleating droplet, we find that the equilibration time
for clusters of the size of the nucleating droplet is τextrap ≈ 30, which is comparable to the
time τnequil ≈ 40 for the nucleation rate to become independent of time.
To determine if our results are affected by finite size effects, we compared the equilibra-
tion time of the clusters for lattices with linear dimension L = 2000 and L = 5000. The
equilibration times of the clusters were found to be unaffected.
IV. STRUCTURE OF THE NUCLEATING DROPLET
Because nucleation can occur both before and after the system is in metastable equilib-
rium, we ask if there are any structural differences between the nucleating droplets formed
in these two cases. To answer this question, we determine the nature of the nucleating
droplets for the one-dimensional (1D) Ising model where we can make R (and hence the
size of the nucleating droplets) large enough so that the structure of the nucleating droplets
is well defined. In the following we take R = 212 = 4096, h = 1.265, and L = 218. The
relaxation time for m(t) is τg ≈ 40, and the time for the distribution of nucleation times to
reach equilibrium is τnequil ≈ 90.
We use the intervention method to identify nucleation [15]. To implement this method, we
choose a time at which a nucleating droplet might exist and make many copies of the system.
Each copy is restarted using a different random number seed. The idea is to determine if
the largest cluster in each of the copies grows in approximately the same place at about
the same time. If the percentage of copies that grow is greater than 50%, the nucleating
droplet is already in the growth phase; if it is less than 50%, the time chosen is earlier than
nucleation. We used a total of 20 trials to make this determination.
Our procedure is to observe the system for a time tobs after the intervention and determine
if the size of the largest cluster exceeds the threshold size s∗ at approximately the same
location. To ensure that the largest cluster at tobs is the same cluster as the original one, we
require that the center of mass of the largest cluster be within a distance r∗ of the largest
cluster in the original configuration. If these conditions are satisfied, the nucleating droplet
is said to grow. We choose tobs = 6, r
∗ = 2R, and s∗ = 2000. (In comparison, the size of the
nucleating droplet for the particular run that we will discuss is ≈ 1080 spins.)
There is some ambiguity in our identification of the nucleation time because the saddle
point parameter is large but finite [9]. This ambiguity manifests itself in the somewhat
arbitrary choices of the parameters tobs, r
∗, and s∗. We tried different values for tobs, r
∗, and
s∗ and found that our results depend more strongly on the value of the parameter r∗ than on
the values of tobs and s
∗. If we take r∗ = R/2, the nucleating droplets almost always occur
one to two Monte Carlo steps per spin later than for r∗ = 2R. The reason is that the linear
size of the nucleating droplet is typically 6 to 8R, and its center of mass might shift more
than R/2 during the time tobs. If such a shift occurs, a cluster that would be said to grow for
r∗ = 2R would not be counted as such because it did not satisfy the center of mass criterion.
This shift causes an overestimate of the time of the nucleating droplet. A reasonable choice
of r∗ is 20% to 40% of the linear size of the nucleating droplet. The choice of parameters is
particularly important here because the rate of growth of the transient nucleating droplets
is slower than the growth rate of droplets formed after metastable equilibrium has been
reached. Hence, we have to identify the nucleating droplet as carefully as possible.
Because nucleation studies are computationally intensive, we used a novel algorithm for
simulating Ising models with a uniform long-range interaction [20]. The algorithm uses a
hierarchical data structure to store the magnetization at many length scales, and can find
the energy cost of flipping a spin in time O((lnR)d), rather than the usual time O(Rd),
where d is the spatial dimension.
Figure 6 shows the fraction of copies for which the largest cluster grows as a function of
the intervention time. For this particular run the nucleating droplet is found to occur at
t ≈ 37.4.
We simulated 100 systems in which nucleation occurred before global quantities such
as m(t) became independent of time, t < τg ≈ 40, and 100 systems for which nucleation
occurred after the nucleation rate became time independent (t > τnequil ≈ 90). We found
that the mean size of the nucleating droplet for t < τg is ≈ 1200 with a standard deviation
of σ ≈ 150 in comparison to the mean size of the nucleating droplet for t > τnequil of ≈ 1270
and σ ≈ 200. That is, the nucleating droplets formed before metastable equilibrium has
been reached are somewhat smaller.
We introduce the cluster profile ρcl to characterize the shape of the largest cluster at the
time of nucleation. For a particular bond realization a spin that is in the stable direction
might or might not be a part the largest cluster due to the probabilistic nature of the bonds.
For this reason bond averaging is implemented by placing 100 independent sets of bonds
between spins with probability pb = 1 − e
−2βJ(1−ρ) in the stable direction. The clusters are
identified for each set of bonds, and the probability pi that spin i is in the largest cluster is
determined. The values of pi for the spins in a particular bin are then averaged using a bin
width equal to R/4. This mean value of pi is associated with ρcl. Note that the spins that
point in the unstable direction are omitted in this procedure. The mean cluster profile is
found by translating the peak position of each droplet to the origin.
Figure 7(a) shows the mean cluster profile formed after metastable equilibrium has been
established (t > τnequil ≈ 90). The position x is measured in units of R. For comparison we
fit ρcl to the form
ρ(x) = A sech2(x/w) + ρ0, (5)
with Acl = 0.36, wcl = 2.95 and ρ0 = 0 by construction. In Fig. 7(b) we show a comparison
of ρcl to the Gaussian form Ag exp(−(x/wg)
2) with Ag = 0.35 and wg = 3.31. Note that
Eq. (5) gives a better fit than a Gaussian, which underestimates the peak at x = 0 and the
FIG. 6: The fraction of copies for which the largest cluster grows for a particular run for a 1D
Ising model with R = 212, h = 1.265, and L = 218. The time for 50% growth is ≈ 37.4. The
largest cluster at this time corresponds to the nucleating droplet and has ≈ 1080 spins. For this
intervention 100 copies were considered; twenty copies were considered for all other runs.
(a) Comparison to Eq. (5). (b) Comparison to Gaussian.
FIG. 7: Comparison of the mean cluster profile (•) in the 1D Ising model after metastable equi-
librium has been established with (a) the form in Eq. (5) and (b) a Gaussian. Note that Eq. (5)
gives a better fit than the Gaussian, which underestimates the peak at x = 0 and the wings. The
x axis is measured in units of R.
wings. Although Unger and Klein [12] derived Eq. (5) for the magnetization saddle point
profile, we see that this form also provides a good description of the cluster profile.
A comparison of the cluster profiles formed before and after metastable equilibrium is
shown in Fig. 8. Although both profiles are consistent with the form in Eq. (5), the transient
nucleating droplets are more compact, in agreement with the predictions in Ref. 19.
We also directly coarse grained the spins at the time of nucleation to obtain the density
profile of the coarse-grained magnetization ρm(x) (see Fig. 9(a)). The agreement between
the simulation and analytical results [24] are impressive, especially considering that the
analytical form is valid only in the limit R → ∞. The same qualitative differences between
the nucleating droplets that occur before and after metastable equilibrium is found (see
Fig. 9(b)), although the magnetization density profile is much noisier than that based on
the cluster analysis.
V. LANGEVIN SIMULATIONS
It is interesting to compare the results for the Ising model and the Langevin dynamics
of the φ4 model. One advantage of studying the Langevin dynamics of the φ4 theory is
that it enables the efficient simulation of systems with a very large interaction range R. If
FIG. 8: The cluster profiles of the nucleating droplets formed before (dashed line) and after (solid
line) metastable equilibrium has been established. Both profiles are consistent with the form given
in Eq. (5), but the transient nucleating droplets are slightly more compact. The fitting parameters
are A = 0.38 and w = 2.67 for the transient droplets and A = 0.35 and w = 2.95 for the droplets
formed after the nucleation rate has become independent of time.
all lengths are scaled by a large value of R, the effective magnitude of the noise decreases,
making faster simulations possible.
The coarse grained Hamiltonian analogous to the 1D ferromagnetic Ising model with
long-range interactions in an external field h can be expressed as
H [φ] = −
+ ǫφ2 + uφ4 − hφ, (6)
where φ(x) is the coarse-grained magnetization. A dynamics consistent with this Hamilto-
nian is given by,
+ η = −M
+ 2εφ+ 4uφ3 − h
+ η, (7)
where M is the mobility and η(x, t) represents zero-mean Gaussian noise with
〈η(x, t)η(x′, t′)〉 = 2kTMδ(x− x′)δ(t− t′).
For nucleation near the spinodal the potential V = εφ2+ uφ4− hφ has a metastable well
only for ε < 0. The magnitude of φ and h at the spinodal are given by hs =
(8|ε|3/27u)
and φs =
(|ε|/6u), and are found by setting V ′ = V ′′ = 0. The distance from the spinodal
is characterized by the parameter ∆h = |hs − h|. For ∆h/hs ≪ 1, the bottom of the
metastable well φmin is near φs, specifically φmin = −φs(1 +
2∆h/3hs).
(a) Comparison with Eq. (5). (b) Comparison of profiles.
FIG. 9: (a) The magnetization density profile of the nucleating droplets formed after metastable
equilibrium has been established. The solid line is the analytical solution [24] which has the form
in Eq. (5) with the calculated values A = 0.085, w = 2.65, and ρ0 = −0.774. (b) Comparison of the
density profile of nucleating droplets formed before (dashed line) and after (solid line) metastable
equilibrium has been established by coarse graining the magnetization. The same qualitative
differences between the nucleating droplets that occur before and after metastable equilibrium are
observed as in Fig. 8, although the magnetization density profile is much noisier than the cluster
density profile.
The stationary solutions of the dynamics are found by setting δH/δφ = 0. Besides
the two uniform solutions corresponding to the minima in V , there is a single nonuniform
solution which approximates the nucleating droplet profile when the nucleation barrier is
large. When ∆h/hs ≪ 1, the profile of the nucleating droplet is described by Eq. (5) with
hs/6∆h/φs, w = (8hs∆hφ
−1/4, and ρ0 = φmin [12].
The dynamics (7) is numerically integrated using the scheme [25]
φ(t+∆t) = φ(t)−∆tM
+ 2εφ+ 4uφ3 − h
η, (8)
where d2φ/dx2 is replaced by its central difference approximation. Numerical stability re-
quires that ∆t < (∆x/R)2, but it is often desirable to choose ∆t even smaller for accuracy.
As for the Ising simulations, we first prepare an equilibrated system with φ in the stable
well corresponding to the direction of the external field h. At t = 0 the external field is
reversed so that the system relaxes to metastable equilibrium. We choose M = 1, T = 1,
ε = −1, u = 1, and ∆h = 0.005. The scaled length of the system is chosen to be L/R = 300.
FIG. 10: Log-linear plot of the distribution P (t) of nucleation times for the one-dimensional
Langevin equation with R = 2000 (×) and R = 2500 (•) averaged over 50,000 runs. The distribu-
tion is not exponential for early times, indicating that the system is not in metastable equilibrium.
Note that the nucleation rate is a rapidly decreasing function of R.
We choose R to be large so that, on length scales of R, the metastable φ fluctuates near its
equilibrium value φmin ≈ −0.44.
After nucleation occurs φ will rapidly grow toward the stable well. To determine the
distribution of nucleation times, we assume that when the value of the field φ in any bin
reaches 0, nucleation has occurred. This relatively crude criterion is sufficient for determin-
ing the distribution of nucleation times if we assume that the time difference between the
nucleation event and its later detection takes a consistent value between runs.
Figure 10 compares the distribution of 50,000 nucleation times for systems with R = 2000
and R = 2500 with ∆x/R = 1 and ∆t = 0.1. The distribution shows the same qualitative
behavior as found in the Metropolis simulations of the Ising model (see Fig. 2). For example,
the distribution of nucleation times is not exponential for early times after the quench. As
expected, the nucleation rate decreases as R increases. Smaller values of ∆x and ∆t give
similar results for the distribution.
To find the droplet profiles, we need to identify the time of nucleation more precisely. The
intervention criterion, which was applied in Sec. IV, is one possible method. In the Langevin
context we can employ a simpler criterion: nucleation is considered to have occurred if φ
decays to the saddle-point profile (given by Eq. (5) for ∆h/hs ≪ 1) when φ is evolved using
FIG. 11: Comparison of the density profile φ(x) of the nucleating droplets found by numerically
solving the Langevin equation after metastable equilibrium has been reached for R = 2000 (×)
and R = 4000 (•) to the theoretical prediction (solid line) from Eq. (5) using the calculated values
A = 0.096, w = 3.58, and ρ0 = −0.44. The numerical solutions are averaged over 1000 profiles.
The results suggest that as R increases, the observed nucleation profiles converge to the prediction
of mean-field theory.
noiseless dynamics [19, 26]. For fixed ∆h these two criteria agree in the R → ∞ limit, but
can give different results for finite R [27].
In Fig. 11 we plot the average of 1,000 density profiles of the nucleating droplets formed
after metastable equilibrium has been established for R = 2000 and R = 4000. Note that
there are noticeable deviations of the averaged profiles from the theoretical prediction in
Eq. (5), but the deviation is less for R = 4000. The deviation is due to the fact that the
bottom of the free energy well in the metastable state is skewed; a similar deviation was also
observed in the Ising model. We also note that the individual nucleating droplets look much
different from their average. It is expected that as R increases, the profiles of the individual
nucleating droplets will converge to the form given by Eq. (5).
In Fig. 12 we compare the average of 1,000 density profiles of nucleating droplets before
and after metastable equilibrium has been established. As for the Ising model, there are
subtle differences consistent with the predictions of Ref. 19. The transient droplets have
slightly lower background magnetization and compensate by being denser and more compact.
FIG. 12: The density profile of the nucleating droplets found from numerical solutions of the
Langevin equation formed before (dotted line) and after (solid line) metastable equilibrium has
been established. Nucleation events occurring before t = 15 are transient, and events occurring
for t ≥ 30 are metastable. Both plots are the result of 1000 averaged profiles with an interaction
range R = 2000.
VI. SUMMARY
Although the time-independence of the mean values of macroscopic quantities such as
the magnetization and the energy is often used as an indicator of metastable equilibrium, we
find that the observed relaxation time of the clusters is much longer for sizes comparable to
the nucleating droplet. This longer relaxation time explains the measured non-constant nu-
cleation rate even when global quantities such as the magnetization appear to be stationary.
By identifying the nucleating droplets in the one-dimensional long-range Ising model and the
Langevin equation, we find structural differences between the nucleating droplets which oc-
cur before and after metastable equilibrium has been reached. Our results suggest that using
global quantities as indicators for metastable equilibrium may not be appropriate in general,
and distinguishing between equilibrium and transient nucleation is important in studying
the structure of nucleating droplets. Further studies of transient nucleation in continuous
models of more realistic systems would be of interesting and practical importance.
Finally, we note a subtle implication of our results. For a system to be truly in equilibrium
would require that the mean number of clusters of all sizes be independent of time. The
larger the cluster, the longer the time that would be required for the mean number to
become time independent. Hence, the bigger the system, the longer the time that would
(a) R = 1. (b) R = 128.
FIG. 13: The relaxation of the magnetization m(t) of the 2D Ising model at T = Tc starting from
T0 = 0. (a) R = 1, Tc = 2.269, L = 5000. (b) R = 128, Tc = 4, L = 1024. The straight line is the
fit to a power law with slope ≈ 0.057 for R = 1 and slope ≈ 0.51 for R = 128.
be required for the system to reach equilibrium. Given that the system is never truly
in metastable equilibrium so that the ideas of Gibbs, Langer, and others are never exactly
applicable, when is the system close enough to equilibrium so that any possible simulation or
experiment cannot detect the difference? We have found that the magnetization and energy
are not sufficient indicators for nucleation and that the answer depends on the process being
studied. For nucleation the equilibration of the number of clusters whose size is comparable
to the size of the nucleating droplet is the relevant indicator.
APPENDIX A: RELAXATION OF CLUSTERS AT THE CRITICAL TEMPER-
ATURE
Accurate determinations of the dynamical critical exponent z have been found from the
relaxation of the magnetization and energy at the critical temperature. In the following
we take a closer look at the relaxation of the Ising model by studying the approach to
equilibrium of the distribution of clusters of various sizes.
We consider the Ising model on a square lattice with L = 5000. The system is initially
equilibrated at either zero temperature T0 = 0 (all spins up) or at T0 = ∞, and then
instantaneously quenched to the critical temperature Tc. The Metropolis algorithm is used.
As a check on our results we first determine m(t) starting from T0 = 0. Scaling arguments
FIG. 14: The evolution of the number of clusters of size s = 100 at T = Tc starting from T0 = 0.
The fit to Eq, (A2) gives ns,∞ = 51.3, C1 = −42, C2 = −15, τ1 = 156, and τ2 = 1070.
suggest that m(t) approaches its equilibrium value as [28]
f(t) = Bt−β/νz + f∞, (A1)
where the static critical exponents are β = 1/8 and ν = 1 for finite R and β = 1/2 and
ν = 1/2 in the mean-field limit. The fit of our results in Fig. 13 to Eq. (A1) yields the
estimate z ≈ 2.19 for R = 1 and z ≈ 1.96 for R = 128, which are consistent with previous
results [29, 30]. Note that no time scale is associated with the evolution of m(t).
We next determined ns(t), the number of clusters of size s at time t after the temperature
quench. Because all the spins are up at t = 0, the number of (down) clusters of size s begins
at zero and increases to its (apparent) equilibrium value ns,∞. The value of the latter
depends on the size of the system.
Figure 14 shows the evolution of clusters of size s = 100 for one run. Because we know
of no argument for the time dependence of ns(t)− ns,∞ except in the mean-field limit [30],
we have to rely on empirical fits. We find that the time-dependence of ns(t) can be fitted
to the sum of two exponentials,
ns(t)− ns,∞ = C1e
−t/τ1 + C2e
−t/τ2 , (A2)
where C1, C2, τ1, and τ2 are parameters to be fitted with τ2 > τ1.
Figure 15(a) shows the relaxation time τ2 as a function of s for R = 1 at T = Tc starting
from T0 = 0. Note that the bigger the cluster, the longer it takes to reach its equilibrium
distribution. That is, small clusters form first, and larger clusters are formed by the merging
(a) T0 = 0. (b) T0 = ∞.
FIG. 15: The relaxation time τ2 versus the cluster size s at T = Tc for R = 1 starting from (a)
T0 = 0 and (b) T0 = ∞. The log-log plot in (a) yields τ2 ∼ s
(a) s = 30. (b) s = 3000.
FIG. 16: The time dependence of the number of clusters of size s = 30 and s = 3000 at T = Tc
for R = 1 starting from T0 = ∞. Note that ns=30 monotonically decreases to its equilibrium value
and ns=1000 overshoots its equilibrium value. (a) C1 = 2367, C2 = 332, ns=30,∞ = 738, τ1 = 16,
and τ2 = 403. (b) C1 = −0.42, C2 = 0.22, ns=3000,∞ = 0.11, τ1 = 130, and τ2 = 1290.
of smaller ones. The s-dependence of τ2 can be approximately fitted to a power law with
the exponent 0.4.
To prepare a configuration at T0 = ∞, the system is randomized with approximately half
of the spins up and half of the spins down. The temperature is instantaneously changed to
T = Tc. As before, we focus on the relaxation of down spin clusters. In contrast to the T0 = 0
case, the evolution of the clusters falls into three classes (see Fig. 16). For small clusters
(1 ≤ s ≤ 40), ns monotonically decreases to its equilibrium value. This behavior occurs
because the initial random configuration has an abundance of small clusters so that lowering
the temperature causes the small clusters to merge to form bigger ones. For intermediate
size clusters (40 < s < 4000), ns first increases and then decreases to its equilibrium value.
The initial growth is due to the rapid coalescence of smaller clusters to form intermediate
ones. After there are enough intermediate clusters, they slowly coalesce to form bigger
clusters, which causes the decrease. For clusters with s > 4000, ns slowly increases to its
equilibrium value. The range of sizes for these different classes of behavior depends on the
system size. In all three cases ns(t) can be fitted to the sum of two exponentials. One of
the two coefficients is negative for 40 < s < 4000 for which ns(t) overshoots its equilibrium
value. The relaxation time τ2 is plotted in Fig. 15(b) as a function of s.
ACKNOWLEDGMENTS
We thank Aaron O. Schweiger for very useful discussions. Bill Klein acknowledges the
support of Department of Energy grant # DE-FG02-95ER14498 and Kipton Barros was
supported in part by the National Science Foundation grant # DGE-0221680. Hui Wang
was supported in part by NSF grant # DUE-0442581. The simulations at Clark University
were done with the partial support of NSF grant # DBI-0320875.
[1] Dimo Kashchiev, Nucleation: Basic Theory with Applications (Butterworths-Heinemann, Ox-
ford, 2000).
[2] N. E. Chayen, E. Saridakis, R. El-Bahar, and Y. J. Nemirovsky, Mol. Biol. 312 (4), 591 (2001).
[3] N. Wang, Y. H. Tang, Y. F. Zhang, and C. S. Lee, Phys. Rev. B 58, R16 024 (1998).
[4] S. Auer and D. Frenkel, Nature 413, 711 (2001).
[5] N. Delgehyr, J. Sillibourne, and M. Bornens, J. Cell Science 118, 1565 (2005).
[6] E. Pechkova and C. Nicolini, J. Cell Biochem. 85 (2), 243 (2002).
[7] A. R. Fersht, Proc. Natl. Acad. Sci. U. S. A. 92, 10869 (1995).
[8] D. Stauffer, A. Coniglio, and D. W. Heermann, Phys. Rev. Lett. 49, 1299 (1982).
[9] W. Klein, H. Gould, N. Gulbahce, J. B. Rundle, and K. Tiampo, Phys. Rev. E 75, 031114
(2007).
[10] J. D. Gunton, M. san Miguel, and P. Sahni, in Phase Transitions and Critical Phenomena,
edited by C. Domb and J. L. Lebowitz (Academic Press, New York, 1983), Vol. 8.
[11] J. S. Langer, Ann. Phys. (NY) 41, 108 (1967).
[12] C. Unger and W. Klein, Phys. Rev. B 29, 2698 (1984).
[13] More precisely, the Gibbs assumption is correct in systems for which the interaction range
R ≫ 1 and which are not too close to the spinodal. See Ref. 9 for more details.
[14] We assume that when the nucleation rate and mean number of clusters of a given size become
apparently independent of time that they have reached their equilibrium values.
[15] L. Monette, W. Klein, and M. Zuckermann, J. Stat. Phys. 66, 117 (1992).
[16] D. W. Heermann and C. Cordeiro, Int. J. Mod. Phys. 13, 1419 (2003).
[17] K. Brendel, G. T. Barkema, and H. van Beijeren, Phys. Rev. B 71, 031601 (2005).
[18] H. Huitema, J. van der Eerden, J. Janssen, and H. Human, Phys. Rev. B 62, 14690 (2000).
[19] A. O. Schweiger, K. Barros, and W. Klein, Phys. Rev. E 75, 031102 (2007).
[20] K. Barros, manuscript in preparation.
[21] E. Luijten, H. W. J. Blöte, and K. Binder, Phys. Rev. E 54, 4626 (1996).
[22] W. Klein in Computer Simulation Studies in Condensed Matter Physics III, edited by D. P.
Landau, K. K. Mon, and H. B. Schuttler (Springer-Verlag, Berlin, Heidelberg, 1991).
[23] A. Coniglio and W. Klein, J. Phys. A 13, 2775 (1980).
[24] The density profile of the nucleating droplet of the Ising model has been calculated analytically
in the limit R → ∞ (K. Barros, unpublished). The result is consistent with the form in
Eq. (5) with A = 3(βJ)−1
∆h/φs, w = 3
−1/2φ
−1/4, and ρ0 = −φs − A/3, where
1− (βJ)−1 is the magnitude of φ at the spinodal. From this analytical solution, the
calculated parameters are found to be A = 0.085, w = 2.65, ρ0 = −0.774 which are very close
to the values fitted to the simulation data, A = 0.084, w = 2.45, ρ0 = −0.764.
[25] J. G. Gaines, in Stochastic Partial Differential Equations, edited by A. M. Etheridge (Cam-
bridge University Press, Cambridge, 1995), pp. 55–71.
[26] A. Roy, J. M. Rickman, J. D. Gunton, and K. R. Elder, Phys. Rev. E 57, 2610 (1998).
[27] Consider a fluctuation that decays to the metastable phase under noiseless dynamics. To
perform the intervention method we make many copies of the configuration and examine the
percentage that grow after a given waiting time. Although the expected drift is a decay to the
metastable phase, every copy has time to sample a path in configuration space. It is possible
that during this waiting time the majority of copies discover and grow toward the stable phase,
contradicting the result from the zero-noise criterion. However, for R ≫ 1 the sampling path
will be dominated by the drift term and the two nucleation criteria agree.
[28] M. Suzuki. Phys. Lett. A 58, 435 (1976) and M. Suzuki. Prog. Theor. Phys. 58, 1142 (1977).
See also A. Linke, D. W. Heermann, P. Altevogt, and M. Siegert, Physica A 222, 205 (1995).
[29] M. Nightingale and H. Blöte, Phys. Rev. B 62, 1089 (2000).
[30] L. Colonna-Romano, A. I. Mel’cuk, H. Gould, and W. Klein, Physica A 209, 396 (1994).
	Introduction
	Distribution of nucleation times
	Relaxation of clusters to metastable equilibrium
	Structure of the nucleating droplet
	Langevin simulations
	Summary
	Relaxation of clusters at the critical temperature
	Acknowledgments
	References
