# 1 Introduction

## 1 Introduction

Metastable states are ubiquitous in systems undergoing first order phase transitions. During their lifetime (which can be very long indeed) these states are practically indistinguishable from equilibrium states. Nevertheless, they cannot be described in the framework of the equilibrium Gibbsian formalism . Their analysis in terms of dynamical models has led to a deeper understanding of metastability by providing detailed descriptions of the “escape routes” from metastable to stable states in certain idealized limiting situations.

Following earlier work on systems with long range interactions , the pathwise approach to metastability was introduced in . It was then used in  and  to study rigorously the escape from metastability in the Ising model with nearest neighbor interactions and a small external magnetic field, evolving via Metropolis Glauber dynamics in a finite periodic domain, in the limit of temperature going to zero. The effect of the boundary conditions on the exit path from the metastable phase was analyzed in . Metastability for spin systems with different interactions was investigated in  and in  the problem was considered in a more general context.

The case of finite temperature, infinite volume and external magnetic field going to zero was studied in ; this situation, very interesting from the physical point of view, is mathematically much more complicated than the zero temperature limit. The finite temperature case also has been studied by means of Monte Carlo simulations (see e.g. ) and by transfer-matrix methods ; a clear discussion of these results can be found in .

In all the above works (except for that of Penrose-Lebowitz  dealing with deterministic continuum systems) the spin systems evolve according to a stochastic continuous time or serial dynamics, for which at most one spin of the system is updated at any time. In this paper we investigate metastable behavior in systems with parallel evolution, i.e., in which all the spins of the system are updated simultaneously, at integer times $`t=1,2,3,\mathrm{}.`$ In particular, we are interested in finding how the escape time and escape path from the metastable phase are influenced by the parallel dynamics. A natural setting for this question is that of probabilistic cellular automata (PCA), specifically those whose stationary measures are Gibbs states of Ising models with short range interactions.

PCA were first studied in the Soviet literature of the early seventies  and since then have been applied in many different contexts; in particular, their connections with statistical mechanics were investigated in . In this paper we will consider a PCA for which the dynamics depends on two parameters, $`\beta `$ and $`h`$, and which has the property that its stationary states are Gibbs measures for a certain Hamiltonian $`H(\beta ,h)`$. Here $`\beta `$ plays the role of an inverse temperature and $`h`$ that of an external magnetic field, but the coupling constants in $`H`$ depend on $`\beta `$ in a complicated way; in particular, $`H(\beta ,h)\beta \stackrel{~}{H}(h)`$. As in the standard Ising model we have that when $`\beta `$ is large enough and $`h`$ is zero there exist two different stationary Gibbs measures for the PCA, characterized by nonzero average magnetizations $`\pm m^{}`$, while for $`h0`$ there is a unique stable phase. We then pose the usual question of metastability: if at large $`\beta `$ the system is prepared in the minus (plus) phase and the magnetic field is chosen positive (negative) and small, how does the system reach the stable phase? For definiteness we will always consider the escape from an initial all minus phase.

In spin models with continuous time dynamics an important role in this transition is played by the stable configurations , which are fixed points of the evolution in the limit of zero temperature. For example, a rectangle of pluses of width greater than one inside a sea of minuses is a stable configuration for the nearest neighbor Ising model with a small positive external field. The tendency of such a rectangle to grow or to shrink by the repeated addition or loss of single sites, as a function of its size, yields the behavior of the exit from the metastable state. In the PCA discussed below all configurations are accessible at a single updating. Nevertheless, we will argue that, for large $`\beta `$, the only configurations relevant for the description of the exit from the metastable phase are those in which the plus phase is inside well separated rectangles in the minus sea. We then discuss quantitatively the growth and shrinkage of such droplets and compare our theoretical prediction with results of Monte Carlo simulations.

In Section 2 we define our model and show that it undergoes a phase transition at low temperature, and in Section 3 describe the specific model on which we will focus in the balance of the paper. We discuss our heuristics on the critical behavior of droplets in Section 4, and compare theoretical and Monte Carlo results in Section 5. Section 6 is devoted to some brief conclusions.

## 2 Description of the General PCA Model

Let $`\mathrm{\Lambda }`$ be a $`d`$-dimensional torus containing $`L^d`$ lattice sites, i.e., $`\mathrm{\Lambda }𝐙^𝐝`$ is a cube containing $`L^d`$ points and having periodic boundary conditions. At each site $`x\mathrm{\Lambda }`$ there is spin variable $`\sigma (x)=\pm 1`$; the space $`\{1,1\}^\mathrm{\Lambda }`$ of configurations is denoted by $`\mathrm{\Omega }`$.

