Discontinuous non-equilibrium phase transition in a threshold Schloegl model for autocatalysis: Generic two-phase coexistence and metastability

Discontinuous phase transitions, and associated phase coexistence, have been studied extensively in equilibrium thermodynamic systems of infinite size. Here, the Gibbs phase rule implies that for a system with n-control parameters, phase coexistence can only occur in an (n − 1)-dimensional region of parameter space where chemical potentials of the two phases are equal. Nonetheless, challenges remain regarding a complete characterization of metastability or hysteresis and associated nucleation phenomena.^1,2 For non-equilibrium systems where the rates for transitions between states do not satisfy detailed-balance, there is even less understanding regarding fundamental behavior given the lack of a thermodynamic free energy framework.^3,4 Characterizing behavior in these systems has been presented as a “scientific grand challenge.”⁵ In particular, the lack of applicability of the Gibbs phase rule opens the possibility for phase coexistence in a finite region of parameter space with non-zero “volume” (i.e., for finite range of a single control parameter or in a region with non-zero area rather than just along a curve in a two-dimensional parameter space). This phenomenon is described as “generic two-phase coexistence” (generic 2PC) or sometimes as generic non-ergodicity.^6–12 

A general strategy to assess generic 2PC is to explore the evolution of a large droplet of one phase embedded in the other phase. Generic 2PC will be taken to imply that in a finite region of parameter space, the radius of the droplet shrinks at roughly constant rate and thus, ultimately disappears for both choices of embedded phase.^7,8,10 It is straightforward to relate this feature of generic 2PC to an orientation-dependence in the propagation of planar interfaces between the phases,^7,8,10 as is discussed further below.

Stochastic systems where such behavior might be revealed are those exhibiting bistability in a mean-field (MF) description, but where there is an underlying discrete spatial lattice structure to break rotational isotropy, as well as strong fluctuations. Such systems might include catalytic CO-oxidation on crystalline surfaces with limited surface mobility at high pressure^13–15 or spatially discrete population dynamics or epidemic models.^16,17

Perhaps, the most simple and natural models in which to explore non-equilibrium phase transitions are threshold versions of single-component Schloegl models for autocatalysis¹⁸ or the equivalent threshold contact processes.^17,19–24 In these Schloegl models, particles, X, reside on the sites of an infinite lattice. Particles spontaneously annihilate, X → ∅, at rate, p. Particles are created autocatalytically at empty sites if N or more particles reside on adjacent nearest-neighbor (NN) sites, ∅ + N X → (N + 1)X. For spatially homogeneous states, we will let C denote the concentration of particles, i.e., the probability that a site is occupied.

The case N = 1 corresponds to the standard contact processes, which is known to undergo a continuous transition from a populated phase or state with 0 < C ≤ 1 to a vacuum phase or state with C = 0 above a critical value of p.³ Based on a mean-field analysis,¹⁸ choices with N ≥ 2 can instead potentially exhibit a discontinuous transition between populated and vacuum states where C jumps discontinuously to zero upon increasing p. We will just consider the case N = 2 (a quadratic contact process) for a square lattice. There are various possible prescriptions of the autocatalytic configurations and associated creation rates. The conventional threshold model for N = 2 considered in this study, and previously by other groups, assigns a single creation rate of unity for empty sites with any configuration of N ≥ 2 particles on NN sites,^20,22–24 see Fig. 1. The Durrett threshold model for N = 2 assigns a creation rate equal to k/4, where k is the number of diagonally neighboring pairs of particles adjacent to the empty site.^9,17,25 There can be k = 4, 2, 1, or 0 such pairs as apparent from the 1st, 2nd, 3rd, or 4th-6th configurations for particle creation in Fig. 1, respectively. Thus, the Durrett creation rate is normalized so that its maximum value is unity for an empty site with all four NN sites occupied.

It is appropriate to mention that both conventional and Durrett N = 2 threshold models on a square lattice incorporate an “anomaly:” the populated state cannot displace the vacuum state separated from it by a vertical or horizontal interface, and a droplet of the populated state embedded in the vacuum state can never expand beyond a rectangular region circumscribing the initial droplet.^{10–12,17,23} For this reason, it will be instructive to also consider more general models “perturbed” from the above basic models to allow spontaneous creation of particles at empty sites with a (small) rate ε > 0. This perturbation removes the above-mentioned anomaly.

Regarding the possible occurrence of generic 2PC in stochastic models, this phenomenon was rigorously demonstrated for Toom’s model of voter dynamics.^6–8 Generic 2PC was also observed in simulations of the Durrett threshold model for N = 2^10–12 and in perturbations thereof,^26,27 as well as in models involving interface pinning.⁹ However, there is a perception that generic 2PC in Toom or Durrett type models derives from the unusual choices of rates in those models.^23,28 This perception is buoyed by the feature that a discontinuous transition at a single point apparently occurs in the classic Ziff-Gulari-Barshad (ZGB) reaction model,²⁹ in variants of ZGB-type models with large diffusion^30,31 or reaction³² lengths, and also in contact processes with large diffusion²⁷ or interaction^33,34 lengths. Furthermore, there appears to be a general expectation that for a conventional threshold Schloegl model or contact process, generic 2PC should not occur, so that the Gibbs phase rule applies.^23,35,36 In fact, for the conventional N = 2 threshold model on a square lattice shown in Fig. 1, a previous simulation study by da Silva and de Oliveira²³ suggested that in the limit of large system size, a discontinuous transition occurs from the populated to the vacuum state and that phase coexistence occurs at a unique value of p = p_e ≈ 0.35 in the spirit of the Gibbs rule. In addition, it was suggested in Ref. 23 that critical behavior with divergent fluctuations³⁷ occurs at this discontinuous transition. We will dispute both claims.

