# Recent Results on the Decay of Metastable Phases

###### Abstract

We review some aspects of current knowledge regarding the decay of metastable phases in many-particle systems. In particular we emphasize recent theoretical and computational developments and numerical results regarding homogeneous nucleation and growth in kinetic Ising and lattice-gas models. An introductory discussion of the droplet theory of homogeneous nucleation is followed by a discussion of Monte Carlo and transfer-matrix methods commonly used for numerical study of metastable decay, including some new algorithms. Next we discuss specific classes of systems. These include a brief discussion of recent progress for fluids, and more exhaustive considerations of ferromagnetic Ising models (i.e., attractive lattice-gas models) with weak long-range interactions and with short-range interactions. Whereas weak-long-range-force (WLRF) models have infinitely long-lived metastable phases in the infinite-range limit, metastable phases in short-range-force (SRF) models eventually decay, albeit extremely slowly. Recent results on the finite-size scaling of metastable lifetimes in SRF models are reviewed, and it is pointed out that such effects may be experimentally observable.

## 1 Introduction

Metastable phases are common in nature, and some of them have such extremely long lifetimes that they are practically indistinguishable from equilibrium phases. For example, even though the stable phase of carbon at room temperature and atmospheric pressure is graphite, the claim that diamonds are forever is not likely to be challenged on physical grounds. Other examples of metastable phases with average lifetimes that are measured in milliseconds to days, rather than in millions or billions of years, are supercooled or supersaturated fluids [1–32], ferroelectrics [33–38], the small single-domain ferromagnetic particles important in paleomagnetism and high-density recording media [39–49], and vortex states in superconductors [50]. Further possible examples are the supercooled quark/gluon plasma associated with the QCD confinement transition [51, 52] and the “false vacuum” associated with the electroweak transition [53–57], both of which may have played important roles in the early development of the universe.

To clarify what we understand by a metastable phase, it is hard to
improve on the empirical descriptions due to
Penrose and Lebowitz [58, 59] and Sewell [60].

(i) The free energy of the system is not fully minimized.

(ii) Only one thermodynamic phase is present, and for sufficiently
small and slow perturbations the usual laws of reversible thermodynamics
apply.

(iii) A system starting in the metastable phase
is likely to take a long time to escape.

(iv) Escape is irreversible: once out of the metastable phase,
the system is extremely unlikely to return.

We concentrate on systems for which the order parameter is a nonconserved scalar (Model A in the Hohenberg-Halperin scheme of dynamic universality classes [61]) so that the metastability is imposed by an applied external field. (The closely related phenomenon of hysteresis imposed by an oscillating field is discussed by Acharyya and Chakrabarti elsewhere in this volume [62].) Even so, many of the phenomena we discuss can be generalized to systems with a multidimensional and/or conserved order parameter .

The wide variety of contexts in which metastability has been studied, often independently, makes it difficult to present a discussion that is equally accessible to all potentially interested readers. In this review we have chosen to use the language of Ising models. The most important reason for this choice is our belief that the present theoretical understanding of metastable phases and their modes of decay is most highly developed in applications to the relatively simple kinetic Ising models and their equivalent formulations as lattice-gas models. A secondary reason is the high symmetry of the Ising formulation. Throughout this paper we therefore use Ising ferromagnets as generic examples. Below we define these models and give the standard transformations that will enable the reader to convert our results to the lattice-gas language appropriate to, e.g., fluids, binary mixtures, and adsorbate systems.

An Ising ferromagnet is defined by a Hamiltonian (i.e., an energy functional),

(1) |

where (“up” or “down”) is a binary variable, or “spin,” at site , is the applied field, and the sums run over all sites on a -dimensional lattice. The interaction energies are positive and symmetric under interchange of the site indices and , and without loss of generality we set =0. By requiring that is independent of we ensure the spatial homogeneity of the energy functional. We also note that is invariant under the transformation . The order parameter conjugate to is the “magnetization” or “polarization”,

(2) |

As is well known [63], this model is equivalent to a two-state lattice-gas model with local concentration variables =0 or 1 (“empty” or “occupied”), for which Eq. (1) takes the form

(3) |

where the quantities appearing in the two equivalent formulations of the Hamiltonian, Eq. (1) and Eq. (3), are linked by the transformations

(4) |

Here are attractive lattice-gas interaction energies and is the chemical potential, whose value at coexistence (i.e., for =0) is =. The chemical potential is related to the (osmotic) pressure as , where is Boltzmann’s constant times the absolute temperature, and is the pressure at coexistence. The order parameter conjugate to is the density,

(5) |

The energy functional is not a quantum-mechanical Hamiltonian and does not impose a unique dynamic on the system. When one uses an Ising or lattice-gas model to describe the kinetics of a physical system, one must therefore define a specific stochastic dynamic to simulate the interactions with the system’s surroundings. Usually it is physically most realistic to choose a dynamic which is local in the sense that it only allows transitions involving a single site or a pair of nearest-neighbor sites, and in this review we mostly limit our attention to such dynamics. Simple examples are the Metropolis [64] and Glauber [65] dynamics, which are among the algorithms we consider in the context of Monte Carlo simulations in Sec. 3.1. (Under somewhat restrictive conditions, the Glauber dynamic has been obtained from the quantum mechanical equations of motion for a system of distinguishable spin-1/2 particles weakly coupled to an infinitely large quantum reservoir [66].)

The Ising (lattice-gas) model below its critical temperature and in the absence of an applied field (with =) has two coexisting, ordered phases of equal free energy: one with positive magnetization (high density) and one with negative magnetization (low density). This degeneracy is lifted by applying a nonzero field (changing away from ): the phase whose magnetization is parallel to the field (for which has the same sign as ) becomes the unique equilibrium phase, and the one with the opposite magnetization (“wrong” density) becomes metastable.

The system can be prepared in the metastable phase by equilibrating it in a nonzero field that is then instantaneously reversed. Although the system is no longer in equilibrium immediately after the field reversal, it is nevertheless stable against small fluctuations, and its thermodynamic properties are similar to those it would have in the equilibrium phase. This is because configurations obtained by flipping small clusters of neighboring spins cost more free energy by introducing an interface between previously parallel spins than they gain by lowering the field energy. In order for this metastable phase to decay, it is therefore necessary that a cluster is created that is sufficiently large for the free energy gained by aligning more spins with the field to just outweigh the cost of breaking the necessary extra bonds. Such a fluctuation corresponds to a local maximum in the free-energy landscape and is usually called a “critical nucleus” or a “critical droplet”. Once randomly created through a thermal fluctuation in the metastable phase, it is likely to grow further. During this growth period the now “supercritical” droplet incorporates spins from the metastable phase into a growing, and soon macroscopic, domain of the equilibrium phase. The timescale for the creation of a critical droplet is in general much longer than that characteristic of the subsequent growth.

Real systems are of course more complicated than the simple picture discussed above. Nevertheless, it covers the essential physics of nucleation in a variety of situations. Some of the additional complications that arise in the study of metastability in fluids with continuous degrees of freedom are briefly discussed in Sec. 4.

The nucleation mechanism described above, which is usually known as homogeneous or thermally activated nucleation, will be the focus of our attention. It is the dominant mechanism in systems that are free of defects, or in which the defect concentration is low, or in which the defects are much smaller than the critical droplet size. Whenever this is not the case, heterogeneous nucleation on defects or interfaces may be the dominant process for producing droplets of the equilibrium phase [67]. However, the subsequent growth process depends little on the nucleation mechanism [33, 68].

The organization of the remainder of this review is as follows. Section 2 contains a brief review of classical nucleation theory and some of the general results obtained in the “post-classical” period after the 1940’s. In Sec. 3 we review numerical Monte Carlo (in Sec. 3.1) and transfer-matrix (in Sec. 3.2) methods to study nucleation and metastable decay, primarily in kinetic Ising models. Following these general sections we consider specific classes of systems: fluids in Sec. 4, Ising models with weak long-range forces (WLRF models) and their relation to mean-field approximations in Sec. 5, and Ising models with short-range forces (SRF models) in Sec. 6. The subsections in Sec. 5 are Sec. 5.1, which contains specific nucleation rates for the ramified critical droplets characteristic of WLRF models, and Sec. 5.2, which contains numerical results. The subsections in Sec. 6 are Sec. 6.1, which contains specific nucleation rates for the compact critical droplets characteristic of SRF models, Sec. 6.2, in which we discuss the interplay between nucleation and the growth of supercritical droplets, Sec. 6.3, in which we consider finite-size effects, and Sec. 6.4, in which we consider effects of the discrete lattice at low temperatures. Finally, in Sec. 7, we give a concluding summary and discussion.

## 2 Droplet Theory of Homogeneous Nucleation

The first scientific description of metastability may well be Fahrenheit’s experiments with supercooled water, published in 1724 [1], and further observations were reported during the eighteenth and early nineteenth centuries by several workers [19]. The basis for a theoretical understanding was laid 150 years after Fahrenheit’s paper with van der Waals’ [2] and Maxwell’s [3] early mean-field approaches and Gibbs’ realization that the reversible work of formation of a droplet of the stable phase in a metastable background is a “measure of the stability” of the metastable phase [4–6]. Although Gibbs’ work clearly states the basic equations of what is today known as “Classical Nucleation Theory” (CNT), further development did not occur until the period from the 1920’s through the 1940’s. Gibbs’ essentially thermodynamic (i.e., static) approach was then complemented by kinetic considerations by Volmer and Weber [7], Farkas [8], Becker and Döring [9], and Zeldovich [10, 11], whereas Bijl [12], Frenkel [13], and Band [14] focused on calculating a constrained, metastable partition function through a cluster-expansion procedure inspired by that of Mayer [69]. The basic result for the nucleation rate can be written in a Van’t Hoff-Arrhenius [70, 71] form,

