Effect of diffusion in simple discontinuous absorbing transition models

Nonequilibrium phase transitions into absorbing states describe several problems, such as wetting phenomena, spreading of diseases, chemical reactions and others [1, 2]. In the last years, a large effort for its characterization including experimental verifications [3–5] and the establishment of distinct universality classes [1, 2, 6] have been undertaken. Generically, continuous phase transitions into an absorbing state for systems without conservation laws nor extra symmetries fall into in the directed percolation (DP) universality class [2, 7, 8]. The best example of this category is probably the contact process (CP) [7]. It is defined on a given d-dimensional lattice and the dynamics comprehends the creation in the presence of at least one adjacent particle and spontaneous annihilation. Since the particle creation requires the presence of adjacent particles, a configuration devoided of species is absorbing.

Conversely, discontinuous absorbing transitions also appear in different systems [9–13], but comparatively they have received less attention than the continuous ones. Recently, they have attracted interest for the search of minimal ingredients for their occurrence, being one possibility the so called 'restrictive' (threshold) contact processes (CPs). These models are variants of the second Schlögl model, in which the particle creation is similar to the usual CP, but instead it requires a minimal neighborhood larger than 1 particle for creating a new species [14–18]. Different studies have stated that the phase transition, induced by this mild change (with respect to the usual CP), remains unaffected under the inclusion of distinct creation [14–18] and annihilation rules [18] as well as for different lattice topologies [19]. In some specific cases [14, 15], in which the particle creation occurs only if one has at least two adjacent diagonal pairs of particles, the phase transition is characterized by a generic two-phase coexistence and exhibits an interface orientational dependence at the transition point. On the other hand, when the transition rates depend only on the number of nearest neighbors particles (and not their orientations) the discontinuous transitions take place at a single point [16–18].

Despite the apparent robustness of first-order transitions for the above mentioned restrictive examples, the effect of some (relevant) dynamics has been so far unexplored in the present context. These dynamics, such as spatial disorder and particle diffusion, can cause drastic changes in continuous phase transitions [2]. One of the few available studies shows that spatial disorder suppresses the phase coexistence, giving rise to a continuous transition belonging to a new universality class [20]. A similar scenario of scarce results also holds for the outcome of diffusion. Very recently, a stochastic differential equation (such as a Langevin equation) reported that different diffusion strengths lead to opposite findings (in two dimensions). Whenever the transition is discontinuous for larger diffusion values, limited rates suppress it, giving rise to a critical phase transition belonging to the DP universality class [21]. With these ideas in mind, we give a further step by tackling the influence of diffusion in lattice systems presenting discontinuous absorbing phase transitions. [16, 17]. Our study aims to answer three fundamental points: (i) what is the effect of strong and limited diffusion in these cases? (ii) does it suppress the discontinuous transition? (iii) How does our results compare with those obtained from the coarse-grained description in [21]? In other words, are there differences between lattice models and Langevin equations?

We consider two lattice models and three representative values of diffusion, in order to exemplify the low, intermediate and large regimes. Models will be studied via MFT and distinct kinds of numerical simulations (explained further). Results suggest that the discontinuous transition is maintained in both models for all diffusion values. Also, a finite-size scaling behavior similar to discontinuous transitions studied in [18] has been found. Since there is no theory for the nonequilibrium case, our results can shed light over a general finite-size scaling for them.

This paper is organized as follows: section 2 presents the models and mean-field analysis. Section 3 shows the numerical results and finally conclusions are drawn in section 4.

Let us consider systems of interacting particles placed on a square lattice of linear size L. Each site has an occupation variable ${{\eta}_{i}}$ that assumes the value 0(1) whenever sites are empty (occupied). The model A is defined by the following interaction rules: particles are annihilated with rate α and are created in empty sites only if their number of nearest neighbor particles nn is larger than 1 ($nn\geqslant 2$), with rate nn/4 [16]. There is no particle creation if $nn\leqslant 1$. The model B is similar to the model A, but the particle creation rate always reads 1, provided $nn\geqslant 2$ [17]. Thus, whenever particles are created with rate proportional to the number of their nearest neighbors for the model A, it is independent on nn in the model B. Besides the above creation-annihilation dynamics, each particle also hops to one of its nearest neighbor sites with probability D, provided it is empty. In the regime of low annihilation parameters, the system exhibits indefinite activity in which particles are continuously created and destroyed. In contrast, for larger α's, the system is constrained in the absorbing phase. The phase transition separates above regimes at a transition point ${{\alpha}_{0}}$. In the absence of diffusion, the phase transitions for both models A and B are discontinuous and occur at ${{\alpha}_{0}}=0.2007(6)$ and 0.352(1), respectively [16, 17].