Next, we provide further discussion of: (i) droplet dynamics as applied to characterize phase stability and generic 2PC; (ii) consequences of and elimination of the above-mentioned anomaly in the N = 2 threshold models; (iii) the relationship between orientation-dependence of planar interface propagation and generic 2PC; and (iv) issues pertaining to metastability and criticality at a discontinuous phase transition.

(i) Dynamics of embedded droplets. Applying the above-mentioned criterion to our threshold Schloegl models, generic 2PC implies that over a finite range of the particle annihilation rate p, the radius of an arbitrarily large droplet both of the populated state embedded in the vacuum state, and of the vacuum state embedded in the populated state, must shrink at roughly constant rate. One might anticipate that outside the region of generic 2PC, an arbitrarily large droplet of the less (more) stable phase would generally shrink (grow). However, the anomaly in our N = 2 threshold models on a square lattice means that the populated droplet can never grow and thus, must ultimately disappear if for no other reason as a result of fluctuations. Indeed, such slow fluctuation-mediated decay of embedded populated droplets must occur for low p. For this reason, we enforce the stronger requirement of shrinking at roughly constant rate rather than a weaker criterion of just shrinking with certainty (even though the former does not apply for curvature-driven shrinkage in thermodynamic systems at phase coexistence). With the weaker criterion, the vacuum state would be regarded as stable for low p, where the populated state is expected to also be stable, i.e., one would have “artificial” generic 2PC in this low p regime.

(ii) Elimination of the anomaly in N = 2 threshold models. A simple strategy to eliminate the feature that the populated droplets cannot expand into the vacuum state is to “perturb” the basic model to incorporate autocatalytic creation of particles at (small) rate ε > 0. (An alternative would be to introduce particle diffusion.) Such perturbations also allow the populated state to displace the vacuum state separated from it from by a horizontal or vertical interface. The perturbed models have a two-dimensional parameters space, and now generic 2PC will exist in a region of finite area. Also outside this region, a sufficiently large droplet of the stable state will now always grow.

(iii) Orientation-dependent interface propagation. A signature feature which can control droplet dynamics and lead to generic 2PC is an orientation-dependence of the propagation and stationarity of planar interfaces separating different states. This feature was already recognized for Toom’s model.^6–8 For threshold models, it is convenient to define an equistability point, p = p_eq, corresponding to a stationary planar interface. Then, the populated state displaces the vacuum state for p < p_eq, and the opposite applies for p > p_eq. Suppose that p_eq increases from p⁻ to p⁺ changing the interface orientation from horizontal (or vertical) to diagonal. Then, for p⁻ < p < p⁺, an embedded vacuum droplet will evolve to a diamond shape and then shrink (and an embedded populated droplet will evolve to a square shape and then shrink) at roughly constant rate (see later examples) demonstrating generic 2PC. Thus, rather than just conventional simulations of homogeneous states in threshold models,^23,31 analysis of interface and droplet evolution in heterogeneous systems is invaluable. (One caveat to the above description for the threshold models is that a vertical or horizontal interface is stationary for all p > p_eq since the populated state cannot displace the vacuum state, but the vacuum state still displaces the populated state for p > p_eq.)

(iv) Metastability versus criticality at discontinuous transitions. For any discontinuous transition, a basic issue is whether there is the usual non-zero effective line tension at the interface between coexisting states,^1,2 or perhaps criticality with divergent fluctuations.³³ A non-zero line tension ensures that fluctuations remain finite at the transition, as in thermodynamic phase transitions, and implies metastability. We propose that the conventional threshold model (as well as Durrett’s version) for N = 2 has a non-zero line tension rather than criticality. We do emphasize the difficulty in developing a rigorous treatment of metastability even in thermodynamic systems^38,39 let alone non-equilibrium systems.^11,27,40,41 For example, it is expected that metastable states cannot be defined as unique analytic extensions of stable states.

In this paper, we present both Kinetic Monte Carlo (KMC) simulations and an approximate analytic treatment for the conventional N = 2 threshold model on a square lattice. Previous simulations in Ref. 23 considered only homogeneous states, but here we exploit both simulations and analytic theory to characterize heterogeneous states (specifically planar interfaces and droplets), which proves invaluable for elucidating generic 2PC. Our previous studies considered exclusively the distinct Durrett version of the N = 2 threshold model but did explore heterogeneous states. It should be noted that simulations exist of other realizations of N = 2 Schloegl models⁴² with distinct behavior from that described here. Model behavior depends on the realization, so results presented below are specific to our realization. In Sec. II, we present the exact hierarchical master equations for the conventional threshold model for both homogeneous and heterogeneous states. We also describe mean-field type truncation approximations. In Sec. III, we will present KMC simulation results which precisely characterize both steady-state and kinetic behavior and also the behavior associated with interfaces and droplets. We present a precise analysis of equistability of planar interfaces, sometimes exploiting finite-size-scaling and refined constant-coverage simulation algorithms. Utilizing new algorithms, we also analyze critical droplet size in the context of nucleation theory. In Sec. IV, we present a heterogeneous mean-field analysis of the above issues (including the first such analysis of critical droplet size), predictions of which are compared with KMC results. Section V provides our conclusions.

First, we present the exact form of the hierarchical master equations for spatially homogeneous states in the conventional N = 2 threshold model on an infinite square lattice, allowing for perturbation to include spontaneous particle creation. Let x (o) indicate a filled (empty) sites and P’s denote probabilities of various configurations, so P[x] = C, P[o] = 1 − C for single-site configurations, P[x x] + P[x o] = P[x] = C for pairs of sites, etc. Accounting for spontaneous annihilation of particles at rate p, possible spontaneous creation at empty sites at (small) rate ε, and autocatalytic creation at empty sites at rate 1 for site with two or more neighboring particles, one has that