(6) |

where is a non-exponential prefactor, is the free-energy cost of a critical droplet (Gibbs’ “measure of stability”), and = is the inverse temperature. (A compact historical sketch is given by Dunning [19]. For more recent general reviews of metastability and the kinetics of first-order phase transitions see, e.g., Refs. [31, 72–74].)

The “post-classical” developments in nucleation theory are mainly efforts to evaluate and/or the prefactor in Eq. (6) for specific systems. This observation in itself provides an important insight: the lifetime of a metastable state depends crucially on the structure and dynamics of the excitations through which it decays [75].

The relevance of the classification scheme for dynamic critical phenomena due to Hohenberg and Halperin [61] to metastable decay is not clear. In addition to those aspects of the physics that influence both the static and the dynamic universality, such as the interaction range and the symmetries and dimensions of both the system itself and the order parameter, the conservation laws that influence the dynamic scaling should also be important for the nucleation rate.

The “metastable thermodynamics” that can be observed during the period before the decay takes place is almost indistinguishable from true equilibrium thermodynamics. This has inspired efforts to treat metastability through suitable generalizations of equilibrium statistical mechanics. Notable among these are three papers by Langer [76–78], in which he shows that for a wide class of models, whose dynamics can be described by a Fokker-Planck equation, the homogeneous nucleation rate (per unit time and volume) may be written as

(7) |

where is the imaginary part of a “metastable” free-energy density , which is obtained by analytically continuing the equilibrium free-energy density into the metastable phase. The “kinetic prefactor” contains all dependence on the specific dynamic. This important work brings together aspects of CNT with results on droplet theory and analytic continuation from Andreev [79] and Fisher [80] and on thermally activated processes from Landauer and Swanson [81] and Kramers [82].

Despite its apparent simplicity and its similarity to the corresponding relation in quantum mechanics (obtained by replacing by the energy and by ) no general proof of Eq. (7) exists, and its domain of validity remains unclear. Analytic verification for specific models has been reported, e.g., by Newman and Schulman [83] for a Curie-Weiss ferromagnet, by Roepstorff and Schulman [84] for an urn model, by Gaveau and Schulman [85] for a class of models including the droplet model of condensation [11, 79, 80], and by Penrose [86] for the droplet model with Becker-Döring dynamics. A different approach is to assume the validity of Eq. (7) and calculate analytically or numerically. Again, analytic calculation requires a specific model for the fluctuations included in the calculation of the analytic continuation . Such field-theoretical calculations were done by Coleman and Callan [53, 54] for the “false vacuum” in quantum-field theory, by Büttiker and Landauer [87, 88] for a one-dimensional (1D) overdamped sine-Gordon chain, by McCraw [89] for a 1D Kac model with algebraically decaying interactions, by Günther, Nicole, and Wallace [90], who generalized Langer’s field-theoretical calculation to arbitrary spatial dimension (see Sec. 6.1), by Zwerger [91] for a field model, by Cottingham et al. [32] for a model of bubble formation in fluids, by Braun [43, 44] for a model of switching in single-domain ferromagnetic particles, and by Klein and Unger [92, 93] and Gorman et al. [94, 95], who used field theories to study WLRF systems near the classical spinodal (see Sec. 5).

As may be understood from the heuristic discussion of the decay of a metastable phase presented in Sec. 1, the metastable lifetime depends not only on the rate of nucleation of critical droplets, but equally importantly on the subsequent rate of growth of the supercritical droplets and on the interplay between these processes of nucleation and growth. This was realized by Kolmogorov [96], Johnson and Mehl [97], and Avrami [98] (KJMA) at about the same time as CNT was being developed, and a few years later by Evans [68].

## 3 Computational Methods

In contrast to analytical calculations of or , nonperturbative numerical methods do not require a priori knowledge about the critical excitations. In this section we briefly discuss two classes of such methods: some varieties of the well-known Monte Carlo (MC) method for numerical simulation, and a recently introduced extension of the equilibrium transfer-matrix (TM) method to also encompass metastable phases.

### 3.1 Monte Carlo Methods

Simulations using standard Metropolis or heat-bath dynamics with spatially local updates (see, e.g., Ref. [99]) have remained the computational methods of choice for studying metastable decay in systems with nonconserved order parameter. The spatial locality of the algorithm preserves the free-energy barriers that dominate the nucleation rate, leading to relatively faithful representations of metastable dynamics in physical systems. (This is not generally true for cluster algorithms such as the one due to Swendsen and Wang [100], which are used to accelerate equilibrium simulations near criticality. However, see further discussion below.)

In equilibrium studies, the choice between the Metropolis [64] and heat-bath algorithms (the latter also known as the Gibbs sampler or, when applied to the Ising model, the Glauber [65] dynamic) is largely a matter of convenience, and it is often considered so trivial that it is not even stated explicitly. This habit sometimes has carried over to dynamical studies where, in our opinion, attention should be given to finding the dynamic best representing the physical system to be simulated. For completeness we give below the probabilities for an allowed transition from state to for these two algorithms [101–103].

In the Metropolis algorithm a candidate state is selected according to a proposal distribution, which vanishes for forbidden transitions. (The decision about which transitions should be allowed is part of the full definition of the algorithm. For an Ising system one could for instance allow single-spin flips, nonlocal cluster flips as in the Swendsen-Wang [100] and related [104–106] algorithms, or nearest-neighbor spin exchanges as in the conserved-order-parameter Kawasaki dynamic [107].) The proposed transition is accepted with probability one if it leads to a reduction in energy, whereas an increase in energy is accepted with probability given by a Boltzmann factor, so that the acceptance probability may be written as

(8) |

In contrast, the heat-bath algorithms accept any state to which a transition from is allowed, with the equilibrium probability over the set of all accessible states :

(9) |

Both the Metropolis and the heat-bath algorithm satisfy detailed balance, and with ergodicity they eventually converge to thermodynamic equilibrium. However, the detailed Markov processes generated are not identical. Beside the choice of the Metropolis or heat-bath transition probability, more subtle differences between algorithms may also influence the results of simulations, and even though data obtained by different local algorithms can often be connected by simply rescaling the time, the rescaling factor may depend nontrivially on temperature and distance from the coexistence curve. For example, in a recent study of metastable decay in the two-dimensional Ising model with the Metropolis dynamic [108], Rikvold et al. found that the manner in which the candidate site for the next spin update was chosen (sequentially or randomly) affected the observed field dependence of the kinetic prefactor in Eq. (7). Similarly, in a recent study of magnetization relaxation in the three-dimensional Ising model at the critical point, Ito demonstrated that the detailed nature of the finite-size effects depends not only on the choice between the Metropolis and heat-bath dynamics, but also on whether a sequential or a checker-board update scheme was used [109]. In choosing the MC algorithm for a particular study, one should therefore consider whether the quantities of interest are universal in the sense that they are the same for all local algorithms, or whether they depend on the dynamical details of the algorithm used.

For low temperatures and close to coexistence, both the local algorithms discussed above and the physical dynamic they simulate spend most of the time creating short-lived, microscopic excitations in the metastable phase. Because of the smallness of the subcritical clusters, this is also true for the modified Swendsen-Wang cluster dynamic used in MC simulations by Ray, Tamayo, and Wang [104, 105] and studied theoretically by Martinelli, Olivieri, and Scoppola [106]. The brute-force way around this problem is to apply more computer power in the form of various kinds of supercomputers or special-purpose machines (as, e.g., in Refs. [108, 110, 111]), but inevitably the slowness of the dynamic makes this approach impracticable. Very recently, two novel MC algorithms have been introduced, which make simulations deep in the metastable region feasible. Both methods utilize absorbing Markov chains [112]. One of them, introduced by Novotny [113–115], generalizes the rejection-free -fold way algorithm of Bortz, Kalos, and Lebowitz [116] and achieves CPU-time savings of many orders of magnitude relative to the local algorithms. This speedup is obtained without changing the underlying dynamic in any way. The other method, introduced by Lee et al. [117], combines absorbing Markov chains with the Multicanonical method [118]. The resulting dynamic depends on the microscopic configurations only through their projection onto the order parameter. Nevertheless, most qualitative and quantitative features of the dependences of the metastable lifetimes on field, temperature, and system size agree with theoretical results and direct simulations for local dynamics. By comparison with such more traditional methods, this algorithm may help deepen our understanding of universality in metastable decay. Preliminary results from both methods are promising (some are presented in Secs. 6.3 and 6.4), and we hope to review further progress in the future.

It is common in dynamical MC simulations of metastable decay (regardless of the particular algorithm employed) to study the relaxation of the order parameter, starting from the metastable phase. This approach is closely related to the use of nonequilibrium relaxation functions, introduced by Binder [119]. It has been used for SRF models in two [37, 104, 108, 110, 111, 120–129], three [38, 105, 130–132], and higher [133] dimensions, and also for WLRF models [129, 131, 134, 135]. Some recent results are presented in Sec. 6.

### 3.2 The Constrained-Transfer-Matrix Method

The Constrained-Transfer-Matrix (CTM) method introduced by one of us [136] is particularly suited for numerical calculation of the analytically continued free-energy density in the field-theoretical droplet theory discussed in Sec. 2, without needing a theoretical model of the critical droplet geometry. The method extends the usual concept of the transfer matrix (TM) [137] to also include constrained nonequilibrium states. Here we give a brief description of the technique. More detailed discussions can be found in Refs. [94, 95, 136, 138, 139].

In a standard TM calculation, an lattice is considered in the limit , and the Hamiltonian (the energy functional) is written as a sum of layer Hamiltonians, , where depends only on the configurations and of two adjacent layers. The TM is a positive matrix,

(10) |