To define the dynamics of the model we introduce the discrete time variable $`n=0,1,\mathrm{}`$ and denote by $`\sigma _n`$ the system configuration at time $`n`$. All the spins are updated simultaneously and independently at every unit time; the conditional probability that the spin at site $`x`$ takes value $`\tau `$ at time $`n`$, given the configuration at time $`n1`$, is

$$Prob(\sigma _n(x)=\tau |\sigma _{n1})p_x(\tau |\sigma _{n1})=\frac{1}{2}\left[1+\tau \mathrm{tanh}\beta \left(\underset{y\mathrm{\Lambda }}{}K(xy)\sigma _{n1}(y)+h\right)\right].$$
(1)   
Thus the time evolution is a Markov chain on $`\mathrm{\Omega }`$ with nonzero transition probabilities $`P_\mathrm{\Lambda }(\eta |\sigma )`$ given by

$$P_\mathrm{\Lambda }(\eta |\sigma )=\underset{x\mathrm{\Lambda }}{}p_x\left(\eta (x)|\sigma \right),\sigma ,\eta \mathrm{\Omega }.$$
(2)   
The coupling is of finite range ($`K(z)=0`$ if $`|z|>z_0`$, with $`z_0<L`$ and typically $`z_0L`$) and the coupling constants $`K(z)`$ will be held fixed throughout our discussion. The parameters $`\beta `$ and $`h`$ play the role of inverse temperature and external magnetic field respectively, as discussed above. Note that for large $`|h|`$, $`\sigma _n(x)=sgn(h)`$ with high probability, while for large positive $`\beta `$, $`\sigma _n(x)=sgn[_yK(xy)\sigma _{n1}(y)+h]`$ with high probability.

We say that a probability measure $`\rho (\sigma )`$ on the configuration space $`\mathrm{\Omega }`$ is stationary for the PCA if and only if it remains invariant under the dynamics, i.e., iff $`_\sigma P_\mathrm{\Lambda }(\eta |\sigma )\rho (\sigma )=\rho (\eta )`$. By the general theory of Markov processes there exists, for any $`\beta `$, $`h`$, and $`\mathrm{\Lambda }`$, a unique stationary measure $`\nu _\mathrm{\Lambda }^{\beta ,h}`$ for the PCA. We say that the PCA is reversible with respect to a measure $`\rho `$ iff

$$P_\mathrm{\Lambda }(\eta |\sigma )\rho (\sigma )=P_\mathrm{\Lambda }(\sigma |\eta )\rho (\eta ),\sigma ,\eta \mathrm{\Omega }.$$
(3)   
Summing (3) over $`\sigma `$ shows that any $`\rho `$ satisfying (3) is stationary for the PCA; the converse need of course not be true. It is however straightforward to check that if $`K(z)=K(z)`$ then the process (2) is reversible with respect to the measure

$$\nu _\mathrm{\Lambda }^{\beta ,h}(\sigma )=Z^1\underset{x\mathrm{\Lambda }}{}e^{\beta h\sigma (x)}\mathrm{cosh}[\beta (\underset{y}{}K(xy)\sigma (y)+h)],$$
(4)   
where $`Z`$ is a normalization constant. To see this one simply notes that $`p_x(\eta (x)|\sigma )`$ can be written as

$$p_x(\eta (x)|\sigma )=\frac{1}{2}\frac{\mathrm{exp}[\beta h\eta (x)+\beta \underset{y}{}K(xy)\eta (x)\sigma (y)]}{\mathrm{cosh}\beta [_yK(xy)\sigma (y)+h]}.$$
(5)   
The measure (4) must of course be the unique stationary measure referred to above. From (4) it is clear that $`\nu _\mathrm{\Lambda }^{\beta ,h}(\sigma )`$ is a Gibbs measure for a Hamiltonian $`H(\beta ,h)`$ with (generally many spin) interactions of finite range, which by our assumptions are independent of $`\mathrm{\Lambda }`$:

$$H(\beta ,h)(\sigma )=\beta h\underset{x}{}\sigma (x)\underset{x}{}\mathrm{log}\left[\mathrm{cosh}\beta \left(\underset{y}{}K(xy)\sigma (y)+h\right)\right].$$
(6)   
Hence taking the limit $`\mathrm{\Lambda }𝐙^𝐝`$ yields a Gibbs measure $`\nu ^{\beta ,h}`$ for $`H`$ which is stationary for the PCA on $`𝐙^𝐝`$, defined by the natural extension of the Markov process (2) to $`𝐙^𝐝`$.