This equation is not closed, but one can develop equations for probabilities of pairs and larger configurations of sites to form an infinity hierarchy. To truncate the hierarchy, the mean-field site approximation neglects all spatial correlations setting $P x x o x x =P [ x ] 4 P [ o ]$⁠, etc. The refined Kirkwood-type pair approximation retains pair probabilities setting $P x x o x x =P [ x o ] 4 /P [ o ] 3$⁠, where the denominator compensates for over-counting of the central empty site, etc.^3,12,25,43

In the site approximation, the time-evolution of the concentration is given by

Our main interest is in steady-state behavior for ε = 0 where in the site-approximation,²³ 

From (3), one determines the upper spinodal boundary for the populated state as p_s(site) = 0.815 423.²³ A corresponding pair approximation analysis is most naturally formulated in terms of the conditional concentration, K, for a site to be occupied given an adjacent empty site.^12,25 Then, steady-state K for ε = 0 satisfies (cf. Ref. 23)

leading to the pair approximation estimate of upper spinodal, p_s(pair) = (10 + 7^2/3)/54 = 0.528 153. We shall see in Sec. III that these estimates are significantly above the value obtained from KMC simulation. However, mean-field site and pair predictions for the steady-state C(p) versus p are actually very close to the KMC results for p up to about 0.25, see Sec. III for a detailed comparison. Interestingly, a more detailed analysis for the steady-state for ε = 0 shows that for the populated state

where it is only in the O(p⁶) term that a difference appears between site (−77p⁶) and pair (−113p⁶) results. The dominant small-p steady-state behavior, C(p) ≈ 1 − p, is clear as in this regime, the model corresponds to random particle annihilation and creation.

The above master equations and approximations naturally generalize to spatially heterogeneous states. Let x_i,j (o_i,j) indicate that site (i,j) is filled (empty), and let P’s denote probabilities for various configurations of sites. Thus, P[x_i,j] = C_i,j and P[o_i,j] = 1 − C_i,j denote the probability that site (i,j) is occupied and empty, respectively, P[x_i,j x_i,j+1] denotes the probability of a filled pair, etc. Accounting for creation and loss of particles at site (i,j), one has that

In (6), for the autocatalytic particle creation (gain) terms at empty site (i,j), we have explicitly listed contributions from the configuration with four particles NN to site (i,j), one out of four terms with three NN particles, but left implicit six additional terms with two NN particles. One can also develop equations for the evolution of probabilities for multisite configurations to generate an infinite hierarchy. The site and pair approximations used to truncate the hierarchy for homogeneous states naturally extend to this heterogeneous case.

The truncated heterogeneous hierarchical master equations reduce to lattice differential equations (LDEs)^44,45 which might also be regarded as discrete reaction-diffusion type equations.^{12,25,46–49} In the mean-field site approximation for ε = 0, from (6) and using R(C) from (2), one obtains

In the expression for the gain terms after the first equality in (7), we first explicitly list contribution for four particles NN to site (i,j), then for one of four terms with three NN particles, but leave implicit six additional terms with two NN particles. The expression after the second equality extracts separately the reaction kinetics term, R(C_i,j), to more explicitly show the nature of the spatial coupling (noting that all the coupling terms vanish for homogeneous states). Clearly, these LDEs do not have the simple and most common discrete Nagumo form, where the spatial coupling has the generic form DΔC_i,j involving a discrete Laplacian ΔC_i,j = C_i+1,j + C_i,j+1 + C_i−1,j + C_i,j−1 − 4C_i,j. The spatial coupling in (7) has more complex non-linear form in part reflecting a concentration-dependent effective diffusivity (cf. Refs. 12 and 25, and 46–49).

For analysis of droplet dynamics, one must use the general form of the above equations. However, some simplification is possible for the analysis of the evolution of planar interfaces. For example, for diagonal interfaces of slope S = 1, one has that C_i,j = C_i−j depends only on i − j and (7) reduces to

It is clear that for planar interfaces, the spatial coupling structure in the LDE depends upon the orientation S and thus so should interface propagation and stationarity.^12,25,43 We shall see that this feature recovers actual behavior for the conventional N = 2 threshold model. In the pair approximation, one has coupled LDEs for single-site and pair quantities. We do not present the complicated form of these equations, but the reader is referred to Refs. 12 and 25, and 49 for analogous examples for the Durrett N = 2 threshold model.

For the conventional threshold model with N = 2 andε = 0 on a square lattice, KMC simulations for specified constant-p can follow the time-evolution of the model to determine stable steady-states and also to assess metastability. Such an analysis of steady-states was performed previously in Ref. 23 to determine the concentration, C(p), of the populated steady-state versus p. Our own simulations starting from a state with C = 1 are entirely consistent with these previous results in the limit of large system size. Results for this steady-state concentration, C(p), are shown as the upper solid black curve and symbols in Fig. 2 for p < p_e ≈ 0.35. If one chooses p just above p_e, the system remains in a metastable populated state with high C-value before slowing evolving to the vacuum state. This metastable extension of the stable state is also indicated in Fig. 2. This constitutes classic metastability behavior expected for a discontinuous phase transition. We also note that simple analytic approximations discussed in Sec. II accurately recover the simulated C(p) versus p at least for p up to about 0.25. The success of these approximations is expected not just because low-p behavior corresponds to essentially random particle creation (annihilation) at rate 1 (p), but also given the observation in Sec. I that site and pair approximations Taylor expansions for C(p) versus p agree to O(p⁶).