where the second equality represents a standard eigenvalue expansion. By the Perron-Frobenius theorem [137], the dominant eigenvalue is positive and nondegenerate, and the corresponding eigenvector is the only one with elements that can all be chosen positive. From one can calculate by standard methods the probability densities, correlation functions, and partition function that fully describe the equilibrium phase [137].

The CTM method generalizes this well-known equilibrium technique by associating with each eigenvalue a “constrained” TM so that “constrained” joint and marginal probability densities are defined in analogy with the equilibrium (=0) case:

(11) |

It was pointed out by McCraw, Schulman, and Privman [127, 140, 141] that the constrained marginal probability densities , as defined above, can be interpreted as actual probability densities over single-layer configurations in a constrained phase. To obtain explicitly the constrained joint probability densities that define the inter-layer correlations in the constrained phase, one must determine the constrained TM, . It is chosen to commute with , and in order to ensure convergence of towards stochastic independence as , each eigenvalue is reweighted to become , so that the dominant eigenvalue of is . A “constrained” free-energy density is defined by

(12) |

where Ln() is the principal branch of the complex logarithm. For =0 this reduces to the standard equilibrium result, . For 0 the eigenvector , corresponding to the dominant eigenvalue of , is orthogonal to and therefore cannot have all positive elements. As a result, cannot in general be a positive matrix, and the second term in Eq. (12) becomes complex-valued. It has been observed by Günther et al. [139] that this second term can be considered as a complex generalization of the Kullback discrimination function (see, e.g., Ref. [142]) for with respect to the divergent ‘probability density’ obtained by substituting for in Eq. (3.2).

## 4 Fluid Systems

As mentioned in Sec. 1, metastable decay in fluid systems, such as condensation of a supersaturated vapor or freezing of a supercooled liquid, presents complications not encountered in Ising or lattice-gas systems. In this section we mention some of these difficulties and give a few references.

Most of the early theoretical work on metastability that led to the formulation of CNT [4–14] was explicitly concerned with fluids, as were the early mean-field approaches of van der Waals [2] and Maxwell [3]. Extensive reviews can be found in Refs. [17–19].

One of the most obvious differences between fluids and lattice-gas systems is that fluids have continuous degrees of freedom associated with droplet translation and rotation. Their effects were considered by Lothe, Pound, and collaborators [20, 21], who evaluated the prefactor in Eq. (6) and obtained an increase in the predicted condensation rate, relative to earlier theories, on the order of 10. (See, however, a recent discussion in Ref. [25].) The long controversy over the Lothe-Pound result emphasizes the importance and difficulty of understanding the pre-exponential factors in the nucleation rate.

An important conceptual problem arises when one tries to formulate a precise droplet definition [25–29]. It was already considered by Gibbs, who introduced a “dividing surface” between the liquid droplet and the vapor, but it does not seem yet to have reached a fully satisfactory solution for fluid systems [145]. Even in the much simpler Ising model, the proper definition of droplets has been realized only relatively recently in terms of the “Fortuin-Kastelyn-Coniglio-Klein-Swendsen-Wang” cluster definition. For general and historical discussions of the Ising droplet definition, see Refs. [146, 147], and for discussions in the context of metastability, see Refs. [104–106, 145].

A serious problem of a numerical nature is the difficulty of accurately locating the coexistence curve in a model in which its position is not given by symmetry, as it is for the Ising or binary lattice-gas model. (Even in more complicated lattice-gas models this is a nontrivial problem [148].) The extremely strong dependence of the nucleation rate on the distance from the coexistence curve will cause even a small error in its location to produce large errors in numerical estimates of the droplet free energy and the pre-exponential factor.

## 5 Ising Models with Weak Long-Range Forces

In the limit of infinite interaction range, ferromagnetic Ising models with weak, long-range interactions (WLRF models) have been rigorously shown to possess infinitely long-lived metastable phases [58, 59, 149]. These phases lie on a hysteresis (or van der Waals) loop that reaches beyond the Maxwell construction in the phase diagram, and whose critical points correspond to sharp spinodals beyond which the metastable phases disappear. In this respect a WLRF model is similar to a mean-field theory [2, 3, 39], except that in a mean-field theory the only allowed form of fluctuation is a spatially uniform change of the order parameter. The free-energy cost of such a change is extensive in the system volume [89, 135], and the mean-field approximation therefore predicts that metastable lifetimes are infinite in the thermodynamic limit.

Since the behavior of a number of physical systems, including superconductors and long-chain polymer mixtures, are often well described by WLRF models, there has been an interest in studying metastability in these models. Analytic approaches include the early field-theory treatment of a model fluid with conserved order parameter by Cahn and Hilliard [15, 16], the Fokker-Planck equation [89, 135, 150], renormalization-group analyses [151, 152], and arguments from random long-range bond percolation [153, 154], in addition to analytic continuation of the free-energy density [83, 92–95]. Numerical approaches include MC simulations [129, 131, 134, 135, 153–155], traditional TM calculations [156], and calculations of by the CTM method described in Sec. 3.2 [94, 95, 136, 143].

As pointed out by Klein and coworkers [92, 93, 153–155], the critical droplet in a WLRF model near the spinodal is a ramified structure whose local order parameter is only slightly different from that of the metastable phase. (Note the contrast with the SRF case, in which the critical droplet is compact. See Sec. 6.1.) Only during the supercritical growth phase does the droplet center compactify, leading to a measurable change in the global order parameter. As a consequence, it is difficult to determine the nucleation time accurately from MC simulations by monitoring the order parameter or, equivalently, the nonequilibrium relaxation function. Therefore, several new techniques have been developed to analyze MC data for metastable decay in WLRF systems, including the analysis of recrossing events [129], monitoring the number of spins in the largest cluster [155], and intervention techniques [154]. Some of the recent work in this area has been reviewed in Ref. [145].

The difficulties associated with accurately determining the nucleation rate in WLRF models make them particularly attractive candidates for study by alternative numerical techniques, such as the CTM method. A simple WLRF system suitable for CTM investigation is the Quasi-One-Dimensional Ising (Q1DI) model, introduced by Novotny et al. [156]. In its simplest form it is a chain of layers, each of which contains Ising spins, =1. Each spin interacts ferromagnetically with each of the 2 spins in the adjacent layers with coupling 0 (these are the nonzero in Eq. (1) for this model), and with an external field . In terms of the single-layer magnetizations, =, the Hamiltonian is

(13) |

### 5.1 Nucleation

In a WLRF system like the Q1DI model, where droplet “interfaces” are difficult to define, the onset of nucleation is detected by a rapid increase in the magnetization of some layer to within a small neighborhood of the equilibrium value. The physical picture is that of a bell-shaped magnetization profile that quickly develops into a radially expanding front separating the two competing phases. It was shown by Gorman et al. [94] that the free-energy cost of nucleation for the Q1DI model can be approximated by mapping the free-energy density functional of the model to a Ginzburg-Landau form and solving the resulting Euler-Lagrange equation. A expansion of the potential around the exactly known [83, 135] mean-field spinodal field , in which , gives the following free-energy cost for a -dimensional WLRF system with force range [92, 93]:

(14) |

where is a nonuniversal function of . For the Q1DI model = and =1 [152, 156]. By solving the Euler-Lagrange equation numerically, one can also take into account the corrections to the expansion, as well as discrete lattice effects [94].

The prefactor to the Van’t Hoff-Arrhenius term in can be calculated by expanding the free-energy density functional about the stationary points corresponding to the metastable phase and the critical droplet. The determinants of the two resulting Schrödinger operators give [94]

(15) |

where is a nonuniversal function of , is the system volume, and is the volume of the subspace in which the droplet itself is free to move without a cost in free energy. For the Q1DI model, == [94].

From Eqs. (14) and (15) we see that for large , unless is extremely close to , the free-energy cost of surmounting the nucleation barrier is large, so the exponential factor sets the scale for the metastable lifetime. However, for small , or for , the lifetime is more strongly dependent on the particulars of the dynamic and on the detailed structure of the saddle point in the free-energy functional. The finite-range-scaling of near the spinodal is found by recasting Eqs. (14) and (15) in terms of a scaling variable =, giving

(16) |

with a scaling function dependent only on and . The dynamic prefactor , taken from the lowest eigenvalue of the Schrödinger operator corresponding to the saddle point, scales as , but is independent of . The finite-range scaling of the nucleation rate is thus

(17) |

with a scaling function dependent only on and .

### 5.2 Numerical Transfer-Matrix Results

The CTM method outlined in Sec. 3.2 was applied by Gorman et al. [94] to Q1DI systems of infinite length and finite cross section , with up to between 100 and 500. Figure 1 shows typical spectra of Re and plotted against at a fixed temperature. These spectra are compared with the analytically continued free energy in the limit , which corresponds to a Curie-Weiss mean-field result (thick, dashed curves). In each calculation, as demonstrated in the figure, a unique was found for which Re computed from Eq. (12) closely approximated the real part of the mean-field metastable free energy. Invariably, the same produced the smallest nonzero value for . As was changed away from towards =0, the value of the metastable was exponentially suppressed.

If the metastable represents the analytically continued free-energy density, then one expects to be a measure of the height of the nucleation barrier, as seen from Eq. (15). Extrapolated CTM estimates of the barrier height are compared in Fig. 2 with the free-energy cost of the saddle-point solution to the Euler-Lagrange equation, which was obtained by numerical integration. The agreement is remarkable and extends over an extremely wide range of fields and temperatures, for most of which the lifetimes are too long for standard Monte Carlo techniques to be useful. Similar results were found in a study of a three-state WLRF model [95], even where two competing metastable states were present, suggesting that Langer’s formula has a wider applicability than previous studies have claimed [85]. These CTM results, as well as the MC results in Ref. [95], are consistent with Ostwald’s empirical “Law of Stages” (Gesetz der Umwandlungsstufen) [19, 157], whereby a metastable phase may decay via intermediate metastable states before reaching equilibrium.