The first inspection of the effect of diffusion can be achieved by performing mean-field analysis. The starting point is to write down the time evolution of relevant quantities from the interaction rules and truncating the associated probabilities at a given level. Since the diffusion conserves the number of particles, it is required to take into account at least correlations of two sites and hence two equations are needed. Designating the symbols $\bullet$ and ∘ to represent occupied (${{\eta}_{i}}=1$) and (${{\eta}_{i}}=0$) empty sites, the system density ρ corresponds to the one-site probability $\rho =P\left(\bullet\right)$. Another quantity to be considered here is the two-site correlation given by $u=P\left(\circ \bullet\right)$. From the above model rules, it follows that

$$\begin{eqnarray}&&\frac{\text{d}\rho}{\text{d}t}=2P\left(\circ \bullet\bullet\circ \circ \right)+P\left(\circ \bullet\circ \bullet\circ \right)+3P\left(\circ \bullet\circ \bullet\bullet\right)+P\left(\circ \bullet\bullet\bullet\bullet\right)-\alpha P\left(\bullet\right),\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}\begin{array}{*{20}{l}}\frac{\text{d}u}{\text{d}t}&=&(1-D)[-\frac{3}{2}P(\circ \bullet\circ \bullet\bullet)-P(\circ \bullet\bullet\bullet\bullet) \\ &&-\alpha \,P(\circ \bullet)+\alpha \,P(\bullet\bullet)]+D[6P(\bullet\circ \bullet)-6P(\circ \bullet\bullet)],\end{array}\end{eqnarray} \tag{ 2 }$$

for the model A and

$$\begin{eqnarray}&&\frac{\text{d}\rho}{\text{d}t}=4P\left(\circ \bullet\bullet\circ \circ \right)+2P\left(\circ \bullet\circ \bullet\circ \right)+4P\left(\circ \bullet\circ \bullet\bullet\right)+P\left(\circ \bullet\bullet\bullet\bullet\right)-\alpha P\left(\bullet\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}\frac{\text{d}u}{\text{d}t}&=&(1-D)[-2P(\circ \bullet\circ \bullet\bullet)-P(\circ \bullet\bullet\bullet\bullet) \\ &&-\alpha \,P(\circ \bullet)+\alpha \,P(\bullet\bullet)]+D[6P(\bullet\circ \bullet)-6P(\circ \bullet\bullet)],\end{array}\end{eqnarray} \tag{ 4 }$$

for the model B. Here the symbol $P\left({{\eta}_{0}},{{\eta}_{1}},{{\eta}_{2}},{{\eta}_{3}},{{\eta}_{4}}\right)$ denotes the probability of finding the central site in the state ${{\eta}_{0}}$ and its four nearest neighbors in the states ${{\eta}_{1}}$, ${{\eta}_{2}}$, ${{\eta}_{3}}$ and ${{\eta}_{4}}$.

The pair mean-field approximation consists of rewriting the n-site probabilities (n  >  2) as products of two-site in such a way that

$$\begin{eqnarray}&&P\left({{\eta}_{0}},{{\eta}_{1}},...,{{\eta}_{n-1}}\right)\simeq \frac{P\left({{\eta}_{0}},{{\eta}_{1}}\right)P\left({{\eta}_{0}},{{\eta}_{2}}\right)...P\left({{\eta}_{0}},{{\eta}_{n-1}}\right)}{P{{\left({{\eta}_{0}}\right)}^{n-2}}}.\end{eqnarray} \tag{ 5 }$$

From this approximation, the above equations read

$$\begin{eqnarray}&&\frac{\text{d}\rho}{\text{d}t}=\frac{3{{u}^{2}}}{\left(1-\rho \right)}-\frac{3{{u}^{3}}}{{{\left(1-\rho \right)}^{2}}}+\frac{{{u}^{4}}}{{{\left(1-\rho \right)}^{3}}}-\alpha \rho ,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}\frac{\text{d}u}{\text{d}t}&=&(1-D)\left[-\frac{3{{u}^{3}}}{2{{(1-\rho )}^{2}}}+\frac{{{u}^{4}}}{2{{(1-\rho )}^{3}}}-2\alpha \,u+\alpha \rho \right] \\ &&+\,D\left[6u-\frac{6{{u}^{2}}}{\rho (1-\rho )}\right],\end{array}\end{eqnarray} \tag{ 7 }$$

for the model A and