The stationary measures for the infinite volume PCA need of course no longer be unique. It is known in general, however, that if one stationary translation invariant (or periodic) measure is Gibbsian then all such measures are Gibbsian for the same Hamiltonian . Hence to find all translation invariant stationary states of our PCA we need only investigate translation invariant Gibbs states for $`H(\beta ,h)`$. Such an investigation begins with the ground states of the Hamiltonian. For the model considered here it is easy to see from (6) that if $`K(z)0`$ for all $`z`$, and if the set of $`K(z)`$ which are nonzero is not chosen in a very special way, then for $`h=0`$ there are exactly two ground states of $`H(\beta ,h)`$, $`+\underset{¯}{1}`$, in which $`\sigma (x)=1`$ for all $`x`$, and $`\underset{¯}{1}`$, in which $`\sigma (x)=1`$ for all $`x`$, while for $`h0`$ there is only one ground state. It then follows from the Pirogov-Sinai theory  that for $`d2`$ and $`\beta `$ sufficiently large there will be in general two extremal translation invariant Gibbs measures for $`h=0`$ and a unique such measure for $`h0`$. By the argument above, the same conclusion holds for stationary states of the infinite volume PCA. We are thus exactly in the same set up as in the familiar ferromagnetic Ising model. We remark that although we are not dealing here with pair interactions, or even exclusively ferromagnetic interactions, it is easy to see that the measures (4) satisfy the both the FKG  and GKS  inequalities.

We may now pose the paradigm question of metastability: if we prepare the system in the starting configuration $`\underset{¯}{1}`$ and take $`h`$ to be small and positive, how quickly and in what manner does the PCA reach its stationary measure? We want to answer this question in the limit of $`\beta \mathrm{}`$, with $`\mathrm{\Lambda }`$ and $`h`$ fixed, in which the stationary state is $`+\underset{¯}{1}`$, and hence may formulate the first part of the problem as that of estimating the first hitting time $`\tau _{+1}=inf\{n0:\sigma _n=+\underset{¯}{1}\}`$, in the limit $`\beta \mathrm{}`$, when the system is prepared in $`\sigma _0=\underset{¯}{1}`$. To answer the second part of the question, we need to describe the path that the system follows to reach $`+\underset{¯}{1}`$; typically, such a path will involve the necessity of passing through one of a small number of critical configurations.

As in the case of continuous (or serial) dynamics, the first problem is to understand the behavior of the “stable” configurations, that is to estimate the probability that a stable configuration grows or shrinks. Rather than discussing this problem in general we shall now focus on a specific model.

## 3 A Special Model

For the rest of this paper we will focus on one special model from among those specified by (1): the two-dimensional model ($`d=2`$) with $`K(z)=1`$ for $`zA_0`$ and $`K(z)=0`$ otherwise, where $`A_0=\{0,\pm e_1,\pm e_2\}`$ is the set consisting of the origin and its four nearest neighbors. Thus the probability distribution of the spin $`\sigma _n(x)`$ is determined by the spins at time $`n1`$ at the five sites in a cross centered at $`x`$. According to (4), the stationary measure $`\nu _\mathrm{\Lambda }^{\beta ,h}`$ of this system will then be

$$\nu _\mathrm{\Lambda }^{\beta ,h}(\sigma )=Z^1(\beta ;\mathrm{\Lambda })\mathrm{exp}[\underset{x\mathrm{\Lambda }}{}U_x(\sigma ;\beta ,h)+\beta h\underset{x\mathrm{\Lambda }}{}\sigma _x],$$
(7)   
where $`U_x(\sigma ;\beta ,h)=U_0(\tau _x\sigma ;\beta ,h)`$ with $`\tau _x`$ the shift operator (with periodic boundary conditions on $`\mathrm{\Lambda }`$) and

$$U_0(\sigma ;\beta ,h)=\underset{AA_0}{}J_{|A|}(\beta ,h)\sigma (A)=\mathrm{log}\mathrm{cosh}(\beta \underset{yA_0}{}\sigma (y)+\beta h),$$
(8)   
with $`\sigma (A)=\mathrm{\Pi }_{yA}\sigma (y)`$ for any $`A\mathrm{\Lambda }`$. The six coefficients $`J_{|A|}(\beta ,h)`$ are determined by the six values which the $`_{yA_0}\sigma (y)`$ can take. For $`h=0`$ only even values of $`|A|`$ occur, and we find

$$U_0(\sigma ;\beta ,0)=J_0(\beta ,0)J_2(\beta ,0)\underset{\{x,y\}A_0}{}\sigma _x\sigma _yJ_4(\beta ,0)\underset{\{x,y,z,w\}A_0}{}\sigma (x)\sigma (y)\sigma (z)\sigma (w),$$
(9)   
with

