Metastability and the Ising Model

Roberto H. Schonmann

Abstract. We present recent results on a classical non-equilibrium statistical mechanics problem, in the context of a well-studied idealized interacting particle system, called kinetic Ising model. The problem concerns the speed and the patterns of relaxation of statistical mechanical systems in the proximity of the phase-transition region, and is related to the problem of understanding the metastable behavior of systems in such regions.

1991 Mathematics Subject Classification: 60K35 82B27

Keywords and Phrases: kinetic Ising model, stochastic Ising model, Glauber dynamics, metastability, relaxation, nucleation, droplet growth, Wulff shape, large deviations, asymptotic expansion

It is well known that a ferromagnetic material which is in equilibrium under a negative external magnetic field relaxes to equilibrium very slowly after the magnetic field is switched to a small positive value. A detailed mathematical analysis of such a phenomenon can only be performed on simplified models. It is widely accepted that an appropriate model for this and many other purposes is a kinetic Ising model: a Markov process which endows the traditional Ising model with a particular stochastic dynamics. On each vertex of an infinite lattice ℤ^{d}, we have variables (called spins) which take the values −1 or +1. The system evolves in continuous time as a Markov process which is time-reversible and has as invariant measures the classical Gibbs measures of statistical mechanics. When the temperature parameter, T, appearing in the definition of the model is small enough, there is a phase transition which takes place when the external field parameter, h, changes sign (this corresponds to the change from a negative to a positive orientation of most spins). The question then arises of how the system relaxes to equilibrium when h is small and positive, and the system is initially in an equilibrium distribution corresponding to a small negative value of h.

Simulations have shown that the relaxation mechanism is driven by the behavior of the clusters (droplets) of +1-spins which form initially in the sea of −1-spins. While small droplets tend to shrink and disappear, large ones tend to grow and are responsible for the relaxation. This phenomenon has long been understood on non-rigorous heuristic grounds, and can be used to predict for instance the order of magnitude, as h ↘ 0, of the relaxation time for the process. The prediction is that the relaxation time grows as exp(λh^{d−1}), where λ is a constant which can be computed. The value of λ is, in particular, related to the equilibrium surface tension of the Ising model through the Wulff construction, which solves a variational problem.

In this note we will overview rigorous results of the type described above and also some important extensions. A thorough review of metastability, even in the context of the kinetic Ising models, is far beyond the scope of this note. Here we will limit ourselves to the main results in the papers [Sch1] and [SS], which concern metastability in the vicinity of the phase transition region. A great deal of recent progress on metastability of the kinetic Ising models stems from the fact that these models also display metastable behavior away from this region, at low enough temperature. For a detailed discussion of relations between the various manifestations of metastability of kinetic Ising models we refer the reader to [Sch2], where further reference to the literature can also be found. More recent progress in this direction is contained in [Nev], [BC], [DS], [CO], [CL] and references therein. For a paper which reports on extensive numerical studies directly related to the mathematical work reviewed here, we refer the reader to [RTMS].

The precise definition of the kinetic Ising models is lengthy and somewhat technical. It can be found, e.g., in [Sch1] and [SS]. For the purpose of this note it is best to just give a somewhat intuitive description. At each site of the lattice ℤ^{d} there is a variable (spin) which can take the value −1 or +1. The configuration on the complete lattice is then an element of the space Ω = {−1, +1}^{ℤ^{d}}. The system evolves in time, with spins flipping back and forth, at rates which depend on the state of nearby spins. The system as a whole is a Markov process with state space Ω. The interaction among spins is driven by an energy function (Hamiltonian) formally defined on Ω by

H_{h}(σ) = −\frac{1}{2} ∑_{x,y n.n.} σ(x)σ(y) − \frac{h}{2} ∑_{x} σ(x),

where “x, y n.n.” means that x and y are nearest neighbors in ℤ^{d}, i.e., they are separated by Euclidean distance 1, h ∈ ℝ is the external field and σ ∈ Ω is a generic configuration.