## 6 Ising Models with Short-Range Forces

In contrast to the situation for WLRF models, metastable phases in systems with short-range forces (SRF models) eventually decay, even though their lifetimes may be many orders of magnitude larger than other characteristic timescales of the system [17]. Thus, whereas mean-field theory provides a qualitatively acceptable description of metastability for WLRF systems in the long-range limit, it gives a quite misleading picture for SRF systems. The following subsections are devoted to a discussion of metastability and metastable decay appropriate for SRF systems.

As a prototype for the metastable dynamics of SRF systems, we consider in detail the decay of the magnetization in an impurity-free kinetic nearest-neighbor Ising ferromagnet in an unfavorable applied field. The critical point of this model is in the same static universality class as the liquid/vapor phase transition [158]. (However, the liquid/vapor transition is in a different dynamic universality class: that of Model H [61].) The Hamiltonian is obtained from Eq. (1) by setting = if and are nearest-neighbor sites and =0 otherwise, explicitly yielding

(18) |

where and run over all nearest-neighbor pairs and over all sites on a -dimensional hypercubic lattice of volume , respectively [159].

In a typical MC study of metastable decay in this model one starts from the metastable phase and follows the relaxation of the magnetization, which is closely related to Binder’s nonequilibrium relaxation functions [119] and is directly obtained as an average over the droplet size distribution [80, 123, 132]. The volume fraction of stable phase at time is and are the bulk equilibrium and metastable magnetizations, respectively. The metastable lifetime is typically estimated as the mean-first-passage-time for to a preset value. Monte Carlo studies using this or similar methodology have been performed in two [37, 104, 108, 110, 111, 120–129], three [38, 105, 130–132], and higher [133] dimensions. Several of these studies were analysed in terms of droplet theory, establishing general agreement between theory and simulations. A potentially more accurate, but also more computationally intensive, method to estimate the lifetimes could use the recrossing-event distribution discussed by Paul and Heermann [129] to determine a field-dependent cutoff value for , as suggested by Rikvold et al. [108]. , where

Metastability in the two-dimensional nearest-neighbor Ising ferromagnet has also been studied by traditional TM methods by McCraw, Schulman, and Privman [127, 140, 141], and by the CTM method by Günther, Rikvold, and Novotny [138, 139, 144].

### 6.1 Nucleation

To obtain quantitative comparisons between numerical results and the droplet-based nucleation theory, one must calculate explicitly the field-theoretical expression for the nucleation rate, Eq. (7). This requires as input the free energy of a compact critical droplet , as obtained by CNT [76, 90]. (Note the contrast with the WLRF case, in which the critical droplet is ramified. See Sec. 5.) Sufficiently far below the critical temperature we can calculate this by a standard droplet-theory argument (see, e.g., [18, 72, 73]) that is modified to consider the nonspherical droplets that appear at low because of the anisotropy of the surface tension [108, 138, 139, 144, 160–163]. The free energy of a -dimensional droplet of radius (defined as half the extent of the droplet along a primitive lattice vector) and volume is

(19) |

The quantity is a temperature- and, in principle, field-dependent proportionality factor which relates the surface contribution to with the droplet volume [162, 164]. The difference in bulk free-energy density between the metastable and stable states is , which takes some account of droplet nesting through the magnetization difference [160, 165]. Maximizing yields the critical radius,

(20) |

where is the surface tension along a primitive lattice vector [162], and the free-energy cost of a critical droplet,

(21) |

In addition to , the critical droplet is characterized by other degrees of freedom, including the critical growth mode, droplet translations, and deformations represented by capillary waves on the droplet surface. The field-theoretical expression for the nucleation rate, Eq. (7), properly accounts for the effects of these additional degrees of freedom, and it is obtained explicitly by a saddle-point calculation that yields [76–78, 90]

(22) |

with

(23) |

The quantity is a nonuniversal function of the temperature only,

(24) |

is a universal exponent related to excitations on the surface of the critical droplet [76, 90], and gives the dependence of the kinetic prefactor [77, 78]. The kinetic prefactor is the only part of that may depend explicitly on the specific dynamic.

The Becker-Döring cluster dynamics is defined in terms of a master equation for the probability distribution of -particle clusters. The original theory only allows to change by 1 [11], but later modifications also allow cluster coagulation and fragmentation [122]. The CNT result for the exponential part of can be obtained from this approach. The existence and uniqueness of a solution that is metastable in the sense of the Penrose-Lebowitz-Sewell criteria have been proven both for the original dynamic [166, 167] and for the coagulation-fragmentation generalization [168]. However, describing the clusters only in terms of their particle numbers does not suffice to obtain the pre-exponential factor in Eq. (6) for a particular system. This was demonstrated explicitly by Binder and Stauffer [125], who obtained a result formally analogous to Eq. (22) by using a droplet model in which the individual clusters were characterized by several coordinates in addition to the particle number .

For =2, there is substantial numerical evidence that =1, as predicted by Eq. (24). This is obtained from calculations that do not involve the dynamics, such as analyses of series expansions [160, 169, 170] and transfer-matrix calculations [138, 139, 144]. These studies, as well as MC work [108, 171–173], also indicate that the free-energy cost of the critical droplet is given by Eq. (21) with the zero-field equilibrium values for and . We therefore adopt the notations =, =, and = to emphasize the lack of field dependence in these quantities. The quantity can be obtained with arbitrary numerical precision by combining a Wulff construction with the exact, anisotropic zero-field surface tension [162], and is obtained from the exact Onsager-Yang equation [174]. These general results, that the surface free energy and bulk magnetization of compact critical droplets are determined by the zero-field equilibrium surface tension and magnetization, respectively, are also supported by MC studies of nucleation rates in three dimensions [105, 130, 132].

For dynamics that can be described by a Fokker-Planck equation, it is expected that the kinetic prefactor is proportional to [77, 78, 90], which by Eq. (20) yields =2. This value has been confirmed numerically for the Metropolis algorithm with updates at randomly chosen sites, but it does not appear to apply if the sites are chosen sequentially [108].

The numerical results cited above indicate that CNT, modified by the post-classical results for the prefactor exponents, gives very good agreement with the observed behavior of kinetic Ising and lattice-gas models, even moderately far away from the coexistence curve. A rough estimate for the field strength beyond which the simple droplet theory discussed here should break down (or at least become suspect) is obtained by requiring that 1 [139]. The resulting crossover field is sometimes called the “mean-field spinodal point”, or MFSP [108, 110], but it is not identical to the sharp spinodal found when the Ising ferromagnet is treated in the mean-field approximation. The MFSP is located at . The field region beyond this limit we call the strong-field region. It will not be discussed further here, but we hope to return to it in the future.

### 6.2 Growth and the KJMA Theory

The reasoning behind the Kolmogorov-Johnson-Mehl-Avrami (KJMA) theory is quite simple [96–98, 175–178]. Critical droplets are assumed to nucleate in the metastable phase and subsequently grow without substantial deformation [179]. (Recent discussions relevant to the limitations of the latter assumption can be found in Refs. [180–183].) The nucleation rate may either be constant (homogeneous nucleation), or all the nuclei may already be present at =0 (heterogeneous nucleation [33]). The KJMA theory is simple to work out for both cases [33, 68], but only the former will be pursued here.

The radial growth velocity, which approaches a constant limit for large droplets, is obtained in an “Allen-Cahn” approximation [72, 73, 175–178, 184–186] as

(25) |

where the coefficient depends on the details of the kinetics. (To avoid confusion, we note that the usual growth law for spinodal decomposition in systems with nonconserved order parameter, , results from setting =0 (i.e., =0) in Eq. (25). The growth law obtained in the case of metastable decay is a consequence of the difference in bulk free energy between the two phases, which acts as a driving force for the growth [73, 184, 185].)

If we set =, neglect the volumes of the critical droplets, and consider an uncorrelated “ideal gas” of freely overlapping domains of the stable phase, then the volume fraction of stable phase at time is

(26) |

where the integration variable is the time at which a particular supercritical droplet was nucleated. This relation is known as Avrami’s law. The timescale is the characteristic time for collisionless growth and sets the basic timescale for the decay of the metastable phase. By performing the integration in Eq. (26), one sees that is given by

(27) |

where is a nonuniversal function of . The second equality in Eq. (27) was obtained by using Eq. (22) for the homogeneous nucleation rate , neglecting for simplicity the correction term in the exponential. By comparing Eq. (27) for with Eq. (22) for , we notice that, apart from the pre-exponential factors, the characteristic time is simply given by the nucleation rate to the power +1) [132].

Associated with is the characteristic lengthscale for collisionless growth, which we loosely call the mean droplet separation. It is given by

(28) |

where is a nonuniversal function of . The lengthscale can be very large for weak fields, and we note that, even though the critical droplet radius , in the same limit. as

The KJMA theory gives a good approximation for sufficiently small (well below the percolation limit [147]) that droplet correlations do not significantly alter the growth, and it quickly became popular for analyzing experimental results on metastable decay in a number of areas. These include situations in which is smaller than the dimension of the physical space and must be interpreted as the dimension of the subspace in which the droplets grow [98]. However, the theory does not seem to have attracted sustained attention from theorists until the 1980’s, when Sekimoto derived exact expressions for the space-time correlation function and the structure factor for the KJMA process [175–178]. These results were later generalized to infinitely [187] and multiply [188] degenerate equilibrium phases, and mean-field results for the multiply degenerate case with multiple nucleation and growth rates have been obtained, both with homogeneous and heterogeneous nucleation [189, 190]. Applications of Eq. (26) have recently been made to theoretical studies of metastable decay in long-range interaction models and ferroelectrics [35–38], to nanometer-sized ferromagnetic particles [49], and to kinetic Ising models [37, 38, 49, 108]. Several of these papers contain kinetic MC simulations [37, 38, 49, 108, 187, 189, 190].