A particularly effective approach to assess this discontinuous transition is the constant-concentration (CC) ensemble simulation approach of Ziff and Brosilow.⁵⁰ One identifies a target concentration, C_t, picks a site at random, and attempts to create (annihilate) a particle at that site if the current C < C_t (C > C_t). Creation only occurs if the site is empty with two or more populated NN and annihilation only occurs if the site is populated. The number of creation attempts divided by the number of annihilation attempts (whether successful or not) averaged over a long simulation is set equal to 1/p_t. Here, p = p_t is the p-value corresponding to C = C_t. CC-simulation should generally match constant-p simulation for large system size. For systems where C jumps discontinuously from C₋ to C₊ between low and high concentration states at a unique p = p_eq, choosing C_t anywhere in the middle of the range C₋ < C_t < C₊ will recover p_t = p_eq. Also, the simulated long-time state should correspond to a strip configuration of one state separated from the other by a planar interface. However, for the N = 2 threshold model, judiciously starting the simulation with a strip geometry reveals that p_t = p_eq depends on the orientation of the strip, as described below. (This feature is obscured starting with a random configuration.²³)

Selection of an initial strip geometry with prescribed orientation is readily achieved for a vertical strip (slope S = ∞) or horizontal strip (S = 0) or a diagonal strip (S = 1) using an L × L site square simulation cell with periodic boundary conditions (PBCs). Strips with other slopes, S, can be generated using a rectangular simulation cell with PBC with the appropriate aspect ratio, S. Note that while the strip edges fluctuate during the simulation, the overall orientation or slope is preserved by the periodic boundary conditions. Performing CC simulations for a range of initial strip orientations allows precise determination of the S-dependence of the equistability p = p_eq(S). We find that $p e = max S p eq ( S ) = p eq ( S = 1 ) =0.351 78 ( 2 )$ is slightly above p_eq(S = 2) = 0.3497(1) and further above $p f = min S p eq ( S ) = p eq ( S = 0 or ∞ ) =0.3161 ( 1 )$⁠, see Appendix  A for details. Stationary interfaces for the two limiting cases S = 1 and S = 0 or ∞ are illustrated in the insets of Fig. 2.

In Sec. III C, we discuss in detail how the orientation-dependence of interface propagation associated with the above S-dependence of p_eq ensures that arbitrarily large droplets of either phase embedded in the other shrink at roughly constant rate for p_f < p < p_e. Thus, our conclusion is that this conventional N = 2 threshold model with ε = 0 exhibits generic 2PC for p_f < p < p_e, where both the populated and vacuum states are regarded as stable. In addition, we note that for p > p_e, the vacuum state is the unique stable steady state. For p < p_f, the populated state is stable. The absorbing vacuum state can also persist indefinitely for p < p_f. However, we do not characterize it as stable in the usual sense as droplets of the populated state are only eliminated very slowly as a result of fluctuations. We have noted above how a discontinuous transition occurs from the populated state to the vacuum state upon increasing p above p_e. It is also tempting to associate the lower boundary, p_f, of the generic 2PC region with a discontinuous transition from the vacuum state to the populated state. However, the absorbing nature of the vacuum state precludes this simple interpretation, an anomaly which is resolved by consideration of the perturbed model below.

A broader perspective on generic 2PC, and further demonstration of its robustness, comes from consideration of the more general perturbed version of the conventional N = 2 threshold model with spontaneous particle creation at empty sites with (small) rate ε > 0. In this model, the discontinuous transition persists for 0 ≤ ε < ε_c, where ε = ε_c ≈ 0.006 corresponds to a non-equilibrium critical point, as does generic 2PC for a range p_f(ε) < p < p_e(ε). (The critical point is expected to be in the Ising universality class.²⁶) For 0 < ε < ε_c, the discontinuous transition now occurs between a highly populated and a relatively unpopulated (near-vacuum) states. Specifically, a discontinuous transition occurs from highly populated to the near-vacuum state upon increasing p above p_e(ε), and a discontinuous transition occurs from the near-vacuum to the highly populated state upon decreasing p below p_f(ε). For p < p_f(ε) [p > p_e(ε)], the highly populated [near-vacuum] state is the unique stable steady state, and for p_f(ε) < p < p_e(ε), both states are stable. Also, just as there is a metastable extension of the highly populated state extending above p_e(ε), analogous to that characterized above for ε = 0, there is also a metastable extension of the near-vacuum state extending below p_f(ε). This extension is indicated in the inset of Fig. 2.

Fig. 3 shows the variation of the equistability p_eq(S) with ε for 0 ≤ ε ≤ ε_c and for different interface orientations S = 1, 2, 3, and ∞ (or 0). For each ε in the range 0 < ε < ε_c (as for ε = 0), p_eq(S) decreases monotonically as S increases from 1 to ∞, where $p e ( ε ) = max S p eq ( S )$ and $p e ( ε ) = max S p eq ( S )$⁠. It follows that the droplets of the vacuum state embedded in the populated state will grow for p > p_e(ε) and shrink for p < p_e(ε); droplets of the populated state embedded in the vacuum state will grow for p < p_f(ε) and shrink for p > p_f(ε). Growth and shrinkage occur at non-zero velocity consistent with generic 2PC existence for p_f(ε) < p < p_e(ε), see Sec. III C. Significantly, one sees that p_f(ε) → p_f(0) = 0.3161 and p_f(ε) → p_f(0) = 0.351 78 (values listed above) smoothly as ε → 0, supporting our identification of generic 2PC for the conventional N = 2 threshold model for ε = 0.