Formally, Gibbs distributions are defined as probability distributions μ over Ω, with

μ(σ) = \frac{exp(−H_{h}(σ)/T)}{Normalization},

where T = 1/β > 0 is the temperature. When h ≠ 0 or T > T_{c} = T_{c}(d) it is known that there is a unique Gibbs distribution, which then describes the system in equilibrium and will be denoted by μ_{T,h}. In d = 1, T_{c} = 0, but for d ≥ 2, T_{c} > 0. The segment {{0} × (0, T_{c})} of the phase diagram h × T corresponds then to the phase-transition region. For these values of the parameters there are multiple Gibbs distributions; one of them corresponds to a limit of Gibbs distributions under h < 0 (resp. h > 0) as h ↗ 0 (resp. h ↘ 0), and is called the (−)-phase (resp. the (+)-phase), represented by μ_{T,−} (resp. μ_{T,+}). Expectations with respect to Gibbs measures will be denoted in the standard fashion

⟨f⟩_{T,h} = ∫ fdμ_{T,h}.

Of particular interest is the magnetization m(T, h) = ⟨σ(0)⟩_{T,h}. Away from the phase-transition region, m(T, ·) is analytic. It is nevertheless believed that for T < T_{c} this function has no analytic continuation from h < 0 to h > 0. This result has been proved indeed for low enough T in [Isa].

The time evolution which defines the kinetic Ising model as a Markov process on Ω is given by a generator L of the form given next. Intuitively, when the configuration is σ, the spin at each site x ∈ ℤ^{d} is flipping at a rate c(x, σ).

(Lf)(σ) = ∑_{x∈ℤ^{d}} c(x, σ)(f(σ^{x}) − f(σ)).

Here f : Ω → ℝ is supposed to be a local observable, i.e., to depend only on the spin at finitely many sites of the lattice, σ^{x} is the configuration obtained from σ by flipping the spin at the site x, and c(x, σ) is called the rate of flip of the spin at the site x when the system is in the state σ. The rates c(x, σ) are supposed to satisfy certain conditions, the main one of them being called detailed balance or reversibility, and formally given by

μ(σ)c(x, σ) = μ(σ^{x})c(x, σ^{x}).

This assures that the Gibbs distributions are invariant for the process. Other conditions are that the rates are invariant under translations of the lattice, are of finite range of dependency, are monotone in the configuration and external field, and are uniformly bounded above and below when T is fixed and |h| is small. Several choices can be made for the rates, satisfying all this conditions. To give a few examples, we introduce

Δ_{x}H_{h}(σ) = H_{h}(σ^{x}) − H_{h}(σ).

Common choices of rates are: Example 1) Metropolis Dynamics

c_{T,h}(x, σ) = exp(−β(Δ_{x}H_{h}(σ))^{+}),

where (a)^{+} = max{a, 0} is the positive part of a.

Example 2) Heat Bath Dynamics

c_{T,h}(x, σ) = \frac{1}{1 + exp(βΔ_{x}H_{h}(σ))}.

Example 3)

c_{T,h}(x, σ) = exp (−\frac{β}{2} Δ_{x}H_{h}(σ)) .

If in the kinetic Ising model the initial configuration is selected at random according to a probability measure ν, then the resulting process is denoted by (σ_{T,h;t}^{ν})_{t≥0}. When ν is concentrated on the configuration with all spins −1, we will denote this process by (σ_{T,h;t}^{−})_{t≥0}. The probability measure on the space of trajectories of the process will be denoted by ℙ, and the corresponding expectation by 𝔼.

The following is the main result of [Sch1].

Theorem 1. For each dimension d ≥ 2 there is T_{0} > 0 such that for every temperature T ∈ (0, T_{0}) the following happens. There are constants 0 < λ_{1}(T) ≤ λ_{2}(T) < ∞ such that if we let h ↘ 0 and t → ∞ together, then for every local observable f