### 6.3 Finite-Size Effects

In an infinitely large system, the number of critical and supercritical droplets is of course infinite, and the decay of the order parameter is well described by Avrami’s law, Eq. (26). In this section we consider how finite system size modifies the KJMA picture. Although finite-size scaling has been extremely useful in the study of equilibrium critical phenomena (see, e.g., Ref. [191]), systematic finite-size analysis of metastable decay seems to have been performed only recently [36, 37, 108, 110]. The discussion below follows closely that of Ref. [108]. For simplicity we retain our exclusive focus on hypercubic systems of volume with periodic boundary conditions.

For temperatures well below the critical temperature and fields inside the MFSP, the correlation lengths in both the stable and the metastable phase are microscopic. We are then left to consider the interplay between three lengths: the system size , the mean droplet separation , and the critical radius .

In the large- limit, where

(29) |

the system can be approximately partitioned into cells of volume . Each cell decays in an independent Poisson process of rate . The volume fraction is then self-averaging, so that Eq. (26) can be inverted to yield the average time it takes for the volume fraction of the equilibrium phase to increase to a given value ,

(30) |

The relative standard deviation of is

(31) |

where 1 is the relative standard deviation of a single Poisson process (not to be confused with the lattice-gas density defined in Eq. (5)).

The regime characterized by has been termed “the deterministic region” [108, 110]. It is subdivided into the strong-field region mentioned in Sec. 6.1, and a region characterized by a finite density of growing droplets, which we call the multi-droplet region [108]. Observations of the deterministic region in MC simulations are also indicated in Refs. [37, 38, 49, 105, 108, 122, 126, 132]. The characteristic absence of -dependence in the multi-droplet regime was noted by Binder and Müller-Krumbhaar, who also derived equations equivalent to Eqs. (26) and (30) [122]. Although the main emphasis was on the nucleation process, the multi-droplet picture was also implied in Langer’s work [76–78].

For smaller , so that

(32) |

the random nucleation of a single critical droplet in a Poisson process of rate is the rate-determining step. This is followed by relatively rapid growth, until this droplet occupies the entire system after an additional time much shorter than the average waiting time before a second droplet nucleates. Therefore, the characteristic lifetime becomes

(33) |

In this case , and depends only weakly on the threshold . This single-droplet region is part of “the stochastic region” identified in Ref. [110], and it was detected in MC simulations in Refs. [37, 38, 49, 105, 108, 126, 127, 132].

The crossover between the deterministic and stochastic regimes is determined by the condition with a proportionality constant of order unity. We identify the crossover field with the “dynamic spinodal point” (DSP) introduced in Ref. [110], and in the limit we explicitly obtain from Eq. (28)

(34) |

This crossover field was observed in MC simulations reported in Refs. [37, 38, 49, 105, 108, 126, 132]. Following Refs. [108, 110], we estimate as the field where the relative standard deviation for the lifetime is =1/2. We emphasize that, although vanishes as , the approach to zero is exceedingly slow, especially for =3 and above. Therefore, may well be measurably different from zero for systems that are definitely macroscopic as far as their equilibrium properties are concerned. (As an illustration, increasing from 100 to 10 for =3 decreases the leading term in only to approximately one-half of its original value!)

Finally, we consider the small- limit,

(35) |

In this case the volume term can be neglected in Eq. (19), and the free-energy cost of a droplet occupying a volume fraction is with a proportionality constant between 1 and 1/, so that the first-passage time to a given is independent of and diverges exponentially with . Since the dynamics in this region of extremely weak fields or extremely small systems is similar to that on the coexistence line, =0 [192], we call it the coexistence region. The crossover field between the coexistence and single-droplet regions, called “the thermodynamic spinodal point” (THSP) in Ref. [110], is determined for a given by , which yields

(36) |

This crossover field was observed in MC simulations reported in Refs. [49, 108, 127].

In summary, by comparing the characteristic lengths and with the lattice constant and the system size , one can identify four different field regions, in which the decay proceeds through different excitations. In order of increasingly strong unfavorable field , these are the “coexistence region,” characterized by subcritical fluctuations on the scale of the system volume; the “single-droplet region,” characterized by decay via a single critical droplet; the “multi-droplet region,” characterized by decay via a finite density of droplets; and the “strong-field region,” in which the droplet picture is inappropriate. The crossover fields between these regions,

(37) |

are accurately predicted by droplet theory. The different regions and crossover fields are illustrated in Figs. 3–6.

In Fig. 3 are shown average metastable lifetimes for square Ising ferromagnets with =128 and 720 at =0.8. The data points were obtained by the Metropolis algorithm with updates at randomly chosen sites [108] and have here been extrapolated past for =128 (given by Eq. (36)), using Eq. (33) for the lifetime in the single-droplet region. The four field regions can clearly be distinguished.

An alternative view of the information contained in Fig. 3 is found in Fig. 4, which shows as a function of the field at which =1/2), plotted for two different values of at =0.8 [49]. This figure corresponds to a contour plot of data like those shown in Fig. 3. In accordance with usage in the experimental literature on small magnetic particles [41–48], we call the “switching field.” For qualitative comparison we have also included in Fig. 4 data digitized from Fig. 5 of Ref. [45], which shows the effective switching field (corrected for the demagnetization field) vs. the particle diameter for single-domain ferromagnetic barium ferrite particles, measured by magnetic-force microscopy at room temperature. Considering the different dimensionalities of the model and the experimental system, and that no particular effort was made to fit the parameters in the Ising model to the experiments, we find the similarity between the simulated and the experimental switching fields striking. This qualitative agreement may indicate that the decay of the metastable magnetization state in the barium ferrite particles proceeds through similar nucleation and growth mechanisms as in the Ising model, and it should be relevant to the current debate over the magnetization reversal mode in single-domain ferromagnetic particles [42–48].

In both the multi-droplet and the single-droplet regions, the metastable lifetime (determined by Eqs. (27) and (33), respectively) has the form of an exponential in multiplied by a power-law prefactor in . In both regions the derivative of with respect to can therefore be written as

(38) |

with

Due to the dramatic increase in the metastable lifetime as is lowered, numerical data confirming the theoretical predictions at lower temperatures must be obtained by other techniques, such as the CTM method or one of the new MC methods discussed in Sec. 3. In Fig. 6 we show the temperature dependence of for systems with =24 and 240, as obtained by the new MC with absorbing Markov Chains (MCAMC) method [113–115], together with our analytic estimate for . The slow decay of with predicted by Eq. (34) can be seen. Also note the dramatic widening of the single-droplet region as is lowered for a system of fixed size, in agreement with recent exact predictions [106, 193, 194].

Confirmation that the quantity is given by its zero-field equilibrium value is given in Fig. 7 (after Ref. [195]), which shows estimates of for between 0.17 (=0.4) and 0.8 based on the CTM method [138, 139], MCAMC simulations [115], and standard Metropolis MC simulations [108]. The CTM estimates were obtained by fitting the logarithmic derivative of the metastable to Eq. (38) with =1, =0, and the correction term from Eq. (22) included. The numerical estimates for obtained by the different methods agree very well with each other, as well as with the exact equilibrium result.

### 6.4 Discrete-Lattice Effects

Recently, a series of rigorous papers [106, 193, 194, 196–198] have appeared that discuss effects of the lattice discreteness on the metastable lifetimes in the limit 0. The approach is to consider a birth-death process describing the time dependence of the fluctuations corresponding to subcritical droplets in the two-dimensional Ising ferromagnet with single-spin-flip Metropolis [193, 194, 196, 197] or modified Swendsen-Wang [106] dynamics. For isotropic interactions [106, 193, 194, 196] and 2 it is shown that the metastable phase almost certainly decays through a single “proto-critical” droplet of spins pointing parallel to the field. This droplet is shaped like a rectangle of 1) overturned spins, with one additional overturned spin attached as a “knob” to one of its long sides. The length = is the smallest integer larger than , where is the zero-temperature limit of the diameter of the critical droplet. In the zero-temperature limit the average metastable lifetime is found to be a piecewise linear function in which diverges asymptotically as for small . It is given by

(39) | |||||

These results are significant for two reasons. First, they provide an explicit dynamical justification for the existence of a critical droplet size in kinetic Ising models. Second, they generalize the standard continuum droplet theory outlined in Sec. 6.1 to consider the discreteness of the lattice. Using the exact zero-temperature value, =, one sees that Eq. (39) is consistent to leading order in with the continuum droplet-theory result for , Eqs. (22) and (23).

The temperatures at which these discrete-lattice results are expected to be valid are so low as to be inaccessible with standard MC techniques. However, they are readily accessible, both with the MCAMC method [113–115] and with the CTM method [138, 139]. Results obtained with these methods are shown in Fig. 8, together with the derivative of the first line in Eq. (39) with respect to and the corresponding quantity in continuum droplet theory, from Eq. (38). The discrete-lattice results are represented by the series of solid, parabolic arcs, and the continuum results for the single-droplet region are represented by the two straight lines. The dashed line corresponds to =1 and =2, which is appropriate for the MC results, whereas the dotted line corresponds to =1 and =0 and is appropriate for the CTM results which do not contain a kinetic prefactor. For relatively strong fields, the MC data agree reasonably well with the discrete-lattice predictions. For weaker fields, where the critical droplets become larger, the oscillations caused by the lattice discreteness become less pronounced, and the MC data points appear to approach the continuum result. The CTM method allows a closer approach to =0 than the MC, but the deviations from the continuum result are larger than for the MC at the same field. A detailed comparison of the manner in which the results obtained by the two methods approach their respective continuum limits, represented by the two straight lines which intersect at the common zero-field limit , has yet to be performed. It is likely to involve both the difference between the geometries of the two systems studied (square for MC and infinite strip for CTM) and the fact that the CTM quantities can be evaluated only at particular fields determined by the lobe structure of the metastable [138, 139], which is similar to the one shown for the WLRF Q1DI model in Fig. 1.