For the conventional N = 2 threshold model for ε = 0, we now turn to the utilization of constant-p KMC simulations to assess the kinetics of “rapid” approach to the absorbing vacuum state for p well above the discontinuous transition. There is evidence of a somewhat ill-defined metastable populated state for some finite range of p just above the generic 2PC region. We denote this range by p_e < p < p_s(eff), where p_s(eff) is an effective spinodal point, see Fig. 2. A useful approach to assess p_s(eff) is to anticipate that the kinetics of evolution to the vacuum state well above the transition is primarily controlled by the distance δp_s = p − p_s(eff) > 0 above the effective spinodal. Specifically, we anticipate kinetics of the form,^11,27,30

with the function f(x) roughly independent of p. The feature that f(x) → 0 as x → ∞ ensures that C → 0, as t → ∞, and furthermore f(x) should vanishes exponentially as for mean-field chemical kinetics. An analysis of the kinetics for 0.378 < p < 0.388 for a 1024 × 1024 site system shown in Fig. 4 confirms this behavior and revealing the form of the universal function f(x) as well as yielding the estimate p_s(eff) = 0.369. We caution again that there is no precisely defined metastable state or spinodal, and correspondingly that the estimate of p_s(eff) does depend on somewhat the range of p-values selected. By analogy with behavior in kinetic versions of ferromagnetic Ising models, the regime δp_s > 0 corresponds to the so-called strong-field regime, where evolution is reminiscent of spinodal decomposition.⁵¹ 

Next, we assess the more complex kinetics for p_e < p < p_s(eff) involving “slow” evolution from the metastable populated state to the vacuum state mediated by nucleation and growth of supercritical droplets of the stable vacuum state inside the metastable state. Below, we let C_m denotes the concentration of the metastable state, J_nuc denotes the nucleation rate for supercritical droplets, and V_grow denotes the characteristic radial growth velocity of such nuclei. Then, consistent with the general Johnson-Mehl-Avrami-Kolmogorov (JMAK) theory for nucleation kinetics,^52–54 as in previous studies of nucleation-mediated kinetics for non-equilibrium systems,^11,27,40 we anticipate that

and τ_nuc is the characteristic time for nucleation-mediated evolution to the vacuum state. Simulation data for kinetics for 0.362 < p < 0.374 in a 4096 × 4096 site system shown in Fig. 5 reveal the above sigmoidal Arvami form for lower p. Fitting the slopes of the ln C versus t³ data allows extraction of τ_nuc which varies from 1.2 × 10⁴ time units for p = 0.362 to 6.8 × 10² time units for p = 0.374 (the latter not strictly reflecting nucleation). Results for τ_nuc are shown in Fig. 6.

For a more detailed analysis, by analogy with thermodynamic systems, we anticipate a nucleation rate with the specific form J_nuc ∼ exp[ − c_nuc/δp] , where δp = p − p_e measures the driving force for nucleation.^11,27,40 It is natural to write c_nuc ∝ (σ_eff)², where σ_eff denotes the effective line tension at interface between coexisting states.^27,41 If one also accounts for the feature that the asymptotic growth velocity for large droplets (with near-planar interfaces) should scale like V_grow ∝ δp,^11,27,30,55 then it follows that τ_nuc ≈ b (δp)^−2/3exp[c_nuc/(3δp)].¹¹ Note that various refinements of the prefactor are possible,²⁷ but the dominant dependence on δp is given by the exponential factor. The appropriateness of this form for τ_nuc is clearly demonstrated in Fig. 6. Assessment of the key parameter, c_nuc, follows precisely from the slope in this near-linear variation of ln(δp^2/3τ_nuc) versus 1/δp. The key point is that the slope is unambiguously positive corresponding to c_nuc ≈ 0.125, and as a result, the effective line tension σ_eff is strictly positive in this conventional N = 2 threshold model. Thus, we conclude that this system does exhibit conventional metastability and nucleation behavior associated with a non-zero effective line tension. Finally, we remark that a previous analogous analysis for the Durrett N = 2 threshold model revealed that c_nuc ≈ 0.024¹¹ corresponding to a finite effective line tension roughly half the magnitude of that in this conventional N = 2 threshold model.

Such metastability behavior associated with a non-zero effective line tension is distinct from that at a critical discontinuous phase transition proposed in Ref. 23. Critical behavior has been associated with vanishing line tension in voter-type models⁵⁶ and also with diverging correlation lengths.³⁷ In Appendix  B, we reassess the fluctuation analysis in Ref. 23 for this conventional N = 2 threshold model.

For the conventional N = 2 threshold model with ε = 0, we consider first the fate of a vacuum droplet with roughly circular shape embedded in the populated state. For p_f < p < p_e, the horizontal and vertical edges (S = 0 and ∞) will grow outwards and ultimately disappear since p > p_f. Thus, the droplet will convert to a diamond shape bordered by diagonal edges with S = 1 and will thus shrink at constant rate since p < p_e. Such vacuum droplet evolution is demonstrated in Fig. 7(a) for p = 0.340. For p < p_f, all edge orientations of the vacuum droplet will shrink, and for p > p_e, all edge orientations will grow at constant rate (not shown). Thus, the populated state is stable for all p < p_e. For a populated droplet with roughly circular shape embedded in the vacuum state for some p in the range p_f < p < p_e, edges with diagonal orientations will grow outwards since p < p_e and ultimately disappear. Thus, the droplet will convert to a square shape bordered by horizontal and vertical edges with S = 0 and ∞ and will thus shrink at constant rate since p > p_f, see Fig. 7(b). For p > p_e, all edge orientations of the populated droplet will shrink at constant rate. Thus, the vacuum state is stable for all p > p_f.