i) 𝔼(f(σ_{T,h;t}^{−})) → ⟨f⟩_{T,−} if lim sup h^{d−1} log t < λ_{1}(T). ii) 𝔼(f(σ_{T,h;t}^{−})) → ⟨f⟩_{T,+} if lim inf h^{d−1} log t > λ_{2}(T).

Explicit estimates on the values of λ_{1}(T) and λ_{2}(T) were also given in [Sch1]. The theorem above was conjectured by Aizenman and Lebowitz in [AL], where they proved a similar result for certain deterministic cellular automata evolving from initial random configurations selected according to translation invariant product measures. Actually they conjectured the stronger result, which states that also λ_{1}(T) = λ_{2}(T) =: λ_{c}(T).

Theorem 1 was greatly improved in [SS] in the case in which d = 2. In particular in this paper the conjecture by Aizenman and Lebowitz was fully vindicated in this case. A somewhat simplified and partial statement of the main result in [SS] is as follows.

Theorem 2. Suppose d = 2 and T < T_{c}. There is a constant λ_{c} = λ_{c}(T) such that for every probability distribution ν = μ_{T,h′}, h′ < 0, the following happens.

i) If 0 < λ < λ_{c}, then for each n ∈ {1, 2, ...} and for each local observable f,

𝔼 (f (σ_{T,h;exp(λ/h)}^{ν})) = ∑_{j=0}^{n−1} \frac{1}{j!} \frac{d^{j}⟨f⟩_{T,ĥ}}{dĥ^{j}} |_{ĥ=0_{−}} h^{j} + O(h^{n})

for h > 0, where O(h^{n}) is a function of f and h which satisfies lim sup_{h↘0} |O(h^{n})|/h^{n} < ∞.

ii) If λ > λ_{c}, then for any finite positive C there is a finite positive C_{1} such that for every local observable f,

|𝔼 (f (σ_{T,h;exp(λ/h)}^{ν})) − ⟨f⟩_{T,h}| ≤ C_{1} ||f||_{∞} exp (−\frac{C}{h} ) ,

for all h > 0.

The value of λ_{c}(T) can be written in terms of other quantities which are related to the equilibrium distributions of the Ising model. This expression and its meaning, which are of great relevance, will be presented later in the paper. Next we compare Theorems 1 and 2 and explain some of their content.

Three of the ways in which Theorem 2 improves on Theorem 1 when d = 2 are: 1) There is a single constant λ_{c} separating the regimes (i) and (ii). 2) The temperature is now only required to be below T_{c}. 3) The initial distribution is much more general than in Theorem 1, where it was supposed to be concentrated on the configuration with all spins down. It is natural indeed to start from an equilibrium state at a small negative h, change it to a small positive h and observe the evolution of the system afterwards.

To illustrate and clarify the main way in which Theorem 2 improves further the statement in Theorem 1, let us take the local observable given by f(σ) = σ(0) and n = 2. We have then, when 0 < λ < λ_{c}

𝔼 (σ_{T,h;exp(λ/h)}^{ν}(0)) = −m^{∗} + χh + O(h^{2}),

when h > 0. Here

m^{∗} = m^{∗}(T) = ⟨σ(0)⟩_{T,+} = −⟨σ(0)⟩_{T,−},

is the spontaneous magnetization, and

χ = χ(T) = \frac{d⟨σ(0)⟩_{T,h}}{dh} |_{h=0_{−}} =

(\frac{β}{2} ) ∑_{x∈ℤ^{2}} \{⟨σ(0)σ(x)⟩_{T,−} − ⟨σ(0)⟩_{T,−}⟨σ(x)⟩_{T,−}\} ,