$`J_0(\beta ,0)`$ $`=`$ $`{\displaystyle \frac{1}{16}}\mathrm{log}[(\mathrm{cosh}5\beta )(\mathrm{cosh}3\beta )^5(\mathrm{cosh}\beta )^{10}]0,`$ (10)   
$`J_2(\beta ,0)`$ $`=`$ $`{\displaystyle \frac{1}{16}}\mathrm{log}[(\mathrm{cosh}5\beta )(\mathrm{cosh}3\beta )/(\mathrm{cosh}\beta )^2]0,`$ (11)   
$`J_4(\beta ,0)`$ $`=`$ $`{\displaystyle \frac{1}{16}}\mathrm{log}[(\mathrm{cosh}5\beta )(\mathrm{cosh}\beta )^2/(\mathrm{cosh}3\beta )^3]0.`$ (12)   
The pair interactions are thus ferromagnetic while the four-spin interactions are not, so the usual conditions for GKS inequalities are not satisfied.

## 4 The Time Evolution

We shall first give a heuristic argument showing that the important configurations for exiting from the metastable state in the PCA are, as for the usual Glauber dynamics, isolated rectangles of pluses of minimum width two. We shall then describe, again on a heuristic level, the growth and shrinkage of one such rectangle. We find a critical value $`l^{}(h)`$, for $`h<1`$, for the length $`l`$ of the smaller side of the rectangle such that, in the limit $`\beta \mathrm{}`$, all rectangles with $`l<l^{}(h)`$ will shrink to zero (except for some special values of $`h`$, for which the condition is $`l<l^{}(h)1`$) while those with $`ll^{}(h)`$ will grow, resulting in an escape from the metastable state.

Let us begin by comparing Glauber dynamics—realized via the Metropolis algorithm, as is usual in questions of metastability—with the PCA dynamics considered above, focussing on differences which are relevant when $`h`$ is small and $`\beta `$ is very large. The former, for a spin system with Hamiltonian $`\stackrel{~}{H}=\stackrel{~}{H}(h)`$ and inverse temperature $`\beta `$, proceeds by spin flips at single sites, with the rate $`c(x;\sigma )`$ at site $`x`$ in configuration $`\sigma `$ given by