For p < p_f, all orientations of the populated droplet with roughly circular shape embedded in the vacuum state will initially grow (not shown). However, the anomaly in the N = 2 threshold model means that this droplet can never expand outside a circumscribing rectangular region and thus, will ultimately disappear due to fluctuation effects. The populated droplet will survive a long time, its lifetime scaling exponentially with droplet radius (rather than linearly with radius) for p > p_f. Thus, on average, the droplet radius decreases more slowly than at roughly constant rate, and we do not regard the vacuum state as stable for p < p_f. All the above results imply that generic 2PC occurs only for p_f < p < p_e.

Finally, we briefly describe an additional droplet analysis related to nucleation-mediated evolution to the vacuum state for p_e < p < p_s(eff) and associated metastability of the populated state. Specifically, we characterize critical vacuum droplets in this regime. Droplet size can be defined using a percolation-theoretic framework as the number, M, of connected NN empty sites and an effective radius, R, then determined by M = πR². The basic picture is that there exists a critical size, R_c, such that smaller vacuum droplets with radius R < R_c should shrink, while larger droplets with R > R_c should grow. In the spirit of a conventional nucleation theory, we expect for this non-equilibrium model that

where again σ_eff is an effective line tension.^11,27,41 (Note that droplet fate can be somewhat ambiguous given the fluctuating background in contrast to droplets embedded in an absorbing state.⁵⁷)

Constant-p simulation can be used to track the evolution of vacuum droplets with different initial sizes and thus, to assess R_c for a specific p. More extensive analysis also yields R_c versus δp. However, a potentially more efficient alternative approach applies CC simulations to directly stabilize critical vacuum droplets.^26,27,30,50 We find that complications arise in this analysis since the size of the vacuum droplet is not fixed in CC simulation. One can average droplet size and annihilation rate, p (determined from the relative numbers of particle creation and annihilation attempts), over a long simulation. However, this ignores large fluctuations in droplet size and shape, which are strongly correlated with variations in p. (Fluctuations in size are positively correlated with those in the concentration in the surrounding metastable state.) Thus, instead, we first performed simulations averaged over a moderate time-period to get a scatter plot of R_c versus p which is centered on behavior corresponding to R_c ∼ σ_eff/δp. As an alternative, we have sampled over a short time periods tracking R_c at the end of the period and the average p over that period. We then average p values for all time-periods with the same R_c. This produces the same basic behavior as the approach samples for moderate time-periods but significantly reduces statistical scatter. Our results will be presented in Fig. 8, which reveal that 1/R_c does vary roughly linearly in δp although there are some deviations from simple proportionality for small δp. The inset shows analogous behavior for the N = 2 Durrett threshold model. The same behavior is also seen in our MF-type analyses in Sec. IV B.

For the conventional N = 2 threshold model with ε = 0, our mean-field analysis in Sec. II for uniform states has already determined site and pair approximation predictions for bistable region with coexisting populated and vacuum states. Our initial goal here is to apply the site and pair approximation versions of the associated LDE (discrete reaction-diffusion equations) to analyze the propagation of planar interfaces of various slopes, S, separating populated and vacuum states in the bistable region. We thereby determine the orientation-dependence of equistability, p_eq = p_eq(S). This yields mean-field predictions for the extent of generic 2PC in N = 2 threshold models. For details of our numerical implementation and analysis, the reader is referred to Refs. 25 and 49.

To determine the upper boundary, p = p_e = p_eq(S = 1), of the generic 2PC regime, we analyze the propagation and stationarity for diagonal interfaces with slope S = 1. LDE analysis at either the site or pair approximation level reveals that the velocity, V, of the interface versus p varies smoothly through zero at p = p_eq(S = 1), just as in the KMC simulations. From this LDE analysis, we obtain p_eq(site) = 0.702 84 and p_eq(pair) = 0.442 53 in the site- and pair-approximations, respectively.

To determine the lower boundary, p = p_f, of the generic 2PC regime, we claim that it is most appropriate to make the identification $p f = lim S → ∞$ p_eq(S) in the LDE analysis. The reason for assessing this limiting behavior rather than just directly evaluating p_eq(S = ∞) is that the analysis of interface propagation via LDE exhibits some anomalies not seen in stochastic or continuum reaction-diffusion equations (RDE) models. Specifically, interface pinning or propagation failure can occur for certain orientations including the principal lattice directions. Such unusual behavior in the LDE treatment might be anticipated for a perfectly vertical interface, S = ∞, which has no kinks in contrast to a fluctuating interface in the stochastic model. Indeed, we find that $lim S → ∞ p eq ( S ) < p eq ( S = ∞ )$⁠, see Refs. 25, 44, and 45, and 49 for further discussion of these subtle issues. A comprehensive LDE analysis of the orientation dependence of the equistability, p_eq(S), for our conventional N = 2 threshold model with ε = 0 reveals clear convergence with increasing S, as indicated in Table I. We conclude from this analysis that p_f(site) = 0.6847 and p_f(pair) = 0.4187 in the site- and pair-approximations, respectively.

We now summarize our mean-field results for p_e, p_f, and for also p_s (obtained previously in Sec. III A) and compare these with KMC predictions. Table II reports these quantities for both the conventional and Durrett N = 2 threshold models with ε = 0. Rates for autocatalytic particle creation in the Durrett model are generally below those in the conventional model as noted in Sec. I, but the spontaneous annihilation rate remains the same. Consequently, the spinodal points and the regime of generic 2PC occur at lower p-values in the Durrett model. Thus, we also report in Table II suitably rescaled widths of the generic 2PC regime and of the metastable regime. We find a significant improvement in the pair over the site approximation in the prediction of p_s, p_e, and p_f. However, the predicted width for the generic 2PC region is lower than precise KMC values, and for the metastable region, the predicted width is higher (the latter being a typical feature of mean-field type approximations).