is the susceptibility at h = 0_{−}. This means that when h > 0 is small the function −m^{∗} + χh is a better approximation to 𝔼 (σ_{T,h;exp(λ/h)}^{ν}(0)) than the constant function identical to −m^{∗} = ⟨f⟩_{T,−}. This function −m^{∗} + χh is the smooth linear continuation into the region h ≥ 0 of the function which to h < 0 associates the equilibrium expectation ⟨f⟩_{T,h}. Similar interpretations can be given for larger values of n and arbitrary f. In this sense Theorem 2 shows that the dynamics gives meaning to arbitrarily smooth metastable continuations of the distributions μ_{T,h}, h < 0, into the region h > 0, inspite of the absence of an analytic continuation.

In the Physics literature (see, e.g., [BM]), one sometimes relates the metastable relaxation of a system to the presence of a “plateau” in the graph corresponding to the time evolution of a quantity of the type of 𝔼 (f (σ_{T,h;t}^{ν})). Of course, strictly speaking there is no “plateau”, and generically the slope of such a function is never 0. Still, from the experimental point of view a rough “plateau” can be seen and described as follows. In a relatively short time 𝔼 (f (σ_{T,h;t}^{ν})) seems to converge to a value close to ⟨f⟩_{T,−}; after this, one sees an apparent flatness in the relaxation curve over a period of time which may be quite long compared with the time needed to first approach this value. But eventually the relaxation curve starts to deviate from this almost constant value and move towards the true asymptotic limit, close to ⟨f⟩_{T,+}. The experimentally almost flat portion of the relaxation curve is referred to as a “plateau”. Theorem 2 can be seen to some extent as giving some precise meaning to such a “plateau”, and we discuss now two ways in which this can be done. First note that if 0 < λ′ < λ′′ < λ_{c}, then from Part (i) of the Theorem we have

𝔼 (f (σ_{T,h;exp(λ′/h)}^{ν})) − 𝔼 (f (σ_{T,h;exp(λ′′/h)}^{ν})) → 0,

faster than any power of h. Observe that we are considering times which are of different order of magnitudes, when h is small, and still we are observing a nearly constant 𝔼 (f (σ_{T,h;t}^{ν})). For a second way in which Theorem 2 can be seen

as expressing the presence of a “plateau”, we can think of plotting 𝔼 (f (σ_{T,h;t}^{ν})) versus log(t), rather than versus t. This is somewhat the natural graph to consider, if one is interested in the order of magnitude of the relaxation time. If the log(t)- axis is drawn in the proper scale, amounting to replacing it with h log(t), then, when h is small, Theorem 2 tells us that the graph should be close to that of a step function which jumps at the point λ_{c}, from the value ⟨f⟩_{T,−} to the value ⟨f⟩_{T,+}.

The relation between the constant λ_{c}(T) and some quantities related to the equilibrium Ising model can best be explained by presenting an heuristic reasoning which lies behind Theorems 1 and 2. The heuristics is presented next in the case d = 2. For more on this heuristics including a different way of approaching it and some of its history see [RTMS].

The first ingredient of the heuristics is the idea of looking at an individual droplet of the stable phase (roughly the (+)-phase, since h is small) in a background given by the metastable phase (roughly the (−)-phase). Let S be the shape of that droplet, which a priori can be arbitrary. Say that l^{2} is the volume (i.e., the number of sites) of the droplet, and let us find an expression for the free-energy of such a droplet. This free-energy may be seen as coming from two main contributions. There should be a bulk term, proportional to l^{2}. This term should be obtained by multiplying l^{2} by the difference in free-energy per site between the (+)-phase and the (−)-phase in the presence of a small magnetic field h > 0. This difference in the free-energy per site of the two phases should come only from the term in the Hamiltonian which couples the spins to the external field and should therefore be given by 2m^{∗}h/2 = m^{∗}h. The other relevant contribution to the free-energy of the droplet should come from its surface, where there is an interface between the (+)-phase and the (−)-phase. This contribution is proportional to the length of the interface, which is of the order of l. It should be multiplied by a constant w_{S} which depends on the shape of the droplet. This constant w_{S} represents the excess free-energy per unit of length integrated over the surface of the droplet when its scale is changed so that its volume becomes 1. Adding the pieces, we obtain for the free-energy of the droplet the expression