The theoretical and numerical results discussed in this subsection raise a number of questions that should be answerable using the new MC and TM algorithms. In particular it is important to understand how the lifetimes cross over to the well confirmed results of continuum nucleation theory at higher temperatures and weaker fields, and also how the lifetimes and the crossover to continuum theory are affected by a change from Metropolis to Glauber dynamics. Furthermore, some rather dramatic effects have been predicted for the case of anisotropic interactions as 0 [197]. No sign of these effects were found in a recent study by the CTM method [144], and a further investigation by the MCAMC method is in progress.

## 7 Summary and Discussion

In this review we have presented some aspects of the current knowledge regarding the mechanisms and rates of decay of metastable phases in statistical-mechanical systems. With its record of 270 years of published research in a variety of basic and applied contexts following Fahrenheit’s 1724 paper [1], this is indeed a venerable field of inquiry. Yet we believe it is a sign of continued vigor and relevance that nearly half of the references we cite have appeared within the last decade.

It is now well understood that metastable decay is a kinetic phenomenon, whose rate depends on the interplay between the (homogeneous and/or heterogeneous) nucleation of critical droplets of the equilibrium phase and on the subsequent growth of the supercritical droplets. Formally, the nucleation rate is given by the Van’t Hoff-Arrhenius relation of classical nucleation theory, Eq. (6), or the field-theoretical relation, Eq. (7) [76–78], and the simultaneous nucleation and growth give rise to “Avrami’s law”, Eq. (26) [96–98]. However, both the free energy of the critical droplet and the pre-exponential factor in Eq. (6) (both of which are contained in the imaginary part of the metastable free energy, Im, in Eq. (7)) depend crucially on the geometry of the critical droplet. Whereas the shape of the critical droplet determines the exponential factor in Eq. (6), the prefactor is determined by the fluctuations around the critical droplet shape. The latter is evaluated in the “post-classical” field-theoretical approach by a steepest-descent calculation around a saddle point representing the critical droplet.

If the order parameter is constrained to be uniform over the entire system, so that becomes extensive in the system size, one recovers the mean-field picture of van der Waals [2], Maxwell [3], and Néel [39]. In this picture the metastable phases are infinitely long-lived in the limit of infinite system size, and the limit of metastability is marked by a sharp spinodal line.

A different class of systems are the weak-long-range-force (WLRF) models discussed in Sec. 5, which may be useful models for, e.g., superconductors, long-chain polymers, and systems with elastic interactions. These models also have infinitely long-lived metastability and a sharp spinodal in the limit of infinite interaction range [58, 59, 149]. However, the critical droplet is a ramified object [92, 93, 153–155] whose free energy is given in terms of the interaction range and the distance from the spinodal by Eq. (14) [92, 93], and the resulting nucleation rate is given by Eq. (17). The agreement between these analytical results and the numerical constrained-transfer-matrix (CTM) method is shown in Fig. 2 [94], which illustrates the dramatic increase in metastable lifetime as the distance from the mean-field spinodal increases.

For short-range-force (SRF) models, which are useful to represent, e.g., anisotropic magnets governed by exchange interactions, and which belong to the same static universality class as the liquid/vapor phase transition, the situation is again different. In this case the critical droplet is a compact object with a radius given by Eq. (20), and is given by Eq. (21). The nucleation rate from the field-theoretical saddle-point calculation is given in Eqs. (22)–(24) [76–78, 90], and the effects of simultaneous nucleation and growth are reflected in the explicit form of Avrami’s law, Eqs. (26) and (27) [96–98]. In contrast to mean-field and WLRF systems, SRF systems do not have a sharp spinodal.

Finite-size effects represent an aspect of metastable kinetics which has attracted increased attention in recent years [36, 37, 108, 110]. For SRF systems, these effects give rise to different dependences of the metastable lifetime on applied field (or chemical potential, pressure, or supersaturation) and system size for different values of these parameters. These regions are separated by size dependent crossover fields: the thermodynamic spinodal field given in Eq. (36), the dynamic spinodal field given in Eq. (34), and the size-independent “mean-field” spinodal field . For single droplet, and the mean lifetime is given by the inverse nucleation rate through Eq. (33). In contrast, for nonzero density of simultaneously nucleating and growing droplets, and the mean lifetime is given by Avrami’s law through Eqs. (27) and (30). The detailed agreement between these analytical results and Monte Carlo (MC) simulations is illustrated in Fig. 3, which shows the dramatic dependence on 1/ of the logarithm of the average metastable lifetime for two-dimensional systems of different sizes. Similar agreement has also been demonstrated for three-dimensional systems [130–132], and a rigorous justification for the nucleation-theory results has been obtained in two dimensions at very low temperatures [106, 193, 194, 196–198]. Further agreement between analytic theory and simulation results is illustrated in Figs. 5, 7, and 8. the decay involves a the decay occurs through a

A number of questions regarding metastable decay still remain to be answered. Probably some of the most difficult ones are related to fluids, which we discussed briefly in Sec. 4. In particular we would like to see developed a droplet definition as unambiguous as the “Fortuin-Kastelyn-Coniglio-Klein-Swendsen-Wang” definition [104–106, 145–147] currently used for Ising and lattice-gas systems.

Some of the most exciting potential for both fundamental and technological progress may lie in combining modern “atomic engineering” techniques, such as atomic and magnetic force microscopies, with computer simulation. Using these experimental techniques, one could both build and study well-characterized metastable systems. Employing state-of-the-art computer hardware and algorithms, one could simulate these same systems with tens to billions [109, 132] of particles. This approach should provide the opportunity to study effects of size, boundary conditions, and impurities by a combination of theory, simulation, and experiment in a carefully controlled manner. A step in this direction was recently taken by Richards et al. [49], some of whose results are shown together with experimental data from Ref. [45] in Fig. 4.

In summary, we think the field of metastability is a healthy 270-year-old with a long, exciting life ahead of it.

Acknowledgments

P.A.R. dedicates this article to the memory of his father, Per Rikvold (1919–1994).

We gratefully acknowledge discussions with M. A. Novotny, C. C. A. Günther, H. L. Richards, and J. Lee, correspondence with J. Lothe, K. Nishioka, H. Reiss, and D. Stauffer, and permission to use published and unpublished data by M. A. Novotny, C. C. A. Günther, H. L. Richards, and S. W. Sides. This work was supported by Florida State University through the Supercomputer Computations Research Institute (U.S. Department of Energy Contract No. DE-FC05-85ER25000) and the Center for Materials Research and Technology, and by the U.S. National Science Foundation Grants No. DMR-9013107 and DMR-9315969.

## References