First, we present a LDE analysis for the conventional N = 2 threshold model with ε = 0 of droplet dynamics in the generic 2PC regime, p_f < p < p_e. We naturally compare results with those from KMC simulation discussed in Sec. III C. We have performed an analysis for both the site-approximation and pair-approximation, and the behavior is qualitatively similar. Thus, in Fig. 9, we just show pair LDE results for the shrinkage of (initially) large droplets of one state embedded in the other state for p = 0.43 in the generic 2PC regime, p_f(pair) = 0.4187 < p < p_e(pair) = 0.4425. This behavior, including the detailed droplet shape evolution, is entirely analogous to that shown in Fig. 7 as assessed by KMC simulation.

For p = 0.43, the vacuum droplet in Fig. 9(a) adopts a diamond shape, the diagonal edges of which shrink since p < p_e. It will also shrink for p < p_f. However, for p > p_e, a sufficiently large vacuum droplet will grow (see below). For p = 0.43, the populated droplet in Fig. 9(b) adopts a square shape, the horizontal and vertical edges of which shrink since p > p_eq(S = ∞) = 0.4213 > p_f = 0.4187. Although less obvious, the populated droplet will also shrink for p_f < p < p_eq(S = ∞) given that the edges are somewhat curved and thus, only nearly horizontal or vertical. It will also clearly shrink for p > p_f. We also note that a sufficiently large populated droplet will survive as a stationary solution of the LDE for p < p_f (rather than undergo slow fluctuation-mediated shrinkage as in KMC simulations).

Finally, we present a LDE analysis for the conventional N = 2 threshold model with ε = 0 of critical vacuum droplets embedded in the populated state for p_e < p < p_s. Specifically, we assess the variation of the critical droplet size, R_c, versus δp = p − p_e. Fig. 10 shows an example obtained from the pair LDE for the growth (shrinkage) of a vacuum droplet which is above (below) the critical size. These results immediately determine upper and lower bounds on the critical size for the selected p. To obtain a precise estimate of the critical size, R_c, we systematically adjust the initial droplet size until we achieve a stationary droplet configuration for long times. Fig. 8 summarizes results of such extensive analysis for R_c for both the site- and pair-approximations. Comparing with precise KMC CC-simulation values, we find a significant improvement in the pair prediction for R_c over the site approximation. Also for comparison, we show in the inset of Fig. 8 corresponding and entirely analogous behavior for the Durrett model. In summary, both KMC analysis and site and pair LDE analyses provide evidence for a finite critical size of vacuum droplets. This observation is compatible with our earlier claim of metastability, associated with a non-zero effective line tension σ_eff > 0, in the regime p_e < p < p_s(eff).

Our analysis of the conventional N = 2 threshold model with ε = 0 provides clear evidence for a discontinuous phase transition between populated and vacuum states together with generic 2PC over the finite range of particle annihilation rate, 0.3161 = p_f < p < p_e = 0.3518. Within this range, a droplet (of any size) of the vacuum state embedded in the active state shrinks at roughly constant rate, as does a droplet of the active state embedded in the vacuum state. This behavior follows immediately from the demonstrated orientation-dependence of equistability for planar interfaces separating for coexisting populated and vacuum states. The latter is most conveniently and precisely assessed via constant-coverage simulations in rectangular unit cells starting with an initial strip geometry. An important additional feature of our analysis is to demonstrate that generic 2PC is not just a feature of the N = 2 threshold model in the special case without spontaneous particle creation. Rather, it is robust in that it persists for non-zero spontaneous particle creation rate ε > 0 (below a finite critical value ε_c ≈ 0.006). The occurrence of generic 2PC for this simple canonical model exhibiting a discontinuous non-equilibrium phase transition contrasts behavior in thermodynamic systems, which is constrained by the Gibbs phase rule.

Despite the above non-thermodynamic behavior, one finds metastability and nucleation phenomena just above the generic 2PC region which are analogous to thermodynamic systems. More specifically, our analysis indicates the existence of a non-zero effective line tension for interfaces separating coexisting states. This implies no divergence of correlations or criticality at the transition. However, precise characterization of metastability is challenging given the ill-defined nature of the metastable state and the associated spinodal point and of critical droplets.

Finally, we note that our application of homogeneous and heterogeneous mean-field analyses was quite effective in capturing the key features of model behavior. In some respects, these approaches are perhaps more tuned to non-equilibrium models as the associated LDE naturally captures any orientation-dependence of equistability. They cannot capture so well the thermodynamic constraints associated with the Gibbs phase rule that enforces a unique equistability or coexistence point.

We also suggest that generic 2PC is perhaps not such a rare anomaly for non-equilibrium systems displaying discontinuous phase transitions. However, it is clearly not universal. A relevant example is the ZGB model²⁹ apparently exhibiting a discontinuous transition at a single point. However, we shall show in separate work that certain perturbations of the ZGB model can induce generic 2PC. Other non-equilibrium models, which have been cited as indicating the prevalence of discontinuous transitions at a single point, include contact models with large interaction^33,34 or diffusion²⁷ lengths. However, for the large characteristic lengths in these models (relative to the lattice constant), the basic behavior is expected to be captured by continuum mean-field type reaction-diffusion equations of the type analyzed exactly for Schloegl 2nd order model.⁵⁸ Here, there is no orientation-dependence of interface propagation and a unique equistability point is assured (unlike for analogous LDEs). Thus, generic 2PC is not expected.