Φ_{S}(l) = −m^{∗}hl^{2} + w_{S}l.

The two terms in this expression become of the same order of magnitude, in case l is of the order of 1/h. Therefore, it is natural to write l = b/h, with a new variable b ≥ 0. This yields

Φ_{S}(b/h) = \frac{ϕ_{S}(b)}{h},

where

ϕ_{S}(b) = −m^{∗}b^{2} + w_{S}b.

This very simple function takes the value 0 at b = 0, grows with b on the interval [0, B_{c}^{S},], where B_{c}^{S} = B_{c}^{S}(T) = \frac{w_{S}}{2m^{∗}}, reaching its absolute maximum ϕ_{S}(B_{c}^{S}) = \frac{(w_{S})^{2}}{4m^{∗}} = A^{S}(T) = A^{S} at the end of this interval. Then it decreases with b on the semi-infinite interval [B_{c}^{S}, ∞), converging to −∞ as b → ∞.

Metastability is then “understood” from the fact that systems in contact with a heat bath move towards lowering their free-energy, so that the presence of a free-energy barrier which needs to be overcome in order to create a large droplet of the stable phase with any shape keeps the system close to the metastable phase. Subcritical droplets are constantly being created by thermal fluctuations, in the metastable phase, but they tend to shrink, as dictated by the free-energy landscape. On the other hand, once a supercritical droplets is created due to a larger fluctuation, it will grow and drive the system to the stable phase, possibly colliding and coalescing in its growth with other supercritical droplets created elsewhere. As a function of h, the linear size of a critical droplet, B_{c}^{S}/h, blows up as h ↘ 0. One can then, in a somewhat circular, but heuristically-meaningful way, say that the macroscopic free-energy of droplets is indeed a relevant object of consideration. One can also hope then that sharp theorems could be conjectured and possibly proven regarding the asymptotic behavior of quantities of interest in the limit h ↘ 0.

Regarding the shape of the droplet, the height of this barrier is minimized by minimizing the value of the constant w_{S}. It is a fact (see [DKS]) that indeed one can introduce a well defined surface tension function between the (+)-phase and the (−)-phase, and that it produces a single convex shape S which minimizes w_{S}. This shape is called the Wulff shape. We will simplify the notation by omitting the subscript S when it is the Wulff shape. In particular,

B_{c} = \frac{w}{2m^{∗}} , A = \frac{w^{2}}{4m^{∗}} .

Based on the expression above for the free-energy barrier, one predicts the rate of creation of supercritical droplets with center at a given place to be exp (\frac{−βA}{h}).

In what follows now we write d instead of 2, to make the role of the dimension clear in the geometric argument which comes next. We are concerned with an infinite system, and we are observing it through a local function f, which depends, say, on the spins in a finite set Supp(f). For us the system will have relaxed to equilibrium when Supp(f) is covered by a big droplet of the plus-phase, which appeared spontaneously somewhere and then grew, as discussed above. We want to estimate how long we have to wait for the probability of such an event to be large. If we suppose that the radius of supercritical droplets grows with a speed v, then we can see that the region in space-time where a droplet which covers Supp(f) at time t could have appeared is, roughly speaking, a cone with vertex in Supp(f) and which has as base the set of points which have time-coordinate 0 and are at most at distance tv from Supp(f). The volume of such a cone is of the order of (vt)^{d}t. The order of magnitude of the relaxation time, t_{rel}, before which the region Supp(f) is unlikely to have been covered by a large droplet and after which the region Supp(f) is likely to have been covered by such an object can now be obtained by solving the equation

(vt_{rel})^{d} t_{rel} exp (−\frac{βA}{h}) = 1.

This gives us