$$\begin{eqnarray}&&\frac{\text{d}\rho}{\text{d}t}=\frac{6{{u}^{2}}}{\left(1-\rho \right)}-\frac{8{{u}^{3}}}{{{\left(1-\rho \right)}^{2}}}+\frac{3{{u}^{4}}}{{{\left(1-\rho \right)}^{3}}}-\alpha \rho ,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}\frac{\text{d}u}{\text{d}t}&=&(1-D)\left[-\frac{2{{u}^{3}}}{{{\left(1-\rho \right)}^{2}}}+\frac{{{u}^{4}}}{{{\left(1-\rho \right)}^{3}}}-2\alpha \,u+\alpha \rho \right] \\ &&+D\left[6u-\frac{6{{u}^{2}}}{\rho \left(1-\rho \right)}\right],\end{array}\end{eqnarray} \tag{ 9 }$$

for the model B.

The steady solutions are obtained by taking $\frac{\text{d}\rho}{\text{d}t}=\frac{\text{d}u}{\text{d}t}=0$ implying that for locating the transition point and the order of transition it is required to solve a system of two coupled equations for a given set of parameters (α, D). Although alternative treatments have been considered [22], here we shall identify the order of transitions by inspecting the dependence of ρ versus α. In similarity with equilibrium transitions, the existence of a spinodal behavior (with ρ increasing by raising α) signals a discontinuous transition. Despite this, no analogous treatments similar to the 'Maxwell construction' are available. For instance, the coexistence points have been estimated by the maximum value of α and the spinodal behavior replaced by a jump in ρ.

Figures 1(a) and (b) show the phase diagram for distinct diffusion rates. For all values of D, mean-field results predict discontinuous transitions separating the absorbing and the active phases. However, as D increases, the transition point moves to larger values of ${{\alpha}_{0}}$ and the active phase becomes less dense. For example, for D  =  0.1 and D  =  0.9, the active phases have densities ${{\rho}_{\text{ac}}}=0.445$ and 0.371 (model A) and ${{\rho}_{\text{ac}}}=0.419$ and 0.324 (model B), respectively. For $D\to 1$, one recovers the limit of fully uncorrelated particles, so that $P\left({{\eta}_{0}},{{\eta}_{1}},...,{{\eta}_{n-1}}\right)=P\left({{\eta}_{0}}\right)P\left({{\eta}_{1}}\right)...P\left({{\eta}_{n-1}}\right)$. In this limit, equations (1) and (3) become equivalent to those obtained from the one-site MFT given by $\frac{\text{d}\rho}{\text{d}t}={{\rho}^{2}}\left(1-\rho \right)\left[3-3\rho +{{\rho}^{2}}\right]-\alpha \rho$ and $\frac{\text{d}\rho}{\text{d}t}={{\rho}^{2}}\left(1-\rho \right)\left[6-8\rho +3{{\rho}^{2}}\right]-\alpha \rho$, for the models A and B, respectively. For such regime, transition points take place at ${{\alpha}_{0}}=0.4724$ with active phase density ${{\rho}_{\text{ac}}}=0.370$ (model A) and ${{\alpha}_{0}}=0.8154$ with ${{\rho}_{\text{ac}}}=0.322$ (model B). Thus, despite the increase of diffusion displaces particles, MFT predicts a discontinuous transition, regardless the strength of the diffusion rate. As a difference between models, active phases are somewhat more dense for model A than in model B.

Zoom In Zoom Out Reset image size

Numerical simulations have been performed for distinct system sizes of a square lattice with periodic boundary conditions. Due to the lack of a general theory for discontinuous absorbing transitions, including the absence of a finite-size scaling (FSS) theory and unknown spreading experiments behavior for $d\geqslant 2$, two sort of analysis will be presented. In the former, we study the time decay of the density ρ starting from a fully occupied initial condition for distinct independent runs. As for critical and discontinuous phase transitions, for small α the density ρ converges to a definite value indicating endless activity, in which particles are continuously created and annihilated. On the contrary, for sufficiently large α's, the system density ρ vanishes exponentially toward a complete particle extinction. The 'coexistence' point ${{\alpha}_{0}}$ is the separatrix between above regimes, whereas at the critical point ${{\alpha}_{c}}$ the density ρ vanishes algebraically following a power-law behavior $\rho \sim {{t}^{-\theta}}$, with θ its associated critical exponent. For the DP universality class it reads $\theta =0.4505(10)$ [26]. Thus, the difference of above behaviors will be used to identify the order of phase transition.