This work was supported by the U.S. Department of Energy (USDOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences through the Ames Laboratory Chemical Physics program. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358.

Our simulations for the conventional N = 2 threshold model with ε = 0 to determine p_e = p_eq(S = 1) for a diagonal interface utilize a square lattice with L = 512 averaging over 10⁵ MCS. Quite consistent results were obtained for different target concentrations, C_t = 0.30, 0.35, 0.40, and 0.45. Averaging over these simulations gives an estimate p_e = p_eq(S = 1) = 0.351 78(2), where the number in parentheses denotes the standard deviation in the last digit. For comparison, we performed slightly less extensive simulations for slope S = 2 using L = 256 and averaging over 10⁴ MCS to obtain p_eq(S = 2) = 0.3497(1) unambiguously differing from the value for S = 1. We find negligible finite-size effects in the above analysis as in analogous studies for the Durrett N = 2 threshold model with ε = 0.^10,11

Accurate determination of p_f = p_eq(S = ∞ or 0) is more demanding. The above-mentioned anomaly in the model means that a vertical strip of the vacuum state can never be populated. Thus, in a finite system with periodic boundary conditions, one-by-one columns of sites adjacent to the empty strip will eventually become depopulated (as this corresponds to an absorbing state for the individual columns). Thus, simulation requires judicious selection of system size and simulation time, and there are significant finite size effects which must be systematically analyzed for reliable estimation of p_f ≈ 0.316. Detailed discussion of these issues has been presented previously for the Durrett N = 2 threshold model with ε = 0.^10,11

For a more precise analysis of p_f, we have developed a refined version of the CC ensemble simulation algorithm. Now, instead of specifying a target for the total concentration or number of particles in the system, we just specify a target for the partially populated column adjacent to the vacuum strip. Various choices of boundary conditions were also explored. By specifying a target column concentration as C_t = 0.2, one obtains estimates p_f = 0.3130(1), 0.3145(1), and 0.315 35(9) for L = 256, 512, and 1024, respectively. Extrapolating these results using a 1/L variation gives a final estimate of p_f = 0.3161(1).

As an aside, we note that an earlier CC ensemble analy sis for this model in Ref. 23 generated an estimate p_eq = 0.3516. However, this analysis presumably used random initial conditions, which results in complex phase-separated states with a range of interface orientations and thus, an uncontrolled average of p_eq(S)-values.

It is appropriate to mention another simulation strategy for analysis of discontinuous non-equilibrium transitions similar to the CC ensemble. This approach involves simulation of a modified dynamics compatible with that of the contact process of interest, but in a strictly conserved constant particle number (CN) ensemble. Rather than independent spontaneous annihilation and autocatalytic creation events, randomly selected particles are moved with appropriate rates to empty sites which are active in the sense that particles can be autocatalytically created at those sites.^35,59 The equivalence of this constant particle ensemble and conventional simulations was shown in the limit of infinite system size.⁵⁹ 

For the populated homogeneous steady-state in a finite L × L site system with periodic boundary conditions (BCs), fluctuations in the number, M, of particles are described by the general fluctuation-correlation relations, 〈(M − 〈M〉)²〉 = L² χ(p), where the susceptibility χ equals the sum of the pair correlations over all separations. If there is critical behavior at some p_c, then setting Δ = |p − p_c|, one typically writes C (concentration) ∼ Δ^β, ξ (spatial correlation length) ∼ Δ^−ν⊥, and χ ∼ Δ^−γ.³ Then, assuming that χ ∼ C² ξ^d yields the hyperscaling relation γ = d ⋅ ν⊥ − 2β. For discontinuous transitions, one sets β = 0. For criticality at a discontinuous transition, one might write³⁷ C ∼ C₀ + AΔ^β*, but where β^∗ is distinct from β.

For the N = 2 threshold model selecting p < p_e, fluctuations in the populated state should exhibit the feature that χ is independent of L for moderately large L. Indeed, this feature is consistent with the simulation data of da Silva and de Oliveira,^23,35 which indicates that ln [〈(M − 〈M〉)²〉] − 2ln[L] = ln[χ(p)] is roughly independent of L > 40. Of particular significance is the additional feature that χ varies smoothly with p < p_e, and χ does not appear to diverge as p → p_e from below, i.e., γ = 0. Thus, ξ remains finite so ν⊥ = 0, as generally expected for a discontinuous transition. This feature is consistent with our assignment of a non-zero effective line tension at this transition.

For p > p_e, the N = 2 threshold model with ε = 0 will ultimately reside in the fluctuation-free absorbing vacuum state. Behavior for finite systems with periodic BC, in such p-regimes for continuous transitions, is commonly analyzed in terms of scaling for populated quasi-stationary states surviving an initial transient.³ da Silva and de Oliveira^23,35 instead judiciously select “active” BC which assigns permanently populated sites to the edges of the square L × L site simulation cell. Here, behavior for p > p_e is controlled by this selection of BC and not by the general theory described above.

For these active BCs with large L and p > p_e, it is clear that only the interior sites near the corners of the simulation cell will have significant particles population. Thus, behavior for large L is determined by analysis of an infinite square lattice system i ≥ 0 and j ≥ 0, with sites i = 0 and j = 0 permanently populated, and where sites near the origin will be populated. It is clear that as L → ∞, one now has that 〈(M − 〈M〉)²〉 ≈ f(p), where f(p) is independent of L. Simulation results support this conclusion and further, indicate that f(p) ∼ (p − p_e)⁻² for p > p_e allowing estimation of p_e.^23,35 For simulations in a finite L × L system with active BC, it is interesting to note that the magnitude of fluctuations just above p_e exceeds that just below.^23,35 The latter is controlled by the finite correlation length in the populated steady-state. For the former, the system cannot become trapped in the vacuum state due to the active BC and thus, presumably meanders between a near vacuum state and a populated steady-state.