t_{rel} = v^{−d/(d+1)} exp (\frac{βA}{(d + 1) h}) .

In order to use this relation to predict the way in which the relaxation time scales with h, one needs to figure out the way in which v scales with h. If we suppose, for instance, that v does not scale with h, or at least that if it goes to 0, as h ↘ 0, it does it so slowly that

(1)  lim_{h↘0} h^{d−1} log v = 0,

then we can predict that

t_{rel} ≃ exp (\frac{βA}{(d + 1) h}) = exp (\frac{λ_{c}}{h} ) ,

where

(2)  λ_{c} = \frac{βA}{d + 1} = \frac{βw^{2}}{(d + 1) 4m^{∗}} = \frac{βw^{2}}{12 m^{∗}} .

The heuristics above may seem extremely crude. Potentially the interaction between droplets could spoil the whole picture and lead to a much faster decay. In the opposite direction, even if the droplet picture makes sense, their speed of growth could be so slow that (1) could fail an therefore the relaxation time would be much larger than predicted above.

One of the major contributions of [SS] is to prove that indeed λ_{c} in Theorem 2 is given by (2). This means that close to the phase transition region the time evolution can be well described in first approximation by the highly simplified droplet dynamics.

Acknowledgements: It is a pleasure to thank Senya Shlosman for the collaboration in [SS] and other related projects. This work was supported by the N.S.F. grant DMS-9703814.

References

[AL] M. Aizenman and J. L. Lebowitz, Metastability effects in bootstrap percolation, J. Phys. A 21 (1988), 3801–3813.

[BC] G. Ben Arous and R. Cerf, Metastability of the three dimensional Ising model on a torus at very low temperatures, Electronic Journal of Probability 1 (1996), 1–55.

[BM] K. Binder and H. Müller-Krumbhaar, Investigation of metastable states and nucleation in the kinetic Ising model, Physical Review B 9 (1974), 2328–2353.

[CL] E. N. M. Cirillo and J. Lebowitz, Metastability in the two-dimensional Ising model with free boundary conditions, J. of Stat. Phys. 90 (1998), 211–226.

[CO] E. N. M. Cirillo and E. Olivieri, Metastability and nucleation for the Blume-Capel model – different mechanisms of transition, J. of Stat. Phys. 83 (1996), 473–554.

[DS] P. Dehghanpour and R. H. Schonmann, Metropolis dynamics relaxation via nucleation and growth, Commun. Math. Phys. 188 (1997), 89–119.

[DKS] R. L. Dobrushin, R. Kotecký and S. B. Shlosman, Wulff construction: a global shape from local interaction, AMS translations series, Providence (Rhode Island), 1992.

[Isa] S. N. Isakov, Nonanalytic features of the first order phase transition in the Ising model, Commun. Math. Phys. 95 (1984), 427–443.

[Nev] E. J. Neves, A discrete variational problem related to Ising droplets at low temperatures., J. of Stat. Phys. 80 (1995), 103–123.

[RTMS]P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides, Metastable lifetimes in a kinetic Ising model: dependence on field and system size, Phys. Rev. E 49 (1994), 5080–5090.

[Sch1] R. H. Schonmann, Slow droplet-driven relaxation of stochastic Ising models in the vicinity of the phase coexistence region, Commun. Math. Phys. 161 (1994), 1–49.

[Sch2] R. H. Schonmann, Theorems and conjectures on the droplet driven relaxation of stochastic Ising models, in Probability theory of spatial disorder and phase transition, G. Grimmett, ed., Kluwer Publ. Co (1994), 265–301.

[SS] R. H. Schonmann and S. B. Shlosman, Wulff droplets and the metastable relaxation of kinetic Ising models, Commun. Math. Phys. (to appear).

Roberto H. Schonmann Mathematics Department UCLA Los Angeles CA 90095 U.S.A.

182

Documenta Mathematica · Extra Volume ICM 1998 · III ·