- [1] D. B. Fahrenheit, Phil. Trans. Roy. Soc. Lond. 33, 78 (1724), in Latin. German translation in Abhandlungen über Thermometrie von Fahrenheit, Réaumur, Celsius, edited by A. J. von Oettingen, Ostwald’s Klassiker der exakten Wissenschaften, No. 57 (Wilhelm Engelmann, Leipzig, 1894), p. 6.
- [2] J. D. van der Waals, Ph.D. Dissertation (University of Leiden, Leiden, 1873).
- [3] J. C. Maxwell (1875), in The Scientific Papers of James Clerk Maxwell, edited by W. D. Niven (Cambridge University Press, Cambridge, 1890), Vol. II, p. 425. Reprinted by Dover, New York, 1965.
- [4] J. W. Gibbs, Trans. Conn. Acad. Sci. 3, 108 (1876); 3 343 (1878). Reprinted in Ref. [5].
- [5] J. W. Gibbs, The Scientific Papers of J. Willard Gibbs (Longmans Green, London, 1906), Vol. I, pp. 252–258. Reprinted by Dover, New York, 1961.
- [6] J. Rice, in A Commentary on the Scientific Writings of J. Willard Gibbs, edited by F. G. Donnan and A. Haas (Yale University Press, New Haven, CT, 1936), Vol. I, pp. 625–631. Reprinted by Arno Press, New York, 1980.
- [7] M. Volmer and A. Weber, Z. Phys. Chem. (Leipzig) 119, 277 (1926).
- [8] Z. Farkas, Z. Phys. Chem. (Leipzig) A125, 236 (1927).
- [9] R. Becker and W. Döring, Ann. Phys. (Leipzig) 24, 719 (1935).
- [10] J. B. Zeldovich, Acta Physicochim (U.R.S.S.) 18, 1 (1943).
- [11] I. Frenkel, Kinetic Theory of Liquids (Oxford University Press, London, 1946), Ch. VII. Reprinted by Dover, New York, 1965.
- [12] A. Bijl, Discontinuities in the Energy and Specific Heat Ph.D. Dissertation (University of Leiden, Leiden, 1938).
- [13] J. Frenkel, J. Chem. Phys. 7, 200 (1939).
- [14] W. Band, J. Chem. Phys. 7, 324; 927 (1939).
- [15] J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 (1958).
- [16] J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 31, 688 (1959).
- [17] J. E. McDonald, Am. J. Phys. 30, 870 (1962); 31, 31 (1963). Reprinted in Ref. [18].
- [18] F. F. Abraham, Homogeneous Nucleation Theory (Academic, New York, 1974).
- [19] W. J. Dunning, in Nucleation, edited by A. C. Zettlemoyer (Marcel Dekker, New York, 1969), p. 1.
- [20] J. Lothe and G. M. Pound, J. Chem. Phys. 36, 2080 (1962).
- [21] J. Feder, K. C. Russell, J. Lothe, and G. M. Pound, Adv. Phys. 15, 111 (1966).
- [22] J. Lothe and G. M. Pound, in Nucleation, edited by A. C. Zettlemoyer (Marcel Dekker, New York, 1969), p. 109.
- [23] K. Nishioka, in Nucleation Phenomena, edited by A. C. Zettlemoyer (Elsevier, Amsterdam, 1977), p. 205.
- [24] K. Nishioka, Phys. Scr. T44, 23 (1992).
- [25] V. Talanquer and D. W. Oxtoby, J. Chem. Phys. 100, 5190 (1994).
- [26] W. C. Swope and H. C. Andersen, Phys. Rev. B 41, 7042 (1990).
- [27] H. M. Ellerby, C. L. Weakliem, and H. Reiss, J. Chem. Phys. 95, 9209 (1991).
- [28] H. M. Ellerby and H. Reiss, J. Chem. Phys. 97, 5766 (1992).
- [29] C. L. Weakliem and H. Reiss, J. Chem. Phys. 99, 5374 (1993).
- [30] J. S. van Duijneveldt and D. Frenkel, J. Chem. Phys. 96, 4655 (1992).
- [31] D. W. Oxtoby, J. Phys.: Condens. Matter 4, 7627 (1992).
- [32] W. N. Cottingham, D. Kalafatis, and R. Vinh Mau, Phys. Rev. B 48, 6788 (1993).
- [33] Y. Ishibashi and Y. Takagi, J. Phys. Soc. Jpn. 31, 506 (1971).
- [34] P. B. Littlewood and P. Chandra, Phys. Rev. Lett. 57, 2415 (1986).
- [35] P. Chandra, Phys. Rev. A 39, 3672 (1989).
- [36] H. Orihara and Y. Ishibashi, J. Phys. Soc. Jpn. 61, 1919 (1992).
- [37] H. M. Duiker and P. D. Beale, Phys. Rev. B 41, 490 (1990).
- [38] P. D. Beale, Integrated Ferroelectrics 4, 107 (1994).
- [39] L. Néel, Ann. Géophys. 5, 99 (1949).
- [40] E. F. Kneller and F. E. Luborsky, J. Appl. Phys. 34, 656 (1963).
- [41] E. Köster and T. C. Arnoldussen, in Magnetic Recording. Volume I: Technology, edited by C. D. Mee and E. D. Daniel (McGraw-Hill, New York, 1987), p. 98.
- [42] A. Aharoni, in Magnetic Properties of Fine Particles, edited by J. L. Dormann and D. Fiorani (Elsevier, Amsterdam, 1992), p. 3.
- [43] H.-B. Braun, Phys. Rev. Lett. 71, 3557 (1993).
- [44] H.-B. Braun, preprint (1994).
- [45] T. Chang, J.-G. Zhu, and J. H. Judy, J. Appl. Phys. 73, 6716 (1993).
- [46] M. Lederman, G. A. Gibson, and S. Schultz, J. Appl. Phys. 73, 6961 (1993).
- [47] M. Lederman, D. R. Fredkin, R. O’Barr, S. Schultz, and M. Ozaki, J. Appl. Phys. 75, 6217 (1994).
- [48] C. Salling, R. O’Barr, S. Schultz, I. McFadyen, and M. Ozaki, J. Appl. Phys. 75, 7989 (1994).
- [49] H. L. Richards, S. W. Sides, P. A. Rikvold, and M. A. Novotny, in preparation.
- [50] Y. N. Ovchinnikov and I. M. Sigal, Phys. Rev. B 48, 1085 (1993).
- [51] K. Kajantie, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 149.
- [52] M. Hackel, M. Faber, M. Markum, and M. Müller, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 189.
- [53] S. Coleman, Phys. Rev. D 15, 2929 (1977).
- [54] C. G. Callan and S. Coleman, Phys. Rev. D 16, 1762 (1977).
- [55] M. Gleiser and E. Kolb, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 1.
- [56] J. Kripfganz, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 11.
- [57] W. Buchmüller and T. Helbig, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 27.
- [58] O. Penrose and J. L. Lebowitz, J. Stat. Phys. 3, 211 (1971).
- [59] O. Penrose and J. L. Lebowitz, in Fluctuation Phenomena, edited by E. W. Montroll and J. L. Lebowitz (North-Holland, Amsterdam, 1979), Chap. 5, p. 293.
- [60] G. L. Sewell, Phys. Repts. 57, 308 (1980).
- [61] P. C. Hohenberg and B. Halperin, Rev. Mod. Phys. 49, 435 (1977).
- [62] M. Acharyya and B. K. Chakrabarti, in this volume.
- [63] See, e.g., H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena (Oxford University Press, Oxford, 1971).
- [64] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
- [65] R. J. Glauber, J. Math. Phys. 4, 294 (1963).
- [66] Ph. A. Martin, J. Stat. Phys. 16, 149 (1977).
- [67] V. A. Shneidman and P. Hänggi, J. Chem. Phys. (1994), in press.
- [68] U. R. Evans, Trans. Faraday Soc. 41, 365 (1945).
- [69] J. E. Mayer, J. Chem. Phys. 5, 67 (1937).
- [70] J. H. Van’t Hoff, Etudes de Dynamiques Chimiques (F. Muller and Co., Amsterdam, 1884), p. 114.
- [71] S. Arrhenius, Z. Phys. Chem. (Leipzig) 4, 226 (1889).
- [72] J. D. Gunton and M. Droz, Introduction to the Theory of Metastable and Unstable States (Springer, Berlin, 1983).
- [73] J. D. Gunton, M. San Miguel, and P. S. Sahni, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic, New York, 1983), Vol. 8.
- [74] K. Binder, Rep. Prog. Phys. 50, 783 (1987).
- [75] M. E. Fisher, in Proceedings of the Gibbs Symposium, Yale University, May 15–17, 1989, edited by G. D. Mostow and D. G. Caldi (American Mathematical Society, Providence, RI, 1990).
- [76] J. S. Langer, Ann. Phys. (N.Y.) 41, 108 (1967).
- [77] J. S. Langer, Phys. Rev. Lett. 21, 973 (1968).
- [78] J. S. Langer, Ann. Phys. (N.Y.) 54, 258 (1969).
- [79] A. F. Andreev, Sov. Phys. JETP 18, 1415 (1964).
- [80] M. E. Fisher, Physics 3, 255 (1967).
- [81] R. Landauer and J. A. Swanson, Phys. Rev. 121, 1668 (1961).
- [82] H. A. Kramers, Physica 7, 284 (1940).
- [83] C. M. Newman and L. S. Schulman, J. Stat. Phys. 23, 131 (1980).
- [84] G. Roepstorf and L. S. Schulman, J. Stat. Phys. 34, 35 (1984).
- [85] B. Gaveau and L. S. Schulman, Lett. Math. Phys. 18, 201 (1989).
- [86] O. Penrose, J. Stat. Phys. (1994), in press.
- [87] M. Büttiker and R. Landauer, Phys. Rev. Lett. 43, 1453 (1979).
- [88] M. Büttiker and R. Landauer, Phys. Rev. A 23, 1397 (1981).
- [89] R. J. McCraw, Phys. Lett. A 75, 379 (1980).
- [90] N. J. Günther, D. A. Nicole, and D. J. Wallace, J. Phys. A 13, 1755 (1980).
- [91] W. Zwerger, J. Phys. A 18, 2079 (1985).
- [92] W. Klein and C. Unger, Phys. Rev. B 28, 445 (1983).
- [93] C. Unger and W. Klein, Phys. Rev. B 29, 2698 (1984).
- [94] B. M. Gorman, P. A. Rikvold, and M. A. Novotny, Phys. Rev. E 49, 2711 (1994).
- [95] T. Fiig, B. M. Gorman, P. A. Rikvold, and M. A. Novotny, Phys. Rev. E (1994), in press.
- [96] A. N. Kolmogorov, Bull. Acad. Sci. USSR, Phys. Ser. 1, 355 (1937).
- [97] W. A. Johnson and P. A. Mehl, Trans. Am. Inst. Mining and Metallurgical Engineers 135, 416 (1939).
- [98] M. Avrami, J. Chem. Phys. 7, 1103 (1939); 8, 212 (1940); 9, 177 (1941).
- [99] H. Müller-Krumbhaar and K. Binder, J. Stat. Phys. 8, 1 (1973).
- [100] R. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
- [101] K. Binder, J. Comput. Phys. 59, 1 (1985).
- [102] K. Binder, in Monte Carlo Methods in Statistical Physics, Second Edition, edited by K. Binder (Springer, Berlin, 1986).
- [103] R. M. Neal, University of Toronto Department of Computer Science Technical Report CRG-TR-93-1 (1993). This review with annotated bibliography is available on the Internet from .
- [104] T. S. Ray and P. Tamayo, J. Stat. Phys. 60, 851 (1990).
- [105] T. S. Ray and J.-S. Wang, Physica A 167, 580 (1990).
- [106] F. Martinelli, E. Olivieri, and E. Scoppola, J. Stat. Phys. 62, 135 (1991).
- [107] K. Kawasaki, in Phase Transitions and Critical Phenomena, edited by C. Domb and M. Green (Academic, London, 1972), Vol. 2.
- [108] P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides, Phys. Rev. E 49, 5080 (1994).
- [109] N. Ito, Physica A 192, 604 (1993).
- [110] H. Tomita and S. Miyashita, Phys. Rev. B 46, 8886 (1992).
- [111] H. Tomita and S. Miyashita, in Computational Approaches in Condensed-Matter Physics, edited by S. Miyashita, M. Imada, and H. Takayama (Springer, Berlin, 1992), p. 278.
- [112] M. Iosifescu, Finite Markov Processes and their Application (Wiley, New York, 1980), p. 99.
- [113] M. A. Novotny, in Computer Simulation Studies in Condensed Matter Physics VII, edited by D. P. Landau, K. K. Mon, and B. Schüttler (Springer, Berlin, in press).
- [114] M. A. Novotny, submitted to Phys. Rev. Lett.
- [115] M. A. Novotny, in preparation.
- [116] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comput. Phys. 17, 10 (1975).
- [117] J. Lee, M. A. Novotny, and P. A. Rikvold, in preparation.
- [118] B. A. Berg and T. Celik, Phys. Rev. Lett. 69, 2292 (1992), and references cited therein.
- [119] K. Binder, Phys. Rev. B 8, 3423 (1973).
- [120] E. Stoll and T. Schneider, Phys. Rev. A 6, 429 (1972).
- [121] K. Binder and E. Stoll, Phys. Rev. Lett. 31, 47 (1973).
- [122] K. Binder and H. Müller-Krumbhaar, Phys. Rev. B 9, 2328 (1974).
- [123] A review of early MC work on metastable decay is given by K. Binder, in Phase Transitions and Critical Phenomena, edited by C. Domb and M. Green (Academic, London, 1976), Vol. 5B.
- [124] K. Binder, Ann. Phys. (N.Y.) 98, 390 (1976).
- [125] K. Binder and D. Stauffer, Adv. Phys. 25, 343 (1976).
- [126] E. Stoll and T. Schneider, Physica 86–88B, 1419 (1977).
- [127] R. J. McCraw and L. S. Schulman, J. Stat. Phys. 18, 293 (1978).
- [128] K. Binder and M. H. Kalos, J. Stat. Phys. 22, 363 (1980).
- [129] W. Paul and D. W. Heermann, Europhys. Lett. 6, 701 (1988).
- [130] D. Stauffer, A. Coniglio, and D. W. Heermann, Phys. Rev. Lett. 49, 1299 (1982).
- [131] D. W. Heermann, A. Coniglio, W. Klein, and D. Stauffer, J. Stat. Phys. 36, 447 (1984).
- [132] D. Stauffer, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 287.
- [133] T. S. Ray, J. Stat. Phys. 62, 463 (1991).
- [134] D. W. Heermann, W. Klein, and D. Stauffer, Phys. Rev. Lett. 49, 1262 (1982).
- [135] W. Paul, D. W. Heermann, and K. Binder, J. Phys. A 22, 3325 (1989).
- [136] P. A. Rikvold, Prog. Theor. Phys. Suppl. 99, 95 (1989).
- [137] C. Domb, Adv. Phys. 9, 149 (1960).
- [138] C. C. A. Günther, P. A. Rikvold, and M. A. Novotny, Phys. Rev. Lett. 71, 3898 (1993).
- [139] C. C. A. Günther, P. A. Rikvold, and M. A. Novotny, submitted to Physica A.
- [140] V. Privman and L. S. Schulman, J. Phys. A 15, L231 (1982).
- [141] V. Privman and L. S. Schulman, J. Stat. Phys. 31, 205 (1982).
- [142] J. N. Kapur and H. K. Kesavan, Entropy Optimization Principles with Applications (Academic, Boston, 1992).
- [143] P. A. Rikvold, B. M. Gorman, and M. A. Novotny, AIP Conf. Proc. Ser. 256, 549 (1992).
- [144] C. C. A. Günther, P. A. Rikvold, and M. A. Novotny, in preparation.
- [145] H. Gould and W. Klein, Physica D 66, 61 (1993).
- [146] J.-S. Wang, Physica A 161, 249 (1989).
- [147] D. Stauffer and A. Aharony, Introduction to Percolation Theory (Taylor and Francis, London, 1992).
- [148] C. Borgs, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 125, and references cited therein.
- [149] P. C. Hemmer and J. L. Lebowitz, in Phase Transitions and Critical Phenomena, edited by C. Domb and M. Green (Academic, London, 1976), Vol. 5B.
- [150] R. B. Griffiths, C.-Y. Wang, and J. S. Langer, Phys. Rev. 149, 301 (1966).
- [151] J. D. Gunton and M. C. Yalabik, Phys. Rev. B 18, 6199 (1978).
- [152] P. A. Rikvold, B. M. Gorman, and M. A. Novotny, Phys. Rev. E 47, 1474 (1993).
- [153] T. S. Ray and W. Klein, J. Stat. Phys. 61, 891 (1990).
- [154] L. Monette and W. Klein, Phys. Rev. Lett. 68, 2336 (1992).
- [155] L. Monette, W. Klein, M. Zuckermann, A. Khadir, and R. Harris, Phys. Rev. B 38, 11607 (1988).
- [156] M. A. Novotny, W. Klein, and P. A. Rikvold, Phys. Rev. B 33, 7729 (1986).
- [157] W. Ostwald, Lehrbuch der Allgemeinen Chemie (W. Engelmann, Leipzig, 1896-1902), Vol. II, ii, p. 444.
- [158] See, e.g., N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group (Addison-Wesley, Reading, MA, 1992).
- [159] Note that, in contrast to the notation in Ref. [108], in the present work the Hamiltonian is not divided by . As a consequence, the definitions of a number of the quantities discussed in Sec. 6 are different from their definitions in Ref. [108]. The notation used in the present work is the same as the one used in Ref. [139].
- [160] C. K. Harris, J. Phys. A 17, L143 (1984).
- [161] C. Rottman and M. Wortis, Phys. Rev. B 24, 6274 (1981).
- [162] R. K. P. Zia and J. E. Avron, Phys. Rev. B 25, 2042 (1982).
- [163] J. E. Avron, H. van Beijeren, L. S. Schulman, and R. K. P. Zia, J. Phys. A 15, L81 (1982).
- [164] Our quantity is identical to in Eq. (3) of Ref. [162].
- [165] A. D. Bruce and D. J. Wallace, Phys. Rev. Lett. 47, 1743 (1981).
- [166] O. Penrose, Commun. Math. Phys. 124, 515 (1989).
- [167] M. Kreer, Ann. Phys. (Leipzig) 2, 398 (1993).
- [168] M. Kreer, Ann. Phys. (Leipzig) 2, 720 (1993).
- [169] M. J. Lowe and D. J. Wallace, J. Phys. A 13, L381 (1980).
- [170] D. J. Wallace, in Phase Transitions, Proceedings of a Summer Institute, Cargèse, Corsica, 1980, edited by M. Levy, J. C. Le Guillou, and J. Zinn-Justin (Plenum, New York, 1982), p. 423.
- [171] G. Jacucci, A. Perini, and G. Martin, J. Phys. A 16, 369 (1983).
- [172] A. Perini, G. Jacucci, and G. Martin, Phys. Rev. B 29, 2689 (1984).
- [173] A. Perini, G. Jacucci, and G. Martin, Surf. Sci. 144, 53 (1984).
- [174] C. N. Yang, Phys. Rev. 85, 808 (1952).
- [175] K. Sekimoto, Phys. Lett. A 105, 390 (1984).
- [176] K. Sekimoto, J. Phys. Soc. Jpn. 53, 2545 (1984).
- [177] K. Sekimoto, Physica 135A, 328 (1986).
- [178] K. Sekimoto, Int. J. Mod. Phys. B 5, 1843 (1991).
- [179] T. Nakamura, J. Phys. Soc. Jpn. 15, 1379 (1960).
- [180] Y. Saito and T. Ueta, Phys. Rev. A 40, 3408 (1989).
- [181] Y. Saito and T. Ueta, J. Crystal Growth 99, 171 (1990).
- [182] E. Brener, K. Kassner, H. Müller-Krumbhaar, and D. Temkin, in Dynamics of First Order Phase Transitions, edited by H. J. Herrmann, W. Janke, and F. Karsch (World Scientific, Singapore, 1992), p. 53.
- [183] L. Jörgenson, H. Guo, R. Harris, and M. Grant, Phys. Rev. E 48, 4592 (1993).
- [184] I. M. Lifshitz, Sov. Phys. JETP 15, 939 (1962).
- [185] S. K. Chan, J. Chem. Phys. 67, 5755 (1977).
- [186] S. M. Allen and J. W. Cahn, Acta Metall. 27, 1085 (1979).
- [187] J. D. Axe and Y. Yamada, Phys. Rev. B 34, 1599 (1986).
- [188] S. Ohta, T. Ohta, and K. Kawasaki, Physica 140A, 478 (1987).
- [189] R. M. Bradley and P. N. Strenski, Phys. Rev. B 40, 8967 (1989).
- [190] Y. A. Andrienko, N. V. Brilliantov, and P. L. Krapivsky, Phys. Rev. B 45, 2263 (1992).
- [191] Finite-Size Scaling and Numerical Simulation of Statistical Systems, edited by V. Privman (World Scientific, Singapore, 1990).
- [192] U.-J. Wiese, Universität Bern preprint (1993).
- [193] E. Jordão Neves and R. H. Schonmann, Commun. Math. Phys. 137, 209 (1991).
- [194] R. H. Schonmann, Commun. Math. Phys. 147, 231 (1992).
- [195] B. M. Gorman and C. C. A. Günther, in Computer Simulation Studies in Condensed Matter Physics VII, edited by D. P. Landau, K. K. Mon, and B. Schüttler (Springer, Berlin, in press).
- [196] E. Scoppola, J. Stat. Phys. 73, 83 (1993).
- [197] R. Kotecký and E. Olivieri, J. Stat. Phys. 70, 1121 (1993).
- [198] E. Scoppola, Physica A 194, 271 (1993), and references cited therein.