Second, the behavior of typical quantities in the steady regime is investigated. For instance, we apply the models dynamics together with the quasi-steady method [23]. Briefly, it consists of storing a list of M active configurations (here we store M  =  2000–3000 configurations) and whenever the system falls into the absorbing state a configuration is randomly extracted from the list. The ensemble of stored configurations is continuously updated, where for each MC step a configuration belonging to the list is replaced with probability $\tilde{p}$ (typically one takes $\tilde{p}=0.01$) by the actual system configuration, provided it is not absorbing. Discontinuous transitions are typically signed by bimodal probability distribution ${{P}_{\rho}}$ (characterizing the 'coexistence' between absorbing and active phases) and for finite systems, a peak in the order-parameter variance $\chi ={{L}^{2}}[\langle {{\rho}^{2}}\rangle -{{\langle \rho \rangle}^{2}}]$ is expected close to the transition point. For equilibrium systems, the maximum of χ and other quantities scale with the system volume and its position ${{\alpha}_{L}}$ obeys the asymptotic relation ${{\alpha}_{L}}={{\alpha}_{0}}-c/{{L}^{2}}$ [24], being ${{\alpha}_{0}}$ the transition point in the thermodynamic limit and c a constant. Although the finite size scaling properties are unknown for nonequilibrium systems, results for some first-order transitions into an absorbing state have shown a similar scaling than the equilibrium case (with the system volume) [18, 27].

Results for the model A and three representative diffusion rates (D  =  0.1, D  =  0.5 (not shown) and 0.9) are summarized in figures 2 and 3. In all cases, panels (a) show that there is a threshold point ${{\alpha}_{0}}$ separating an active state (signed by the convergence to a definite value of ρ) from an exponential vanishing of ρ. Such values increase by raising D, yielding at ${{\alpha}_{0}}\sim 0.2590$, 0.360 and 0.436 for D  =  0.1, 0.5 and 0.9, respectively. Since no power-law behavior is presented, such analysis provides the first evidence of a discontinuous transition for all diffusion strengths. In order to confirm this, panels (b) show the probability distribution ${{P}_{\rho}}$ for distinct system sizes. In all cases ${{P}_{\rho}}$ presents a bimodal shape, whose position of equal peaks deviate mildly as L increases. Also, the peaks corresponding to the active and absorbing phases present distinct dependencies on L. Whereas active phase densities ${{\rho}_{\text{ac}}}$ converge to well defined values, ${{\rho}_{\text{ab}}}$ vanishes with 1/L²(insets). For D  =  0.1, 0.5 and 0.9 the ${{\rho}_{\text{ac}}}$'s converge to 0.603(1), 0.501(2) and 0.434(2), respectively. These features are similar to the restrictive models studied in [18]. For D  =  0.99 (not shown) a bimodal probability distribution with active phase centered at ${{\rho}_{\text{ac}}}=0.405(1)$ is observed. In the third analysis, the behavior of the system density ρ and its variance χ is presented. Note that ρ vanishes in a short range of α followed by a peak of the variance χ (panels (c) and their insets). For each D, the positions ${{\alpha}_{L}}$'s, in which χ presents a maximum, scale with L⁻² whereas the maximum of χ increases with L² (panels (d) and insets). These features are similar to equilibrium discontinuous phase transitions [24, 25] and with the scalings in [18]. From this behavior, we obtain the values ${{\alpha}_{0}}=0.2600(1)$ (D  =  0.1), 0.3624(2) (D  =  0.5) and 0.442(1) (D  =  0.9), which agree with previous estimates, obtained from the time decay of ρ. Small discrepancies between estimates can be attributed to the lattice simulated for the time decay be finite, due to the uncertainties in the position of peaks or both. Thus, steady analysis reinforces above conclusions concerning the phase transition be discontinuous regardless the diffusion rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we extended above analysis for model B, with results summarized in figures 4 and 5 for D  =  0.1 and D  =  0.9, respectively. As for the model A, both values of D show that the separatrix between active and absorbing regimes are signed by the absence of a power-law behavior and thus the transitions seem to be discontinuous. In particular, the separatrix points (panels (a)) yield close to ${{\alpha}_{0}}\sim 0.5015$ (D  =  0.1) and 0.770 (D  =  0.9). In order to confirm, we also examine the probability distribution as well as the behavior of ρ and χ. In the former, ${{P}_{\rho}}$ also presents bimodal shapes exhibiting two well defined peaks (panels (b)). Whenever ${{\rho}_{\text{ac}}}$'s (insets) saturate as L increases, the ${{\rho}_{\text{ab}}}$'s vanish as 1/L². In addition, active phase also becomes less dense as D increases, reading ${{\rho}_{\text{ac}}}=0.482(1)$ and 0.372(1) for D  =  0.1 and 0.9, respectively. As for the model A and those studied in [18] the positions of peaks ${{\alpha}_{L}}$'s (in which χ presents a peak) as well the peaks also scale with L⁻² and L² (panels (d)), respectively. From this scaling behavior, we obtain the values ${{\alpha}_{0}}=0.5027(2)$ and 0.7844(2) for D  =  0.1 and D  =  0.9, respectively, in agreement with above estimates (from the time decay of ρ).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thus, in similarity with MFT, numerical simulations suggest that the diffusion does not change the order of phase transition, although clusters become less compact as D increases. Similar conclusions are found for extremely low diffusion strengths. For example, for the model A with D  =  0.01 and 0.05, all above features are verified at $\alpha \sim 0.215$ and 0.239, respectively. Moreover, results of both models contrast partially with those obtained by Villa-Martin et al [21] in the regime of low diffusion rates, although they agree qualitatively in the limit of intermediate and large diffusion regimes. A possible explanation for such differences is presented as follows: As shown in [21], a coarse-grained description of a discontinuous absorbing transition is the differential equation ${{\partial}_{t}}\rho =-\alpha \rho -b{{\rho}^{2}}-c{{\rho}^{3}}+D{{\nabla}^{2}}\rho +\eta (x,t)$, where $D{{\nabla}^{2}}\rho$ corresponds to the diffusion term and $\eta (x,t)$ is the (white) noise. The parameter b is responsible for the density discontinuity, since their signs b  <  0 (>0) provide two (one) stable solutions. On the other hand, the parameter c is required to be positive (c  >  0) for ensuring finite densities. Following the ideas from [20], the deterministic part of above equation ($\alpha \rho -b{{\rho}^{2}}-c{{\rho}^{3}}$) can be obtained for example by taking a fully connected lattice with the second Schlögl transition rates ${{W}^{-}}\left(\rho \to \rho -1/{{L}^{2}}\right)=\alpha \rho$ and ${{W}^{+}}\left(\rho \to \rho +1/{{L}^{2}}\right)={{\rho}^{2}}\left(1-\rho \right)$, where at each time instant the system density ρ changes by a factor $\pm$1/L². A similar reasoning can be extended for the present studied models, but with different W ⁺ 's. In particular, they read ${{W}^{+}}\left(\rho \to \rho +1/{{L}^{2}}\right)={{\rho}^{2}}\left(1-\rho \right)\left(3-3\rho +{{\rho}^{2}}\right)$ and ${{W}^{+}}\left(\rho \to \rho +1/{{L}^{2}}\right)={{\rho}^{2}}\left(1-\rho \right)\left(6-8\rho +3{{\rho}^{2}}\right)$ for the models A and B, respectively and leading to the (deterministic) terms ${{\partial}_{t}}\rho =-\alpha \rho +3{{\rho}^{2}}-6{{\rho}^{3}}+4{{\rho}^{4}}-{{\rho}^{5}}$ (model A) and ${{\partial}_{t}}\rho =-\alpha \rho +6{{\rho}^{2}}-14{{\rho}^{3}}+11{{\rho}^{4}}-3{{\rho}^{5}}$ (model B). Note that the above equations are in agreement with those obtained from the one-site MFT. Thus, different coarse grained descriptions can explain the difference between results in the regime of low diffusion rates. It is worth remarking that further studies are still required to confirm above points.

