Generic phase coexistence in the totally asymmetric kinetic Ising model

TL;DR

This paper investigates phase coexistence in a specific asymmetric kinetic Ising model, combining analytical and numerical methods to understand droplet dynamics and phase boundaries near critical points.

Contribution

It provides a detailed analysis of phase coexistence and droplet dynamics in a nonequilibrium kinetic Ising model, including exact and numerical results on interface behavior.

Findings

Exact shape of ballistic interfaces at zero temperature via TASEP mapping

Temperature-dependent droplet shrinking velocity and critical exponents

Phase boundary fields vanish with specific critical exponents near T_c

Abstract

The physical analysis of generic phase coexistence in the North-East-Center Toom model was originally given by Bennett and Grinstein. The gist of their argument relies on the dynamics of interfaces and droplets. We revisit the same question for a specific totally asymmetric kinetic Ising model on the square lattice. This nonequilibrium model possesses the remarkable property that its stationary-state measure in the absence of a magnetic field coincides with that of the usual ferromagnetic Ising model. We use both analytical arguments and numerical simulations in order to make progress in the quantitative understanding of the phenomenon of generic phase coexistence. At zero temperature a mapping onto the TASEP allows an exact determination of the time-dependent shape of the ballistic interface sweeping a large square minority droplet of up or down spins. At finite temperature, measuring the mean lifetime of such a droplet allows an accurate measurement of its shrinking velocity $v$, which depends on temperature $T$ and magnetic field $h$. In the absence of a magnetic field, $v$ vanishes with an exponent $\Delta_v\approx2.5\pm0.2$ as the critical temperature $T_c$ is approached. At fixed temperature in the ordered phase, $v$ vanishes at the phase-boundary fields $\pm h_{\rm b}(T)$ which mark the limits of the coexistence region. The latter fields vanish with an exponent $\Delta_h\approx3.2\pm0.3$ as $T_c$ is approached.

Tables

Table 1. Table 1: List of spin-flip moves for the NEC Toom model, with the convention ( σ N , σ , σ E subscript 𝜎 N 𝜎 subscript 𝜎 E {\sigma}_{\rm N},{\sigma},{\sigma}_{\rm E} ) for local spin configurations, and expressions of the corresponding transition probabilities (discrete time) or rates (continuous time).

spin-flip move	prob./rate
$+++\to +-+$	$p$
$++-\to +--$	$p$
$-++\to --+$	$p$
$-+-\to ---$	$1-q$
$+-+\to +++$	$1-p$
$+--\to ++-$	$q$
$--+\to -++$	$q$
$---\to -+-$	$q$

Table 2. Table 2: List of spin-flip moves for the totally asymmetric kinetic Ising model, with the convention ( σ N , σ , σ E subscript 𝜎 N 𝜎 subscript 𝜎 E {\sigma}_{\rm N},{\sigma},{\sigma}_{\rm E} ) for local spin configurations, and expressions of the corresponding rates ( 1.7 ).

spin-flip move	rate
$+++\to +-+$	$\frac{\alpha }{2}{\left(1-\gamma \right)}^2\left(1-\kappa \right)$
$++-\to +--$	$\frac{\alpha }{2}\left(1-{\gamma }^2\right)\left(1-\kappa \right)$
$-++\to --+$	$\frac{\alpha }{2}\left(1-{\gamma }^2\right)\left(1-\kappa \right)$
$-+-\to ---$	$\frac{\alpha }{2}{\left(1+\gamma \right)}^2\left(1-\kappa \right)$
$+-+\to +++$	$\frac{\alpha }{2}{\left(1+\gamma \right)}^2\left(1+\kappa \right)$
$+--\to ++-$	$\frac{\alpha }{2}\left(1-{\gamma }^2\right)\left(1+\kappa \right)$
$--+\to -++$	$\frac{\alpha }{2}\left(1-{\gamma }^2\right)\left(1+\kappa \right)$
$---\to -+-$	$\frac{\alpha }{2}{\left(1-\gamma \right)}^2\left(1+\kappa \right)$

Equations

$$w=\frac{1}{2}\left(b-c\sigma ({\sigma }_{\mathrm{N}}+{\sigma }_{\mathrm{E}})+c{\sigma }_{\mathrm{N}}{\sigma }_{\mathrm{E}}-h\sigma \right),$$

$$\mathcal{H}=-J\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j,$$

$$w$$

$$\gamma =\tanh 2K,$$

$$w={\text{rme}}^{-K\sigma ({\sigma }_{\mathrm{N}}+{\sigma }_{\mathrm{E}})},$$

$$\mathcal{H}=-J\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j-h\sum _i{\sigma }_i.$$

$$w=\frac{\alpha }{2}(1-\gamma \sigma {\sigma }_{\mathrm{N}})(1-\gamma \sigma {\sigma }_{\mathrm{E}})(1-\kappa \sigma ),$$

$$\kappa =\tanh (h/T).$$

$$T_c=\frac{2}{\ln (1+\sqrt{2})}\approx \mathrm{2.269.}$$

$$w=\frac{\alpha}{2}\big{(}1-{\sigma}\tanh[K({\sigma}_{\rm N}+{\sigma}_{\rm E}+{\sigma}_{\rm W}+{\sigma}_{\rm S})]\big{)},$$

$$w=\frac{\alpha }{2}\left(1-\frac{\gamma }{2}\sigma ({\sigma }_{\mathrm{N}}+{\sigma }_{\mathrm{E}})\right).$$

$$w=\frac{\alpha }{2}(1-\sigma {\sigma }_{\mathrm{N}})(1-\sigma {\sigma }_{\mathrm{E}})(1-\kappa \sigma ).$$

$$w=2\alpha (1-\kappa );$$

$$w=2\alpha (1+\kappa ).$$

$$10\to 01,$$

$$w=2\alpha (1\pm \kappa ),$$

$$\rho (\theta )=\frac{\cos \theta }{\sin \theta +\cos \theta }.$$

$$v(\theta )=w\frac{\sin \theta \cos \theta }{\sin \theta +\cos \theta }$$

$$v=\frac{w}{2\sqrt{2}},$$

$$t^{\star }\approx \frac{N\sqrt{2}}{v},$$

$$t^{\star }\approx \frac{4N}{w}.$$

$$v_0=\frac{\alpha }{\sqrt{2}},t_0^{\star }\approx \frac{2N}{\alpha }.$$

$$v=v_0(1\pm \kappa ),t^{\star }=\frac{t_0^{\star }}{1\pm \kappa }.$$

$$wt^{\star}\approx\bigl{(}\sqrt{M}+\sqrt{N}\bigr{)}^{2}=M+N+2\sqrt{MN}$$

$$\underset{N}{\underbrace{1\dots 1}}\underset{N}{\underbrace{0\dots 0}}⟶\underset{N}{\underbrace{0\dots 0}}\underset{N}{\underbrace{1\dots 1}}.$$

$$\frac{\partial \rho }{\partial t}+(1-2\rho )\frac{\partial \rho }{\partial x}=0.$$