$$c(x;\sigma )=\{\begin{array}{cc}1,\hfill & \text{if }\stackrel{~}{H}(\sigma ^x)\stackrel{~}{H}(\sigma )\text{,}\hfill \\ \mathrm{exp}\beta [\stackrel{~}{H}(\sigma ^x)+\stackrel{~}{H}(\sigma )],\hfill & \text{if }\stackrel{~}{H}(\sigma ^x)>\stackrel{~}{H}(\sigma )\text{,}\hfill \end{array}$$
(13)   
where $`\sigma ^x`$ is the configuration obtained from $`\sigma `$ by flipping the spin at site $`x`$:

$$\sigma ^x(y)=\{\begin{array}{cc}\sigma (y),\hfill & \text{for }y=x\text{,}\hfill \\ \sigma (y),\hfill & \text{for }yx\text{.}\hfill \end{array}$$
(14)   
Since $`c(x;\sigma )`$ depends only on $`\beta [\stackrel{~}{H}(\sigma ^x)\stackrel{~}{H}(\sigma )]`$ and is independent of $`\beta `$ if $`\stackrel{~}{H}(\sigma ^x)\stackrel{~}{H}(\sigma )`$, the dynamical landscape is determined entirely by the function $`(\beta H)(\sigma )`$, and the stable configurations, i.e., those invariant under the dynamics in the limit $`\beta \mathrm{}`$, are the local minima of $`\stackrel{~}{H}`$.

In contrast, the PCA dynamics permits transitions from one configuration to any other in a single updating; we will see, however, that this distinction will play only a minor role in the analysis of metastability. Recall that the probability of a transition from $`\sigma `$ to another configuration $`\eta `$ in one time step is given by the product of the probabilities of spin flips at sites where $`\sigma `$ and $`\eta `$ differ with the probabilities of non-flips at sites where they agree. Probabilities of all possible single site flips are shown in Figure 1; it is clear that, at large $`\beta `$, certain flips are almost sure to take place, while all others have exponentially small probability. Thus from an arbitrary initial condition we expect a very rapid evolution to a stable configuration, with further change taking place on an exponentially slow time scale, and involving primarily flips of single spins—in fact, on six different exponential time scales, well separated for $`h<1`$ and very large $`\beta `$, corresponding to the six slow spin flip processes of Figure 1. The parallelism of the dynamics is of relevance during the first, rapid, phase, but the analysis of metastability involves the slow, essentially serial, second phase—although in certain cases we must consider the effect of a small number of unlikely single-site events occurring simultaneously.

A second difference is that, although in the PCA dynamics transitions which lower the energy are generally favored over those which do not, the probability of a transition from $`\sigma `$ to $`\eta `$ is not specified entirely by the energy difference $`H(\eta ;\beta ,h)H(\sigma ;\beta ,h)`$. In particular, for the specific model introduced in Section 3, the single-step probability of flipping a spin which agrees with two of its nearest neighbors is exponentially small, even when such a flip is (energetically) favored by the magnetic field, and hence there are pairs of configurations $`\sigma `$ and $`\eta `$ which differ at a single site but are such that the probability of jumping between them (in either direction) goes to zero as $`\beta \mathrm{}`$. This is illustrated in Figure 1. Consequently there are many more stable configurations for the PCA than for the Glauber dynamics. In fact it is easy to see that any configuration in which the value of the spin at every site agrees with that of at least two of its neighboring sites is stable.

Despite the large number of stable configurations, however, it is rectangular droplets which are important for exit from the metastable state, due to the effect of the most rapid of the “slow” single flip processes of Figure 1. We formalize this as follows. Let $`\widehat{\mathrm{\Omega }}`$ be the set of configurations in which every plus spin agrees with at least two of its nearest neighbors, and define $`T:\widehat{\mathrm{\Omega }}\widehat{\mathrm{\Omega }}`$ so that $`T\sigma `$ is the configuration obtained from $`\sigma `$ by flipping all the minus spins with at least two pluses among their nearest neighbors. For any $`\sigma \widehat{\mathrm{\Omega }}`$ the sequence of configurations $`T^k\sigma `$, $`k=0,1,\mathrm{},`$ is nondecreasing, in the sense that $`(T^{k+1}\sigma )(x)(T^k\sigma )(x)`$ for all $`x`$, and hence must reach a fixed point $`\sigma ^{}`$, in which the set of plus spins forms well separated rectangles (or bands around the torus) inside the sea of minuses. Moreover, if we take $`\sigma `$ as the initial condition $`\sigma _0`$ of the PCA dynamics, and let $`E_\sigma `$ be the event that for some $`n`$, $`\sigma _n=\sigma ^{}`$ and the sequence $`(\sigma _m)`$ is increasing for $`0mn`$, then it is clear from Figure 1 that

$$\underset{\beta \mathrm{}}{lim}ProbE_\sigma =1,$$
(15)   
and that the time to reach $`\sigma ^{}`$ is typically of order $`\mathrm{exp}2\beta (1h)`$. Thus the path for escape from metastability must pass through configurations in which all the pluses are inside well separated rectangles.

We do not attempt to discuss the most general situation but instead consider the fate of a single rectangle; we expect that, as for Glauber dynamics , this is the key element in an analysis of metastability. Let us consider, then, a configuration $`\eta `$ for which all spins inside a rectangle of sides $`l`$ and $`m`$ are up and all other spins are down, and suppose for definiteness that $`lm`$. We will say that such an $`l\times m`$ droplet is supercritical if, starting from $`\eta `$, the system will reach the configuration $`+\underset{¯}{1}`$ before it reaches $`\underset{¯}{1}`$, with probability which approaches 1 in the limit $`\beta \mathrm{}`$; the droplet is subcritical if the reverse is true. We will argue heuristically that if $`l<2/h`$ then the droplet is subcritical and if $`l>2/h`$ the droplet is supercritical, while if $`2/h`$ is an integer and $`l=2/h`$ then the droplet may either shrink or grow. The critical length $`l_h^{}`$ is the smallest integer such that the droplet is supercritical if $`ll_h^{}`$; thus $`l_h^{}=2/h+1`$.

A very rough estimate of $`l_h^{}`$ may be based on energy considerations, with the assumption that if the system starts from a rectangular droplet then the next droplet reached will be one with lower energy. For example, if $`\eta `$ contains an $`l\times l`$ droplet and $`e(l)=H(\eta )`$, one may approximate $`e(l)`$ at very large $`\beta `$ by writing $`\mathrm{log}\mathrm{cosh}\beta x\beta |x|`$ in (8). In this approximation, $`e(l)`$ is a parabola whose maximum is achieved at $`l=2/h`$, supporting the result $`l_h^{}=2/h+1`$ described above.

For a correct calculation of the critical length we must analyze in detail the mechanisms of growth and shrinkage of a rectangular droplet; these are in general similar to those for Glauber dynamics , although the details are different. We will comment below on the possibilities of making the following heuristic discussion rigorous.

Consider first growth. From Figure 1 it is clear that a single plus protuberance on one of the four sides of the rectangular is not stable; growth proceeds through the formation of a double protuberance, which then grows “quickly” (i.e., on the time scale $`\mathrm{exp}2\beta (1h)`$) to complete the additional side. The parallel dynamics permits the double protuberance to form in one time step, as shown in Figure 2(a); the typical time for this process is

$$\tau _{\mathrm{double}}\mathrm{e}^{4\beta (3h)}.$$
(16)   
Alternatively, the protuberance can grow in two consecutive time steps, as shown in Figure 2(b); the parallel character of the dynamics enters here as well, since after formation of a single protuberance at the first step there must occur, at the second step, both the persistence of this protuberance, with probability $`\mathrm{exp}(2\beta (1h))`$, and the flip of a minus spin adjacent to the protuberance, also with probability $`\mathrm{exp}(2\beta (1h))`$ (see Figure 1); the typical time for this growth process is thus

$$\tau _{\mathrm{single}}\underset{\mathrm{first}\mathrm{step}}{\underset{}{\mathrm{e}^{2\beta (3h)}}}\times \underset{\mathrm{second}\mathrm{step}}{\underset{}{\mathrm{e}^{4\beta (1h)}}}$$
(17)   
Clearly $`\tau _{\mathrm{single}}\tau _{\mathrm{double}}`$ for $`\beta `$ large and hence the most efficient growth mechanism is the two-step one.

Again from Figure 1 it is clear that the most efficient shrinking mechanism is the usual corner erosion, shown in Figure 2(c); the shrinking is performed via a sequence of stable configurations. We estimate the time needed for the loss of one of the shorter sides of the rectangle, which requires the erosion of $`l1`$ sites (after which the remaining single protuberance vanishes rapidly). When $`\beta `$ is large, such a process will typically occur without backtracking. The rate at which the entire process occurs is thus estimated as the rate for one erosion, $`\mathrm{exp}^{2\beta (1+h)}`$, times the probability that $`l2`$ further erosions occur within the lifetime $`\mathrm{exp}^{2\beta (1h)}`$ of a stable configuration, which is of order $`\left[\mathrm{exp}^{2\beta (1+h)}\mathrm{exp}^{2\beta (1h)}\right]^{l2}`$. Thus the shrinking time is estimated as

$$\tau _{\mathrm{shrink}}\mathrm{e}^{2\beta (1+h)}\times \left[\frac{\mathrm{e}^{2\beta (1+h)}}{\mathrm{e}^{2\beta (1h)}}\right]^{l2}.$$
(18)   
The estimate of the shrinking time $`\tau _{\mathrm{shrink}}`$ can be supported by considering a random walk which models what happens on one edge of the droplet. Consider a Markov chain $`X_t`$, $`t=0,1,2,\mathrm{}`$, taking values in the nonnegative integers and with transition probabilities

$$P(k,l)=\{\begin{array}{cc}\frac{1}{2}\mathrm{e}^{\beta c}& \mathrm{if}l=k+1,\hfill \\ & \\ \frac{1}{2}\mathrm{e}^{\beta b}& \mathrm{if}l=k1,\hfill \\ & \\ 1\frac{1}{2}\mathrm{e}^{\beta c}\frac{1}{2}\mathrm{e}^{\beta b}& \mathrm{if}l=k,\hfill \\ & \\ 0& \mathrm{otherwise},\hfill \end{array}$$

and

$$P(0,1)=\frac{1}{2}\mathrm{e}^{\beta c},P(0,0)=1\frac{1}{2}\mathrm{e}^{\beta c},$$

where $`k1`$, $`c>b>0`$, and $`\beta >0`$. This chain, with $`b=2(1h)`$ and $`c=2(1+h)`$, is an approximate description of the behavior of the edge of a droplet (see Figure 3), if one thinks of $`X_t`$ as representing the number of minus spins on this edge at time $`t`$. In order to estimate $`\tau _{\mathrm{shrink}}`$ one should calculate the typical time to see $`l1`$ minus spins on the edge, starting from zero minus spins. Such an estimate, which agrees with (18), is provided by the following lemma, whose proof is parallel to that of Lemma 1 of .

###### Lemma 4.1

For $`k1`$, define the hitting time $`\tau _k^0`$ for the Markov chain $`X_t`$ with $`X_0=0`$ by

$$\tau _k^0\stackrel{\mathrm{def}}{=}\{t1:X_t=k\}.$$
(19)   
Then for any $`\epsilon >0`$,

$$P\left(\mathrm{e}^{\beta ck\beta b(k1)\beta \epsilon }<\tau _k^0<\mathrm{e}^{\beta ck\beta b(k1)+\beta \epsilon }\right)\stackrel{\beta \mathrm{}}{}1.$$
(20)   
To complete the derivation of the critical length for rectangular droplets we compare (17) and (18): growth occurs with probability one in the zero temperature limit, that is, $`ll_h^{}`$, if $`lim_\beta \mathrm{}\tau _{\mathrm{single}}/\tau _{\mathrm{shrink}}=0`$. This again leads to

$$l_h^{}=\frac{2}{h}+1.$$
(21)   
We believe that the above argument could be made rigorous along the lines of the corresponding arguments in Section 2 of . The main complicating factor appears to be that, because the growth time $`\tau _{\mathrm{single}}`$ is so large, processes beyond simple corner erosion must be accounted for when evaluating the shrinking time. For example, several corners may disappear in one step; more complicated processes, such as the two-step shrinkage by two sites shown in Figure 4, are also relevant. These modifications appear to be technical only and should not affect the estimate (18) of $`\tau _{\mathrm{shrink}}`$, but we have not carried out a complete analysis.

## 5 Monte Carlo results

The critical length $`l_h^{}`$ introduced in the previous section characterizes the behavior of the system in the limit $`\beta \mathrm{}`$. In this section we define a critical length $`l_{\beta ,h}^{}`$ at finite $`\beta `$ and describe the results of Monte Carlo simulations evaluating $`l_{\beta ,h}^{}`$ numerically for several values of $`\beta `$ and $`h`$. We find that when $`\beta `$ is large enough the resulting estimates of $`l_{\beta ,h}^{}`$ are close the theoretical estimate of $`l_h^{}`$ given in the previous section.

Let $`p_{\beta ,h}(l)`$ denote the probability that a square droplet of side $`l`$ grows and covers the whole lattice, that is, that in evolving from this initial configuration the system reaches the state $`+\underset{¯}{1}`$ before the state $`\underset{¯}{1}`$. Clearly $`p_{\beta ,h}(l)`$ is a nondecreasing function of $`l`$ with $`p_{\beta ,h}(0)=0`$ and $`p_{\beta ,h}(L)=1`$, so that the differences

$$d_{\beta ,h}(l)=p_{\beta ,h}(l)p_{\beta ,h}(l1)$$
(22)   
form a normalized probability distribution. In the limit $`\beta \mathrm{}`$ (assuming for simplicity that $`2/h`$ is not an integer), $`p_{\beta ,h}`$ reduces to a step function,

$$p_{\mathrm{},h}(l)=\{\begin{array}{cc}1,& \mathrm{if}ll_h^{},\hfill \\ 0,& \mathrm{if}l<l_h^{},\hfill \end{array}$$
(23)   
and $`d_{\beta ,h}(l)`$ to a unit mass on the critical length $`l_h^{}`$. At finite temperature, then, we define the critical length to be the mean of the distribution $`d_{\beta ,h}`$:

$$l_{\beta ,h}^{}=\underset{l}{}ld_{\beta ,h}(l).$$
(24)   
From (23) it follows that $`l_{\mathrm{},h}^{}=l_h^{}`$. We note that this approach in the numerical estimate of the critical length is different from that used in .

In the above discussion we have suppressed the dependence of $`p_{\beta ,h}`$, $`d_{\beta ,h}`$, and $`l_{\beta ,h}^{}`$ on the lattice size $`L`$, since for large $`\beta `$ and $`L`$ we expect $`l_{\beta ,h}^{}`$ as defined by (24) to be essentially independent of $`L`$.

We have carried out numerical experiments to estimate the function $`p_{\beta ,h}(l)`$, and hence $`l_{\beta ,h}^{}`$, for $`\beta =0.9,1.1,1.3`$ and $`h=0.05,0.1,0.2`$; we varied $`l`$ over the range of values in which $`p_{\beta ,h}(l)`$ changes rapidly and made $`N=100`$ runs for each value of $`l`$. For each run, we first prepared our system in a starting configuration characterized by a single square droplet of plus spins of size $`l`$, placed in a lattice whose size $`L`$ was chosen large enough to avoid boundary effects. We then followed the evolution of the system and decided, by means of lower and upper cutoffs on the total system magnetization, whether the droplet would ultimately grow or shrink. We also introduced a cutoff on the total length of each run, chosen as a function of $`\beta `$ so that for most runs the fate of the droplet was determined before the cutoff was reached; in the case $`\beta =1.3`$, the highest value of $`\beta `$ we have considered, this cutoff was $`200,000`$ iterations. Letting $`G_{\beta ,h}(l)`$ denote the number of times that the droplet grew and $`S_{\beta ,h}(l)`$the number of times that it shrank, and assuming that the fraction of the remaining $`NG_{\beta ,h}S_{\beta ,h}`$ cases (in which we did not determine the behavior) in which the droplet would have grown if we had waited long enough is the same as for the cases in which the behavior was determined, we are led to the estimate

$$p_{\beta ,h}(l)=\frac{G_{\beta ,h}(l)}{G_{\beta ,h}(l)+S_{\beta ,h}(l)}.$$
(25)   
From the estimated values of $`p_{\beta ,h}(l)`$ we computed the finite temperature critical length, via (24), and the standard deviation of the distribution $`d_{\beta ,h}`$; the values are recorded in Table 1. The results are in a very good agreement with our theoretical prediction: when the temperature is lowered, the numerical measure of the critical length tends to the zero temperature theoretical prediction $`l_h^{}=2/h+1`$. We did not consider higher values of $`\beta `$ because too long runs would have been needed, but the values we have considered seem to be sufficient to see the zero temperature limit behavior. Note that the standard deviation of the distribution $`d_{\beta ,h}(l)`$ decreases when $`\beta `$ is increased. This good behavior is clearest in the case of small external magnetic field; presumably, higher values of $`\beta `$ should be considered at higher $`h`$ to approach the limiting behavior.

We observed that the typical time for growth of the initial square depended strongly on $`\beta `$ and $`h`$, but not on $`l`$; while the typical shrinking time increase sensibly when $`l`$ is increased. This is qualitatively in agreement with theoretical estimates (17) and (18).

## 6 Conclusions

In this paper we have studied the problem of metastable states in probabilistic cellular automata, viewing the latter as the simplest instance of models evolving under parallel dynamics. All detailed work has been focussed on a particular case: a two dimensional model on the square lattice in which the probability of a spin flip at site $`x`$ depends only on the total magnetization of the set of five spins in a cross centered at $`x`$ (see (1)).

We conclude that the general pattern of analysis which has been used for similar models evolving under Glauber dynamics applies here as well, since events in which the system makes a one-step transition to a configuration significantly different from the current one can be neglected in the low-temperature limit. In particular, we argue that the path of escape from metastability passes through a critical rectangular droplet. On the other hand, the parallel nature of the dynamics does influence the details of the analysis of the escape time and, in particular, adds enough complications to make a rigorous analysis more difficult than in the Glauber case.

For the model in question we have shown, through heuristic arguments and Monte Carlo simulations, that the critical length of a rectangular droplet is $`l_h^{}=2/h+1`$. Our theoretical prediction is valid only in the limit of zero temperature, but our simulations confirm estimates close to the theoretical ones even at finite temperature.

It is natural to ask whether escape from the metastable state is facilitated or hindered by the use of parallel (as opposed to serial) dynamics. It is not clear that this question has a universal answer, but as a preliminary approach we may ask what would happen in the model of this paper if a serial evolution rule were adopted, so that at each time step one spin is chosen at random, with uniform probability, and then updated with probability given by (1).

As remarked in Section 4, we used the parallel character of the dynamics only in the estimate of $`\tau _{\mathrm{single}}`$, so that to estimate of the critical length in the serial case one should compare the shrinking time (18) with a new growing time $`\tau _{\mathrm{growth}}\mathrm{exp}[2\beta (3h)+2\beta (1h)]`$, obtained by noting that in the second step of the double protuberance growth the persistence probability of the single protuberance need be taken into account. Comparison of these two times shows that the critical length in the serial case is given, for $`h`$ very small, by $`3/2h+1`$. Thus, for these models, the parallel rule leads to a larger critical droplet and a slower exit from the metastable phase.

Acknowledgments

One of the authors (E.C.) wishes to express his thanks to the Mathematics Department of Rutgers University for its very kind hospitality and to Enzo Olivieri for useful discussions. E.C. also thanks Istituto Nazionale di Fisica Nucleare - Sezione di Bari and Dipartimento di Fisica dell’Università degli Studi di Bari for their financial support. Work at Rutgers was supported by NSF Grant DMR 95 - 23266. JLL would also like to thank DIMACS and its supporting agencies the NSF under contract STC–91–19999 and the NJ Commission on Science and Technology. We thank Alex Mazel for very useful discussions.