Despite the similarities between models A and B, some differences are clearly observed. In particular, as predicted by the MFT, compact clusters are somewhat less dense for model B than for model A. Also, the dependence between D and α in both models is somewhat different.

Extending the aforementioned analysis for distinct values of D, we obtain the phase diagram shown in figure 6. For both models, in the limit $D\to 1$ the transition points approach their values predicted by the MFT.

Zoom In Zoom Out Reset image size

To sum up, we investigated the influence of diffusion in two simple models presenting discontinuous absorbing phase transitions. They have been studied via mean-field analysis and distinct numerical simulations. All results suggest that transitions are discontinuous irrespective the strength of the diffusion, in contrast to recent findings in which limited diffusion induces a critical transition [21]. Thus our results indicate not only an additional feature of diffusion, but also the possible differences between lattice model and Langevin approaches. It is worth emphasizing that the study of other models variants (taking into account the inclusion of distinct annihilation rates) are required for checking whether the disagreement between approaches is also verified in other cases. Other remarkable result is that all transitions presented a finite-size scaling (with the system volume) similar to the discontinuous transitions studied in [18]. Although further investigations are still required, the obtained results reinforces the possibility of a general scaling for nonequilibrium phase transitions.

The authors wish to thank Brazilian scientific agencies CNPq, INCT-FCx and FAPESP for the financial support.