$$\rho(x,0)=\left\{\matrix{1\quad&(0<x<N),\hfill\cr 0\hfill&(N<x<2N).\hfill}\right.$$

$$X={\int }_0^x\rho (x^{\prime },t)\mathrm{d}{x}^{\prime },Y=N-x+X.$$

$$\rho(x,t)=\left\{\matrix{1\hfill&(0<x<N-t),\hfill\cr(N+t-x)/(2t)\quad&(N-t<x<N+t),\hfill\cr 0\hfill&(N+t<x<2N)\hfill}\right.$$

$$(X-Y{)}^2+2t(X+Y-2N)+t^2=0,$$

Full text

1 Introduction and summary

As was demonstrated by Toom [1], an irreversible probabilistic cellular automaton can exhibit generic phase coexistence, i.e., two stable stationary phases can coexist in a whole domain of parameter space. This property is alternatively referred to as generic bistability or generic nonergodicity. Generic phase coexistence in nonequilibrium systems stands in contrast with the case of equilibrium systems, where coexistence requires the equality of the thermodynamical free energies of the two phases, in agreement with Gibbs phase rule, and therefore some fine tuning of parameters. Taking the example of the two-dimensional ferromagnetic Ising model in its low-temperature phase, the two ordered phases only coexist at zero magnetic field.

A prototypical instance of a nonequilibrium model exhibiting generic bistability is the two-dimensional North-East-Center (NEC) Toom model, originally analysed in simple physical terms by Bennett and Grinstein [2], then further investigated in [3, 4, 5]. This model consists of a square lattice of Ising spins ${\sigma}=\pm 1$, evolving according to a synchronous discrete-time dynamics, where the spin ${\sigma}$ located at a generic site, named the central spin, is only influenced by its North and East neighbours ${\sigma}_{\rm N}$ and ${\sigma}_{\rm E}$. At each time step all spins are updated synchronously, according to the following rules:

• •

First, a deterministic majority rule is applied: if two or more of the spins ${\sigma}$, ${\sigma}_{\rm N}$ and ${\sigma}_{\rm E}$ point up (down), then the value ${\sigma}$ of the central spin at time $t+1$ is set to $+1$ ($-1$).

• •

Then, a probabilistic rule is applied: if ${\sigma}(t+1)=+1$, then it is flipped to $-1$ with probability $p$, while if ${\sigma}(t+1)=-1$, it is flipped to $+1$ with probability $q$.

The corresponding transition probabilities are given in table 1. The noise parameter $r=p+q$ is analogous to a temperature, while the symmetry-breaking parameter $h=q-p$ is analogous to a magnetic field. Indeed, when $q>p$, up spins are favored over down spins. The remarkable property of this model is that its stationary-state phase diagram shows coexistence in a whole region of parameter space. The two stable phases have magnetizations with opposite signs. As discussed in [2, 3, 4, 5], the same holds true if, instead of being updated synchronously, the spins are updated in a random sequential manner. The rate at which the central spin is flipped reads

[TABLE]

with $b=(1+r)/2$ and $c=(1-r)/2$. The values given in table 1 are now interpreted as rates.

The gist of the phenomenon of generic bistability relies on the dynamics of interfaces. Let us summarize the argument put forward originally in [2], then further expanded in [3, 4, 5]. Consider a ferromagnetic system of Ising spins in its low-temperature phase in the absence of a magnetic field. If the microscopic dynamical rules are reversible, a flat interface will not move on average, and a curved one will move only under the effect of surface tension. For dynamical rules such as those of the NEC Toom model, which owe their irreversibility to a spatial anisotropy, a flat interface generically moves ballistically with a velocity depending on the orientation of the interface. In the presence of a weak magnetic field $h$, this velocity will be corrected by a term proportional to $h$, and therefore remains nonzero in a whole domain of parameter space. This picture concerning flat interfaces can now be completed by an argument concerning droplet dynamics. If an ordered phase contains a large droplet of the other phase, this droplet will shrink ballistically, i.e., in a time proportional to its diameter, as a consequence of the orientation dependence of the interface velocity. So, such a dynamics is able to spontaneously heal defects in both ordered phases, with large defects shrinking on a ballistic time scale. The simultaneous stability of the two ordered phases takes place in a whole coexistence region, where both noise (measured by the temperature $T$) and bias (measured by the magnetic field $h$) are small enough. This coexistence region is shown schematically in figure 1. It is delimited by curves of the form $h=\pm{h_{\rm b}}(T)$, where the phase-boundary field vanishes according to the power law ${h_{\rm b}}(T)\sim(T_{c}-T)^{\Delta_{h}}$ as the critical temperature is approached.

In addition to the NEC Toom model, other models exhibiting generic bistability were subsequently analysed by several groups, mostly in the two-dimensional setting [6, 7, 8, 9, 10, 11, 12], but also in three [13] and even higher [14] dimensions. Besides its major conceptual role recalled above, the ballistic shrinking of large droplets (or islands) of prescribed form has been instrumental in the numerical determination of the boundaries of the nonergodic region performed in [2, 4, 6, 12]. Let us mention a review putting the phenomenon of generic coexistence in a very broad perspective [15], as well as a very recent work [16] pinpointing various subtleties related to spatial aspects of phase coexistence in nonequilibrium systems.

In the present work we revisit the investigations made in [2, 3, 4, 6], with the aim of giving a more quantitative analysis of generic bistability in two-dimensional spin systems, putting a special emphasis both on the zero-temperature limit and on the vicinity of the critical point. To this end we focus on a special model in the class of totally asymmetric NEC dynamics. This model, introduced initially by Künsch [17] in a slightly different form (see (1.5) below), possesses the remarkable property [17, 18, 19] of having the same stationary measure as the two-dimensional Ising model at temperature $T$ in the absence of a magnetic field, with the usual nearest-neighbour ferromagnetic Hamiltonian

[TABLE]

where the sum runs over pairs $\langle i,j\rangle$ of neighbouring spins. The rate defining the model reads

[TABLE]

where $\alpha$ is an arbitrary frequency fixing the time scale, ${\sigma}$ is the flipping spin, ${\sigma}_{\rm N}$ and ${\sigma}_{\rm E}$ are the North and East neighbouring spins, and

[TABLE]

where $K=J/T$ is the reduced coupling constant. This model is the totally asymmetric version of the two-dimensional directed Ising model, investigated in [18, 20, 21], whose definition encompasses the case where the dynamics is partially asymmetric, i.e., where the South and West spins keep an influence on the central spin. A number of physical properties of this model have been described in these references. In the one-dimensional situation too, the general directed Ising model has rich dynamical properties [22, 23].

The expression (1.3) is invariant under up-down spin symmetry and satisfies the condition of global balance, that is to say, leads to a Gibbsian stationary measure with respect to the Hamiltonian (1.2). The dynamics thus defined is irreversible, since the condition of detailed balance is not satisfied [17, 18, 19]. This model is actually the unique Gibbsian two-dimensional spin system with totally asymmetric NEC dynamics. Imposing the condition of global balance with respect to the Hamiltonian (1.2) on a NEC dynamics indeed uniquely determines the rate (1.3) [18, 19]. Fixing the time scale in (1.3) by the choice $\alpha=2\cosh^{2}2K$ yields the exponential form $w=\exp(-2K{\sigma}({\sigma}_{\rm N}+{\sigma}_{\rm E}))$. Note that the expression of the rate originally given by Künsch [17],

[TABLE]

leads to a ferromagnetic Gibbs measure at temperature $T/2$.

At finite temperature and in the presence of an external magnetic field $h$, it can be shown, along the lines of [19], that there is no NEC dynamics yielding global balance with respect to the full Hamiltonian

[TABLE]

In this work we choose to introduce the effect of the magnetic field $h$ in a multiplicative form, as originally done by Glauber for the one-dimensional kinetic Ising model [24]. We shall thus consider the dynamics defined by the rate function

[TABLE]

with

[TABLE]

The values taken by this rate function for the eight possible local configurations $({\sigma}_{N},{\sigma},{\sigma}_{E})$ are given in table 2. As analysed in detail in the present work, this model exhibits generic phase coexistence. Though its stationary measure is not known in closed form for generic values of the parameters $\gamma$ and $\kappa$, it is however known exactly in two special cases. In the absence of a magnetic field ($\kappa=0$), as said above, the stationary measure is the same as that of the equilibrium model with Hamiltonian (1.2), which is (1.6) restricted to its first (ferromagnetic) term. In the absence of interactions ($\gamma=0$), the dynamics obeys detailed balance with respect to the second (paramagnetic) term of the Hamiltonian (1.6). These special cases demonstrate that temperature $T$ and magnetic field $h$ are natural and properly normalized extensions of what they are in the usual Ising model.

The main motivation behind our choice of the rate (1.7) is that it is more advantageous, especially from the standpoint of numerical simulations, to study generic bistability in a model whose stationary properties are known exactly in the absence of a magnetic field. In particular, the critical temperature is that of the two-dimensional ferromagnetic Ising model, given by $\sinh 2K_{c}=1$. Henceforth we choose energy units such that $J=1$, and so $K=1/T$. We have thus in particular

[TABLE]

Coming back to the NEC Toom model (1.1), let us point out that this model is just a generic member of the class of NEC dynamics, with no specific property. In particular its stationary-state measure is not known in closed form, even in the absence of a bias ($p=q$). Let us note however that the rate (1.1) with $p=q=0$ is the same as (1.7) with $\gamma=1$, $h=0$ and $\alpha=1/2$. This can also be checked by comparing tables 1 and 2.

Then, a natural question to ask is: what are the minimal conditions required to produce generic nonergodicity? It was argued in [4] that the only indispensable ingredient seems to be irreversibility in the specific form of spatially asymmetric rules. Restricting ourselves to spin models, this statement is confirmed by the behaviour of the NEC Toom model, of the model studied in [6], and of the directed kinetic Ising model investigated in the present work. There are however exceptions to this picture. One notable case is the following one. Starting with the usual rate function

[TABLE]

which is the direct generalisation to two dimensions of the Glauber rate function [24], and just suppressing the influence of the West and South spins, results in the following expression, after reduction on a basis of spin operators:

[TABLE]

This rate function describes a truncated noisy voter model with no West and South influence [18, 19]. Just as the isotropic two-dimensional voter model, this model has no phase transition [18, 25].

We now give an overview and a discussion of the main results obtained in the present work. In section 2 we study the zero-temperature limit of the model, allowing the presence of a vanishingly small magnetic field, proportional to temperature. The only flippable spins are located in the NE corners of interfaces separating domains of up and down spins (see figure 2). As a consequence, a descending interface moves ballistically, while an ascending one is immobile. A mapping onto the totally asymmetric exclusion process (TASEP) allows to determine the angle-dependent velocity $v(\theta)$ of a flat interface. This velocity takes its maximal value $v$ when the interface is normal to the bias of the dynamics, which is oriented from NE to SW. In the absence of a magnetic field, we have $v=v_{0}=\alpha/\sqrt{2}$ (section 2.1). We then focus our attention on the fate of a large square droplet of linear size $N$, embedded in a sea of opposite spins. This droplet is swept by a ballistic interface propagating at velocity $v$ in the direction of the bias. Accordingly, its lifetime $t^{\star}$ grows proportionally to $N$ (section 2.2). This interface is however not straight, even if the droplet is very large. Its time-dependent shape is obtained by means of the hydrodynamic approach to the TASEP, based on the inviscid Burgers equation (section 2.3). It is given by an arc of parabola (see figure 5). We finally consider the fluctuations of the droplet’s lifetime $t^{\star}$ around the ballistic law. According to the KPZ theory, these fluctuations scale as $N^{1/3}$ times a dimensionless random variable $\chi$ with a non-trivial distribution (section 2.4). Numerical simulations strongly corroborate the ansatz that $\chi$ is proportional to the variable $\xi_{\rm GUE}$ distributed according to the Tracy-Widom $F_{2}$ law (see figure 7).

In section 3 we investigate the model at finite temperature and zero field by means of numerical simulations. The dynamics of a large droplet concomitantly exhibits two distinct features: its center of mass is advected ballistically in the direction of the bias with velocity $v_{\rm cm}$, whereas the droplet shrinks ballistically with velocity $v$, and therefore disappears after a time $t^{\star}$ which grows linearly with $N$. As temperature increases, $v_{\rm cm}$ increases steadily, while the droplet becomes progressively rounder, and its lifetime $t^{\star}$ increases very fast as $T_{c}$ is approached (section 3.1). The quantitative determination of the shrinking velocity $v$ relies on a precise definition of $t^{\star}$. In analogy with the zero-temperature situation, where $t^{\star}$ is the time where the system reaches its ground state, in the finite-temperature situation we define the stopping time $t^{\star}$ as the first time where the mean magnetization exceeds the equilibrium spontaneous magnetization (section 3.2). The outcome of this study is shown in figure 12, depicting the ratio $v/v_{0}$ as a function of $T/T_{c}$. At low temperature, this ratio departs from unity as $\rme^{-2/T}$. At criticality, a finite-size scaling analysis yields $v\sim(T_{c}-T)^{\Delta_{v}}$, with ${\Delta_{v}}\approx 2.5\pm 0.2$ (section 3.3).

Section 4 is devoted to a numerical investigation of the generic situation where the temperature $T$ and the magnetic field $h$ are both non-zero. The system exhibits generic bistability in a whole coexistence region of parameter space, delimited by the phase-boundary fields $\pm{h_{\rm b}}(T)$. For small enough values of the magnetic field, the situation is similar to what happens for $h=0$, i.e., a large droplet of down spins disappears in a ballistic time. The corresponding shrinking velocity $v$ depends on both $T$ and $h$, vanishing as $h\to-{h_{\rm b}}(T)$, where the up phase loses its stability (section 4.1). An efficient numerical determination of the phase-boundary field again requires an operational definition of a stopping time $t^{\star}$. We have chosen to define it as the first time where the magnetization of a system with droplet initial condition exceeds that of an initially homogeneous system. We have also measured the splitting probability $p^{\star}$ of the up phase, i.e., the probability for the system with droplet initial condition to end in the up phase (section 4.2). The finite-size scaling analysis of $p^{\star}$ yields the most accurate determination of the phase-boundary field. The outcome of this study is shown in figure 17, depicting ${h_{\rm b}}(T)$ as a function of $T/T_{c}$. This figure therefore provides a quantitative version of the phase diagram shown schematically in figure 1. At low temperature we have ${h_{\rm b}}\approx 2-BT$, with $B\approx 1.63$, whereas a critical finite-size scaling analysis yields ${h_{\rm b}}\sim(T_{c}-T)^{\Delta_{h}}$, with ${\Delta_{h}}\approx 3.2\pm 0.3$. Our predictions for the critical exponents ${\Delta_{v}}$ and ${\Delta_{h}}$ are put in perspective with existing analytical and numerical results. We obtain in particular the estimate $y_{i}\approx-1.3\pm 0.2$ for the irrelevant scaling dimension entering the heuristic approach of [6] (section 4.3). Finally, the magnetization curves of the up phase are shown in figure 20 for three temperatures in the ordered phase. These curves are observed to be very regular near the phase-boundary field (section 4.4).

2 Zero-temperature dynamics

We start by considering the zero-temperature limit of the dynamics ($K=1/T\to\infty$, and so $\gamma\to 1$). The parameter $\kappa=\tanh(h/T)$ is kept fixed, thus allowing to consider vanishingly small magnetic fields $h$, proportional to temperature $T$.

For $\gamma=1$ and $\kappa$ arbitrary, the rate function (1.7) reads

[TABLE]

It takes non-zero values for only two configurations, corresponding to lines 4 and 5 of table 2:

• •

If ${\sigma}_{\rm N}={\sigma}_{\rm E}=-1$, ${\sigma}$ is flipped from $+1$ to $-1$ with rate

[TABLE]

• •

If ${\sigma}_{\rm N}={\sigma}_{\rm E}=+1$, ${\sigma}$ is flipped from $-1$ to $+1$ with rate

[TABLE]

The only flippable spins are those located in the NE corners of interfaces between domains of up and down spins (see figure 2). As a consequence, a descending interface is mobile, while an ascending one is not.

2.1 Mapping onto the TASEP

The dynamics of a descending interface can be exactly mapped onto that of the totally asymmetric exclusion process (TASEP) [26, 27]. To this end, horizontal segments of the interface are identified with particles (denoted by 1) and vertical ones with holes (denoted by 0). Flipping a spin corresponds to the displacement of a particle to the right, according to

[TABLE]

with rate

[TABLE]

where the plus (minus) sign applies if up spins are above (below) the interface.

The mapping onto the TASEP provides a useful tool to investigate the dynamics. For the time being, we consider a single macroscopically flat descending interface, making an angle $\theta$ with the horizontal ($0<\theta<\pi/2$). The particle density of the associated TASEP is equal to the fraction of horizontal segments, i.e.,

[TABLE]

The current across the TASEP in its stationary state reads $j(\theta)=w\rho(\theta)(1-\rho(\theta))$. Accordingly, the interface moves with a constant normal velocity $v(\theta)$, that can be derived as follows. A useful intermediate quantity is the time $\tau$ it takes for each particle to move on average one step forward, so that the integrated current over a time lapse $\tau$ equals the particle density. This yields $\tau=\rho(\theta)/j(\theta)$. During the same time, the interface moves by a distance $d=\cos\theta$ (measured perpendicularly to the direction of the interface) [21]. This yields $\tau=d/v(\theta)$. Putting everything together, we obtain the expression

[TABLE]

for the angle-dependent normal velocity of a flat descending interface. This velocity takes its maximal value,

[TABLE]

for $\theta=\pi/4$, i.e., when the interface is normal to the bias of the dynamics, and moves in the direction of this bias, i.e., from NE to SW.

2.2 Lifetime of a large square droplet

Along the lines of earlier works [2, 4, 6, 12, 21], we shall focus our attention on the fate of a single large droplet embedded in a sea of the opposite phase. For definiteness and simplicity, we shall consider square droplets. Other shapes could be considered as well (see [21]). In this section, devoted to an analytical study of zero-temperature dynamics, we consider a large square droplet of linear size $N$, initially consisting of up spins, immersed in an infinite sea of down spins, or vice versa. In sections 3 and 4, devoted to numerical simulations, the sea itself will be modelled as a larger square of size $N_{\rm sea}$, with periodic boundary conditions.

As soon as it is prepared and subjected to zero-temperature dynamics, the droplet is swept by an interface propagating at the maximal velocity $v$ given by (2.8) along the direction of the bias, from its NE corner to its SW corner. The distance between these two points being $N\sqrt{2}$, the lifetime $t^{\star}$ of a large droplet is given by the ballistic formula

[TABLE]

i.e.,

[TABLE]

In the absence of a magnetic field ($\kappa=0$) we have $w=2\alpha$, both for up and down droplets, and so the velocity (2.8) and the droplet lifetime (2.10), denoted by $v_{0}$ and $t^{\star}_{0}$, read

[TABLE]

In the presence of a magnetic field, we have

[TABLE]

Plus (minus) signs hold for droplets of down (up) spins, whose lifetime decreases (increases) as a function of the magnetic field. The lifetimes of the two types of droplets are nevertheless finite for all values of the reduced field $h/T$, i.e., for all parameters $\kappa$ in the range $-1<\kappa<1$.

The simplicity of the ballistic formula (2.10) however hides many interesting features of zero-temperature dynamics, to be described in the sequel:

• •

The ballistic interface is not straight, even if the droplet is very large. The time-dependent shape of the interface will be derived in section 2.3 for a large square droplet, using a hydrodynamic approach for the TASEP based on the Burgers equation.

• •

The hydrodynamic approach can be applied to large droplets with arbitrary initial shapes. Ballistic laws involving simple geometric factors however only hold for a class of convex droplet shapes which are symmetric with respect to the direction of the bias, such as those considered in [21]: square, circle, right isosceles triangle. (The droplet depicted schematically in figure 12 of this reference should be more correctly represented as convex and symmetric.) For more general shapes, the hydrodynamic approach yields more involved expressions for the lifetime $t^{\star}$. Skipping details, let us mention that it yields the prediction

[TABLE]

for the lifetime of a large rectangular droplet of size $M\times N$.

• •

The lifetime $t^{\star}$ of an individual droplet is a random quantity which fluctuates around the ballistic law (2.10). The size of these fluctuations scales as $N^{1/3}$, according to the KPZ theory (see section 2.4).

2.3 Shape of the interface sweeping a large square droplet

In this section we give an analytical description of the time-dependent shape of the interface sweeping a large square droplet of size $N$, using the exact mapping of the zero-temperature dynamics onto the TASEP, as described above. Throughout sections 2.3 and 2.4, the time unit is fixed by setting $w=1$.

A square droplet of size $N$ corresponds to the TASEP on a closed interval of length $2N$, whose initial configuration consists of $N$ particles followed by $N$ holes. The dynamics brings the system to the absorbing configuration where the particles have moved to the right of the holes:

[TABLE]

The lifetime $t^{\star}$ of the droplet is the random duration of this process.

On the ballistic scale, i.e., when spatial distance and time are large and proportional, the appropriate level of description is the hydrodynamic one, which consists in introducing the coarse-grained particle density $\rho(x,t)$, which obeys the inviscid Burgers equation [27, 28, 29, 30]

[TABLE]

The initial condition ($N$ particles followed by $N$ holes) corresponds to a step for the density:

[TABLE]

The time-dependent density profile $\rho(x,t)$ evolved from this initial condition yields a parametric representation of the interface sweeping the droplet in the spin model:

[TABLE]

The solution to the Burgers equation (2.15) with initial condition (2.16) involves two successive stages.

Stage I ($\,0<t<N$).

In this first stage, the density profile evolves as if the system were infinite. The system is thus equivalent to an infinite TASEP with a step initial condition, where all the sites to the left of an ‘origin’ (here, $x=N$) are occupied, while all the sites to its right are empty. In the language of the spin model, the initial configuration consists of an infinite quadrant of up spins surrounded by down spins, or vice versa. The corresponding solution to the Burgers equation is the so-called rarefaction (or fan) scaling solution [31]

[TABLE]

(see figure 3, left). According to (2.17), the corresponding interface is an arc of the parabola with equation

[TABLE]

or, equivalently,

[TABLE]

a result first derived by Rost [28]. The interface is tangent to the North and East sides of the droplet (see figure 3, right). Its apex moves at a constant velocity along the diagonal, according to

[TABLE]

At the end of the first stage $(t=N$), the density profile is linear over the whole system: $\rho(x,N)=1-x/(2N)$, while the interface touches the NW and SE corners of the droplet.

Stage II ($N<t<4N$).

In this second stage, the finite system evolves differently from the infinite one. Two shocks indeed invade the system from either end. This corresponds to the fact that the West ($X=0$) and South ($Y=0$) sides of the droplet begin to be eroded. We have

[TABLE]

(see figure 4, left). In the central region, the dynamics is the continuation of the previous stage. The width $2\lambda$ of that active region is ruled by the differential equation

[TABLE]

The expression on the right-hand side represents the velocity of the left shock, given by the value of $1-\rho$ to its immediate right, which can be read off from (2.22). An integration of (2.23) with initial condition $\lambda=N$ for $t=N$ yields

[TABLE]

The active region shrinks in the course of time and disappears at the end of the process, as $\lambda=0$ for $t=t^{\star}=4N$. Throughout Stage II, the shape of the interface is still given by the parabolic arc (2.20). Its apex still drifts ballistically according to (2.21), while its endpoints now sit on the West and South sides of the droplet (see figure 4, right).

As a summary, figure 5 shows the successive interface shapes at 16 equally spaced instants of time $t$ between 0 and $t^{\star}=4N$. The blue curve shows the interface at the end of the first stage $(t=N=t^{\star}/4$), when it touches the NW and SE corners of the droplet. Even though the interface is macroscopically curved at all times, its apex (red symbols) drifts along the diagonal with a constant velocity, according to (2.21). This explains the simplicity of the ballistic law (2.10), in spite of the non-triviality of the interface shape.

Near the end of the droplet’s life ($t^{\star}-t\ll t^{\star}$), its size is much smaller than its initial size, and it has the shape of a right isosceles triangle. This limiting triangular shape turns out to be universal, in the sense that it is observed (at least) for all initial shapes which are convex and symmetric with respect to the direction of the bias. This universality has been noticed in several works devoted to generic coexistence caused by spatial asymmetric dynamics [4, 6, 10], as well as in [21]. The right isosceles triangular shape is indeed singled out by the zero-temperature dynamics: it is the only droplet shape which shrinks in a self-similar way, i.e., keeping its shape, with its West and South sides being immobile, and its hypotenuse moving at constant velocity $v$. Universal properties of limiting domain shapes have also been discussed for other classes of kinetic Ising models [32, 33, 34, 35, 36].

2.4 Fluctuations around the ballistic law

We now consider the fluctuations of the lifetime of individual droplets around the ballistic law (2.10). As discussed in section 2.3, a square droplet of size $N$ corresponds to the TASEP on a closed interval of length $2N$, and its lifetime $t^{\star}$ is the duration of the process (2.14). This lifetime is a random quantity which depends on the whole stochastic history of the system. In order to describe its fluctuations, we have to go one step beyond the Burgers equation (2.15), and to consider the Kardar-Parisi-Zhang (KPZ) theory (see [37] for an overview). This approach predicts that the deterministic ballistic law (2.10) is complemented by a subleading fluctuating part of the form

[TABLE]

where $\chi$ is a dimensionless random variable with a non-trivial limiting distribution.

Numerical simulations (see below) corroborate the following ansatz for the variable $\chi$:

[TABLE]

where $\xi_{\rm GUE}$ is distributed according to the Tracy-Widom $F_{2}$ law. This distribution, originally derived for the largest eigenvalue of a large Hermitian matrix in the Gaussian Unitary Ensemble (GUE) [38], has then been shown to describe fluctuations in a great deal of problems, including the KPZ equation in various situations [39, 40, 41, 42]. One of these cases is of particular relevance to the present context, as it concerns the infinite TASEP with the step initial condition where all sites to the left of the origin are occupied, and all sites to the right are empty. Johansson [39] has shown by rigorous means that the time $t_{N}$ at which the $N$-th rightmost particle arrives at the origin scales as

[TABLE]

for $N$ large enough. An identical result has been obtained more recently [43] for the interface propagating on a quantum chain of free fermions prepared in a step initial condition. A similar scaling law has been found for the anchored Toom interface [44]. There, stationary height fluctuations obey a scaling law analogous to (2.25) or (2.27), involving the variable $\xi_{\rm GOE}$ distributed according to the Tracy-Widom $F_{1}$ law, characteristic of the Gaussian Orthogonal Ensemble (GOE) [45].

The ansatz (2.26) therefore amounts to stating that the lifetime fluctuations of a large droplet are distributed as if the full process (2.14) were described by Stage I. During Stage I, the active region of the system indeed behaves exactly as an infinite TASEP launched from a step initial condition. In this form, the above ansatz is made plausible by the observation, already made in section 2.3, that the dynamics in Stage II is in several respects the continuation of Stage I. In particular, the same parabolic interface (2.20) holds throughout both stages.

The above ansatz is strongly corroborated by the results of numerical simulations on the TASEP. Figure 6 shows plots of $4N-\langle t^{\star}\rangle$ and $({\mathop{\rm var}\nolimits}\,t^{\star})^{1/2}$ against $N^{1/3}$ for TASEP sizes $2N$ ranging from 10 to 800. Both datasets are rather well represented by linear fits of the form $aN^{1/3}+b$, but fits of the form $AN^{1/3}+BN^{-1/3}$ are almost perfect. The latter fits, shown as full lines, yield the estimates

[TABLE]

for the mean and variance of $\chi$. These numbers are in perfect agreement with the predictions of the ansatz (2.26), i.e. [46],

[TABLE]

Figure 7 shows a plot of the probability distribution $f(\chi)$ of the variable $\chi$, as measured on the TASEP with $N=100$ and $N=200$. Data are gathered over $10^{7}$ samples for each system size. The distribution predicted by the ansatz (2.26) is shown for comparison. It reads $2^{-4/3}f_{2}(2^{-4/3}\chi)$, where $f_{2}$ is the density of the Tracy-Widom $F_{2}$ law. A table of numerical values of $f_{2}$ is available at the web site [47]. The very smooth convergence of the data to the proposed limiting form provides a strong validation of the ansatz (2.26).

3 Finite-temperature zero-field dynamics

We now turn to the finite-temperature situation. In this section we consider the dynamics in the absence of a magnetic field, defined by the flipping rate (1.3), while the finite-temperature dynamics in a field, defined by the rate (1.7), is studied in section 4.

Finite-temperature dynamics creates and annihilates bulk excitations. Therefore, at variance with the zero-temperature situation considered in section 2, it is impossible to derive a closed effective dynamics for single interfaces. We are thus led to study the full spin model by means of numerical simulations. All data will be presented as a function of physical time $\alpha t$, the time scale $\alpha$ being small enough that all rates are smaller than unity, and can be used as flipping probabilities in a random sequential updating scheme.

3.1 Overall picture

We shall focus our attention on the fate of a large square droplet of linear size $N$ consisting initially only of down spins, with $N$ multiple of 3, embedded in a square sea of size

[TABLE]

of up spins, with periodic boundary conditions, as depicted in figure 8. The initial mean magnetization is therefore $M(0)=1-2(N/N_{\rm sea})^{2}=7/25=0.28$. This value is much smaller than the spontaneous magnetization $M_{0}$ (see (3.3)) for all temperatures considered hereafter in numerical simulations.

Figure 9 shows the average shape of a droplet of size $N=90$ embedded in a sea of size $N_{\rm sea}=150$, for $T/T_{c}=0.5$, at times $\alpha t$ multiples of 20. The plotted shapes are obtained by averaging over 1000 thermal histories, in order to practically eliminate any effect of thermal noise. The dynamics of a large droplet concomitantly exhibits two distinct features:

• •

It is advected ballistically in the SW direction, i.e., along the direction of the bias. Let $v_{\rm cm}$ be the velocity of this center-of-mass motion.

• •

It shrinks and disappears after a lifetime $t^{\star}$, whose mean value also grows ballistically, i.e., linearly with $N$. We denote by $v$ the associated asymptotic shrinking velocity, defined according to (2.9). In order to compare $v$ with its exactly known zero-temperature limit $v_{0}$ (see (2.11)), we rewrite (2.9) as

[TABLE]

In order to picture later stages of the dynamics, it is convenient to advect back the droplet in the NE direction, according to some estimated approximate value of $v_{\rm cm}$, so that the center of mass of the droplet remains essentially immobile during its whole lifetime. Figure 10 show such plots of a droplet of size $N=90$. The left panel corresponds to $T/T_{c}=0.5$, with $v_{\rm cm}\approx 0.78\,v_{0}$ and times $\alpha t$ multiples of 60. The right panel corresponds to $T/T_{c}=0.7$, with $v_{\rm cm}\approx 0.935\,v_{0}$ and times $\alpha t$ multiples of 200.

As temperature increases from zero to the critical point, the center-of-mass velocity $v_{\rm cm}$ increases steadily, while the contours of the average droplet shape become progressively rounder, and the mean droplet lifetime $\langle t^{\star}\rangle$ increases very fast.

3.2 Asymptotic velocity

The goal of this section is to investigate the temperature dependence of the asymptotic shrinking velocity $v$. In the analysis of zero-temperature dynamics (see section 2), a key role was played by the droplet lifetime $t^{\star}$. At zero temperature, this lifetime is defined unambiguously as the stopping time where the TASEP reaches the absorbing configuration where all the particles are to the right of all the holes, i.e., where the spin system reaches one of its ferromagnetic ground states. In the present finite-temperature situation, it will again be very convenient to measure a stopping time in numerical simulations. The latter time must however be defined with care, because the system keeps fluctuating. We have chosen to define the stopping time $t^{\star}$ of a given run as the first time where the mean magnetization per spin $M(t)$ exceeds (in absolute value) the spontaneous Onsager magnetization $M_{0}$ of the infinite square lattice, given by [48]

[TABLE]

The stopping time $t^{\star}$ thus defined is a fluctuating quantity. The linear growth of its mean value with $N$, according to (3.2), allows us to measure the velocity $v$.

Consider first the situation where $T/T_{c}=0.5$. Figure 11 shows a plot of the ratio $\alpha\langle t^{\star}\rangle/(2N)$ against $N^{-2/3}$, for $N$ ranging from 12 to 120, i.e., $N_{\rm sea}$ ranging from 20 to 200. The choice of the abscissa axis is motivated by the zero-temperature result (2.25), implying that the first relative correction to the ballistic law is proportional to $N^{-2/3}$. The finite-temperature data are observed to exhibit the very same finite-size correction. They are indeed well fitted by a quadratic function of $N^{-2/3}$, shown as a full line, whose intercept yields $v_{0}/v\approx 2.529$ (see (3.2)), i.e., $v/v_{0}\approx 0.395$.

Repeating the above analysis for other values of temperature in the ordered phase, we have obtained estimates for the ratio $v/v_{0}$ which are plotted against $T/T_{c}$ in figure 12. The above procedure yields accurate values of $v$ up to $T/T_{c}=0.75$. At low temperature, $v$ is observed to depart from its zero-temperature value $v_{0}$ as

[TABLE]

with $C\approx 4$. A similar low-temperature correction in $\rme^{-2/T}$ is shared by the free energy of an interface parallel to the axes,

[TABLE]

another result due to Onsager [49], whereas for bulk thermodynamical quantities the leading low-temperature correction is usually in $\rme^{-4/T}$ (see (4.3)). In the vicinity of the critical temperature, the fast fall-off of the velocity is well described by the power law (3.7), shown as a full line. The critical exponent ${\Delta_{v}}\approx 2.5$ will be determined more accurately in section 3.3 by means of a finite-size scaling analysis.

3.3 Critical region

The critical region deserves a more careful analysis. It is clear from figure 12 that the velocity $v$ vanishes as $T_{c}$ is approached. We expect that this slowing down is described by a power law of the relative distance to the critical point,

[TABLE]

and we denote by ${\Delta_{v}}$ the corresponding critical exponent:

[TABLE]

At the critical temperature ($T=T_{c}$), the velocity $v$ vanishes, and so the mean stopping time $\langle t^{\star}\rangle$ grows faster than ballistically with the droplet size $N$. It is expected to diverge as

[TABLE]

where $z$ is the dynamical critical exponent of the model. Numerical values of $z$ reported in [50] for the anisotropic kinetic Ising model exhibiting generic coexistence studied in [6] suggest that the spatial anisotropy of the dynamics is irrelevant, i.e., that $z$ coincides with the dynamical critical exponent of the two-dimensional Ising model with conventional Glauber dynamics, for which the more accurate estimate $z\approx 2.167$ is available [51].

The power laws (3.7) and (3.8) can be embedded into the following finite-size scaling law:

[TABLE]

The critical behavior (3.8) of the stopping time corresponds to $F(0)$ being some finite number, while the critical slowing down (3.7) of the velocity corresponds to the power-law growth $F(x)\sim x^{\Delta_{v}}$ for $x\gg 1$.

In figure 13 we have plotted $y=N^{z}/(\alpha\langle t^{\star}\rangle)$ against $x={\varepsilon}N^{(z-1)/{\Delta_{v}}}$ for several $N$ and ${\varepsilon}$, in order to draw the finite-size scaling function $y=F(x)$. Taking the value $z\approx 2.167$ for granted, a convincing data collapse is observed for ${\Delta_{v}}=2.5$, which has been chosen to draw the figure. This value of ${\Delta_{v}}$ is corroborated by the good agreement between the data and the three-parameter fit of the form $F(x)=a+b(x+c)^{2.5}$ (full line). The latter fit is however more of a guide to the eye than an independent determination of the exponent ${\Delta_{v}}$. We have checked that the above analysis is robust in several respects. First, relaxing the hypothesis that $z\approx 2.167$, we consistently obtain a best data collapse for $z\approx 2.2$. Moreover, allowing a non-linear dependence of the scaling variable $x$ on the distance to the critical point, i.e., changing ${\varepsilon}$ either to ${\varepsilon}(1+A{\varepsilon})$ or to ${\varepsilon}/(1+B{\varepsilon})$, can possibly accommodate values of ${\Delta_{v}}$ ranging from 2.3 to 2.7, but hardly beyond these bounds. We can therefore confidently give the estimate

[TABLE]

We shall come back to this value of the critical exponent at the end of section 4.3.

4 Finite-temperature dynamics in a field

We now turn to the generic situation where the temperature $T$ and the magnetic field $h$ are non-zero, so that the flipping rate is given by (1.7), with parameters $\gamma$ and $\kappa$. Our main goal is to obtain an accurate determination of the phase-boundary fields $\pm{h_{\rm b}}(T)$ delimiting the region of bistability (see figure 1), especially in the low-temperature regime and in the vicinity of the critical point.

4.1 Overall picture

We again focus our attention on the fate of a large square droplet of down spins embedded in a sea of up spins (see figure 8). The following picture is expected to hold:

• •

For small enough values of the magnetic field ($|h|<{h_{\rm b}}(T)$), the situation is qualitatively similar to what happens for $h=0$ (see section 3). A large droplet disappears in a ballistic time. The corresponding shrinking velocity $v$, defined according to (3.2), now depends both on $T$ and $h$. It vanishes as $h\to-{h_{\rm b}}(T)$, where the up phase loses its stability, and therefore its ability to heal defects spontaneously. The dynamics around $h=-{h_{\rm b}}(T)$ is accordingly expected to be very slow.

• •

For large negative values of the magnetic field ($h<-{h_{\rm b}}$), the up phase of the sea has lost its stability, and so the system ends up in a homogeneous down phase.

• •

For large positive values of the magnetic field ($h>{h_{\rm b}}$), the down phase inside the droplet has lost its stability, and so the droplet may disappear faster than ballistically. In any case the system ends up in a homogeneous up phase.

4.2 Boundary field

The goal of this section is to study the temperature dependence of the phase-boundary field ${h_{\rm b}}(T)$ delimiting the region of bistability. Here again, it will be very convenient to measure a stopping time $t^{\star}$ by means of numerical simulations. At variance with the situation in the absence of a field, the magnetizations of the two nonequilibrium steady states corresponding to the up and down phases are not known exactly. This difficulty can be overcome by monitoring simultaneously two samples of size $N_{\rm sea}$ with the same thermal noise, albeit with two different initial conditions:

• (1)

a droplet of down spins of size $N$ in a sea of up spins, as depicted in figure 8;

• (2)

the uniform configuration where all spins are up.

We define the stopping time $t^{\star}$ of a given run as the first time where the magnetization per spin $M^{(1)}(t)$ of the system with droplet initial condition exceeds in absolute value the magnetization $M^{(2)}(t)$ of the initially homogeneous system. The comparison of the absolute values of the two magnetizations has the advantage of keeping the operational power of the chosen definition of $t^{\star}$ for magnetic fields (slightly) below $-{h_{\rm b}}$, where the up phase has turned from stable to long-lived metastable. We also introduce the splitting probability $p^{\star}$ [52] of the up phase, defined as the probability for the system with droplet initial condition to end in the up phase, i.e., in practice, to have $M^{(1)}(t^{\star})>0$.

Consider first the situation where $T/T_{c}=0.5$. Figure 14 shows the mean stopping time $\alpha\langle t^{\star}\rangle$ (left) and of the splitting probability $p^{\star}$ (right), as functions of the magnetic field $h$, for several droplet sizes $N$. As $N$ increases, the stopping time exhibits peaks which get higher and narrower, whereas the splitting probability crosses over more and more steeply from 0 to 1. These observations are fully consistent with the overall picture of section 4.1. For a given droplet size $N$, the value $h^{(N)}$ of the magnetic field where the stopping time is maximal is virtually identical with the middle point where $p^{\star}=1/2$. The values $h^{(N)}$ thus defined shift slowly to the right as $N$ increases, toward a limit which is nothing but $-{h_{\rm b}}(T)$. This yields ${h_{\rm b}}\approx 0.35$ for $T/T_{c}=0.5$.

In order to look at the problem from a different perspective, we have measured the dependence of the velocity $v$ on the magnetic field $h$. To do so, data for $\alpha\langle t^{\star}\rangle$ have been extrapolated, along the lines of section 3.2, by means of quadratic fits in the variable $N^{-2/3}$. Figure 15 shows plots of the ratio $2N/(\alpha\langle t^{\star}\rangle)$ against $h$ for several $N$, as well as the extrapolated values yielding a plot of the ratio $v/v_{0}$ against $h$. The dashed straight line suggests that $v/v_{0}$ vanishes linearly with the distance $h+{h_{\rm b}}$ to the boundary of the region of bistability. We thus consistently recover ${h_{\rm b}}\approx 0.35$ for $T/T_{c}=0.5$.

Finally, the splitting probability $p^{\star}$ is observed to obey a finite-size scaling law of the form

[TABLE]

near the phase-boundary field $-{h_{\rm b}}$, where the crossover exponent 2/3 is again suggested by the relative size of the fluctuating term in (2.25), stemming from the KPZ theory. The data of figure 14 for $p^{\star}$ are replotted against $x=N^{2/3}(h+{h_{\rm b}})$ in figure 16, where the value ${h_{\rm b}}=0.35$ has been chosen in order to provide the best data collapse. This value of the phase-boundary field agrees with the estimates drawn from the observations made above on figures 14 and 15. This finite-size scaling method is however more accurate, as the estimated uncertainty on ${h_{\rm b}}$ is smaller than 0.005. The sigmoidal fit shown as a full line is meant as a guide to the eye.

The finite-size scaling ansatz (4.1) explains the slow convergence of the $h^{(N)}$ observed in figure 14, as it predicts $h^{(N)}+{h_{\rm b}}\sim N^{-2/3}$. The corresponding maxima of the stopping time are observed to obey an effective power law of the form $t^{\star}_{\rm max}\sim N^{1.2}$, with an apparent dynamical exponent much smaller than the exponent $z\approx 2.167$ entering the zero-field critical law (3.8). In any case, the dependence of $t^{\star}_{\rm max}$ on the system size $N$ is expected to be more complex than a single power law. A good deal of works [53, 54, 55, 56, 57], elaborating on the well-known droplet nucleation theory [58, 59], have indeed revealed that several time scales enter the dynamics of metastable phases in finite systems.

The finite-size scaling analysis of the splitting probability, yielding the most accurate determination of the phase-boundary field ${h_{\rm b}}(T)$, has been repeated for other values of temperature in the ordered phase. This procedure yields accurate values of ${h_{\rm b}}$ up to $T/T_{c}=0.75$. Figure 17 shows a plot of the estimates thus obtained against $T/T_{c}$, providing a quantitative version of the phase diagram shown schematically in figure 1. At low temperature, the phase-boundary field is observed to behave linearly as

[TABLE]

with $B\approx 1.63$, with the next correction being presumably exponentially small. The rationale behind this formula is that the two low-temperature expansion variables

[TABLE]

are proportional to each other along the phase boundary line, as $b\approx qa$ translates to (4.2) with $B=(\ln q)/2$. The fast fall-off of the phase-boundary field as $T_{c}$ is approached is described by the power law (4.4).

4.3 Critical region

The critical region again deserves a more careful analysis. As observed long ago [2, 6], the phase-boundary field ${h_{\rm b}}$ vanishes as a power law as $T_{c}$ is approached, as

[TABLE]

The critical exponent ${\Delta_{h}}$ is larger than its counterpart ${\Delta_{v}}$ by a sizeable amount. This is evidenced in figure 18, showing logarithmic plots of the zero-field velocity $v$ and the phase-boundary field ${h_{\rm b}}$ (data taken from figures 12 and 17) against ${\varepsilon}$. The data are observed to obey two different power laws. The slopes ${\Delta_{v}}=2.5$ and ${\Delta_{h}}=3.2$ of the straight lines are the outcomes (3.10) and (4.10) of more accurate finite-size scaling analyses. Least-square fits to the three to five leftmost data points however essentially yield the same exponents, especially so for ${\Delta_{h}}$.

In order to measure the exponent ${\Delta_{h}}$ with good accuracy, we again make use of finite-size scaling theory. An appropriate observable is provided by the slope at $h=0$ of the data plotted in figure 15, i.e.,

[TABLE]

This quantity becomes the reduced mobility

[TABLE]

in the thermodynamic limit. It has several advantages from a numerical viewpoint: it is very sensitive to a small magnetic field, and can be measured accurately for droplet sizes up to $N=120$. The nearly linear behavior of the data shown in figure 15 suggests that $\mu$ scales as the ratio $v/{h_{\rm b}}$, and therefore that it diverges as

[TABLE]

as the critical is approached after taking the limit of a very large system. The mobility $\mu$ therefore gives direct access to the exponent difference ${\Delta_{h}}-{\Delta_{v}}$.

The power law (4.7) extends to the finite-size scaling law

[TABLE]

where $x={\varepsilon}\,N^{(z-1)/{\Delta_{v}}}$ is as in (3.9), and the scaling function falls off as $G(x)\sim x^{-({\Delta_{h}}-{\Delta_{v}})}$ as $x\gg 1$. This finite-size scaling law is demonstrated by figure 19, where $y=N^{(z-1)({\Delta_{h}}-{\Delta_{v}})/{\Delta_{v}}}/\mu$ is plotted against $x_{1}=(1-A{\varepsilon})x={\varepsilon}(1-A{\varepsilon})N^{(z-1)/{\Delta_{v}}}$. In the definition of $x_{1}$, ${\varepsilon}$ has been changed to ${\varepsilon}(1-A{\varepsilon})$ in order to empirically improve convergence. Taking $z\approx 2.167$ and ${\Delta_{v}}\approx 2.5$ for granted, the best data collapse is observed for $A\approx 0.84$ and ${\Delta_{h}}-{\Delta_{v}}\approx 0.7$. The latter value is consistent with the rather good agreement between the data and the three-parameter fit of the form $y=1/G(x)=a+b(x+c)^{0.7}$ (full line). Here again, the latter fit is more of a guide to the eye than an independent determination of the exponent difference ${\Delta_{h}}-{\Delta_{v}}$.

We can state the outcome of the above analysis as

[TABLE]

Combining this with ${\Delta_{v}}\approx 2.5\pm 0.2$ (see (3.10)), we obtain finally

[TABLE]

The above predictions can be compared with various results in the literature. On the numerical side, we know of no previous measurement of ${\Delta_{v}}$, and of only two estimates of ${\Delta_{h}}$, namely ${\Delta_{h}}\approx 3.0\pm 0.4$ [2], and ${\Delta_{h}}\approx 2.7\pm 0.2$ [6]. Our estimate (4.10) agrees well with the first of those results, albeit only marginally with the second. It should be noticed that the error bar given in (4.10) is not reduced with respect to those given in earlier works, in spite of the good apparent quality of the finite-size scaling plots of figures 13 and 19.

It is also worth putting our results for the critical exponents ${\Delta_{h}}$ and ${\Delta_{v}}$ in perspective with a heuristic renormalization-group approach put forward in [6]. This line of thought yields

[TABLE]

where $y_{h}=2\nu-\beta=\beta+\gamma=\beta\delta=15/8$ and $y_{t}=1/\nu=1$ are the usual scaling dimensions of the two-dimensional Ising model, whereas $y_{i}<0$ is the scaling dimension of some irrelevant operator left undetermined in [6]. The same line of thought yields

[TABLE]

taking again $z\approx 2.167$ for granted, and so the exponent difference is predicted to be

[TABLE]

irrespective of the scaling dimension $y_{i}$. Our numerical result (4.9) is in perfect agreement with this prediction. Moreover, (3.10) yields the estimate

[TABLE]

for the irrelevant scaling dimension entering the heuristic approach of [6].

4.4 Magnetization curves

To close our numerical investigation, we show in figure 20 magnetization curves of the model for three temperatures in the ordered phase ($T/T_{c}=0.4$, 0.5 and 0.6). For each temperature, the mean magnetization $M$ of the nonequilibrium steady state of the up phase is plotted against the magnetic field $h$. Data are taken for samples of size $N_{\rm sea}=200$, with uniform initial condition (2), so that size effects are entirely negligible. The spontaneous magnetization $M_{0}(T)$ (see (3.3)) is recovered at $h=0$ (open symbols), as should be. Full symbols show the values of the phase-boundary field $-{h_{\rm b}}(T)$ and of the corresponding magnetization ${M_{\rm b}}(T)$. From a merely numerical perspective, the magnetization curves can easily be continued slightly below $-{h_{\rm b}}(T)$, where the up phase has turned from stable to metastable but very long-lived.

The magnetization curves are observed to be very regular near the phase-boundary field. This corroborates the commonly accepted expectation that nothing dramatic occurs to the pure phases near the limits of their region of stability.

It is pleasure to thank Paul Krapivsky, David Mukamel, Tomohiro Sasamoto and Gunter Schütz for interesting discussions.