One-dimensional discrete aggregation-fragmentation model
We study an one-dimensional model of aggregation and fragmentation of clusters of particles obeying a stochastic discrete-time kinetics of the generalized Totally Asymmetric Simple Exclusion Process (TASEP) on open chains. Isolated particles and the first particle of a cluster of particles hop one site forward with probability ; when the first particle of a cluster hops, the remaining particles of the same cluster may hop with a modified probability , modeling a special kinematic interaction, or remain in place with probability . The model contains as special cases the TASEP with parallel update () and with sequential backward-ordered update (). These cases have been exactly solved for the stationary states and their properties thoroughly studied. The limiting case of , which corresponds to irreversible aggregation, has been recently studied too. Its phase diagram in the plane of injection () and ejection () probabilities was found to have a different topology.
Here we studied the generic case of attraction . We found that the topology of the phase diagram at changes sharply to the one corresponding to as soon as becomes less than . Then a maximum current phase appears in the square domain and , where are parameter-dependent critical values. By extensive computer simulations, we estimated the singular behavior of , at fixed . The properties of the phase transitions between the three stationary phases at were assessed by computer simulations and random walk theory.
Keywords: non-equilibrium phenomena, one-dimensional processes, generalized TASEP, stationary states, phase transitions
Here we study the stationary properties of particles hopping stochastically in discrete time along finite one-dimensional chains with given boundary conditions at the ends. We believe that the model can be helpful in understanding various types of systems in Nature, such as: kinetics of protein synthesis MacDonald et al. (1968); Zhang et al. (2006), molecular motors on a single track Tipathi and Chowdhury (2008), colloid particles moving in narrow channels Guelich et al. (2007); Kolomeisky (2007); Zilman et al. (2009), and vehicles on a single-lane road Nagle (1996); Chowdhury et al. (2000); Helbing (2001), etc. In the model under consideration, the particles obey the dynamics of the generalized Totally Asymmetric Simple Exclusion Process (TASEP) with backward-sequential update and two hopping probabilities: and . The modified hopping probability describes a kinematic attraction between neighboring particles which hop during the same integer-time moment. In principle, the model admits the study of aggregation-fragmentation phenomena, fluctuations and finite-size effects in nonequilibrium stationary states induced by the boundary conditions.
We remind the reader that the original TASEP was defined as a continuous-time Markov process (random-sequential update in the Monte Carlo simulations) and was solved exactly with the aid of recurrence relations by Domany et al Derrida et al. (1992) for special values of the model parameters, and by Schütz and Domany Schütz and Domany (1993) in the general case. A breakthrough in the methods for solving TASEP on open chains marks the matrix-product representation of the steady-state probability distribution, found in Derrida et al. (1993). Different versions of this approach, known as the Matrix Product Ansatz (MPA), were used also to obtain exact solutions for the stationary states of TASEP and ASEP under several types of discrete-time stochastic dynamics: sublattice-parallel Hinrichsen (1996); Honecker and Peschel (1997), forward-ordered and backward-ordered sequential Rajewski et al. (1996); Rajewski and Schreckenberg (1997), and fully parallel (simultaneous updating of all sites) Evans et al. (1998), de Gier and Nienhuis (1999). The above results show that the only dynamics that can move forward clusters as a whole entity is the backward-ordered one. Then, the probability for translation of a cluster of particles one site to the right is , while such a cluster is broken into two parts with probability . Under the generalized TASEP dynamics, the probability for translation of a cluster of particles one site to the right becomes , and its fragmentation into two parts happens with probability .
We note that the above generalized backward-ordered dynamics was suggested as exactly solvable one by Wölki Wölki (2005), and studied on a ring in Derbyshev et al. (2012, 2015); Aneva and Brankov (2016). The limit case of corresponds to irreversible aggregation, or jam formation in the case of vehicles, suggested by Bunzarova and Pesheva Bunzarova and Pesheva (2017) and further studied in Bunzarova et al. (2017); Brankov et al. (2018).
Ii The model
We consider an open one-dimensional lattice of sites. An occupation number is associated with a site , where , if site is empty and , if site is occupied. The dynamics of the model follows the discrete time backward-sequential rules Rajewski et al. (1996); Rajewski and Schreckenberg (1997). During each discrete moment of time , an update of the configuration of the whole chain with sites, labeled by , takes place in steps, passing through successive updates of the right boundary site , all the pairs of nearest-neighbor sites in the backward order , and, finally, the left boundary site is updated. According to the generalized backward-sequential rules:
1. Each integer time moment (configuration update) starts with the update of the last site of the chain: if site is occupied, the particle at it is removed from the system with probability and stays in place with probability . If the last particle has left the system, then the particle at site takes its place at site with modified probability ; otherwise, it remains immobile with probability .
Next, a particle at site hops to an empty site with probability or , depending on the update of the next nearest neighbor on the right-hand side at the current moment of time.
2. If site is occupied and site is empty at the beginning of the current update, then the particle from site jumps to site with probability and stays immobile with probability . Alternatively, if site is occupied at the beginning of the current moment of time , and became empty after the particle from jumped to the empty , then the particle at site jumps to site with probability and stays immobile with probability .
3. The left boundary condition also depends on the occupation history of the right nearest-neighbor. If site was empty at the beginning of the current update, a particle enters the system with probability or site remains empty with probability . Alternatively, if site was occupied at the beginning of the current moment of time, but became empty under its current update, then a particle enters the chain with probability or the site remains empty with probability , where
We note that the above left boundary condition was introduced first by Hrabák in Hrabák (2014) and, independently, in Bunzarova and Pesheva (2017). It is necessary for consistency with the special cases of backward-ordered sequential update, when and , as well as with the parallel one, when and .
ii.1 Known results in particular cases
Here we summarise the known results about the phase diagrams and the phase transitions between the stationary phases in the particular cases of (the ordinary backward-sequential update) and (irreversible aggregation). The corresponding phase diagrams are shown in Fig. 1.
In Fig. 1(a) the subregions AI (BI) and AII (BII) differ by the shape of the local density profiles. The nonequilibrium phase transitions between phases LD and MC, as well as between HD and MC, are continuous, while the transition between LD and HD is discontinuous, with a finite jump in the local density. The exact critical values are . In our case of , .
In Fig. 1(b) the many-particle phase MP contains a macroscopic number of particles or clusters of size (1) as ; MPI and MPII differ only by the shape of the local density profile, which results from the different type of evolution of the configuration gaps. In the region of MPI, the inequality leads to growing average width of the rightmost gap, hence the profile bends downward near the chain length. In the complementary region MPII, the opposite inequality holds and the rightmost gap is short-living, while the gaps on the left-hand side of it have a critical type of evolution with mean lifetime of the order , see Ref. Brankov et al. (2018). The phase MP+CF is mixed in the sense that the completely filled configurations are perturbed by short living gaps entering the chain from the first site. The configurations of the stationary nonequilibrium phase CF represent a completely filled chain with current . The unusual phase transition, found in Bunzarova and Pesheva (2017), takes place across the boundary between the MPI and CF phases.
Iii The generic case
Firstly, we aim here to analytically approach the question of how the completely filled phase (CF) at , see Fig. 1 (b), is destroyed, when , and transformed into new phases typical for , see Fig. 1(a). One of the methods, developed in Refs. Bunzarova et al. (2017); Brankov et al. (2018), is based on the study of the time evolution of single gaps in different regions of the CF phase.
Secondly, we present results of extensive computer simulations, which suggest the topology of the modified phase diagram, the shift of the triple point under the change of at fixed , and the nature of the phase transitions between the stationary nonequilibrium phases.
iii.1 Time evolution of configuration gaps
We begin with finding out the probability of a single gap appearance under boundary conditions corresponding to the CF phase. Then we consider the first step in the time evolution of the gap width. The problem is rather complicated because the probability of appearance of a gap is position dependent when . In contrast to the case of , here we show that when , the gap width performs a special, position dependent random walk.
Let denote the probability of appearance of a single gap at site in a completely filled configuration with , when . For brevity of notation, we do not show the explicit dependance on the injection and ejection rates of that probability. Since we exclude the appearance of a second gap, under the generalized backward-sequential update we obtain,
Under the assumption , the left boundary condition yields for all , which means that at the beginning of each update. Therefore, we focus on the case when the gap appears at sites . Then, the right edge of the gap, positioned at site can move one site to the right, provided that site is empty. The latter event occurs only if the particle at site leaves the system with probability , and the remaining cluster of particles at sites moves as a whole entity one site to the right, which happens with probability . Thus, the total probability for the particle at the right edge to hop to the right is , and to remain at its place is . On the other hand, the particle at the left edge , being either the rightmost particle of a cluster or isolated, may hop to the right with the position independent probability , and stay immobile with probability . As a result, the gap width increases by one site with probability
decreases by one site with probability
and remains the same with probability
As expected, at these expressions reduce to equalities (4) in Ref. Brankov et al. (2018). As is readily seen, in the alternative case of several coexisting gaps, the above probabilities apply exactly to the rightmost one.
Now we have to average the gap width evolution over the initial probabilities given by (2). The probability normalization factor under the condition of a single gap opened at sites is
Then, the changes in the gap width at the first time step, averaged over all events of gap appearance at sites , become as follows:
The gap width increases by one site with probability
decreases by one site with probability
and remains the same with probability
Notably, at the above results reduce again to equalities (4) in Ref. Brankov et al. (2018).
By comparing the expressions for and , we conclude that on the average a single-site gap will grow after the first time step of its evolution when
When and is fixed, or so that , this condition simplifies to . However, for fixed values of close to 1, will decrease to zero as . For example, in our computer simulations we used: and , which yields and the criterion becomes much stronger, . However, with each time step the right edge of the gap will hop forward by one site with increasing probability , while its left edge may hop to the right with the position independent probability . Thus, the value of will increase and the value of will decrease in the course of time.
Without going into the involved details of the complete gaps evolution, we conjecture that the simple criteria for growing gaps, and for decreasing gaps, holds true. Thus, our expectation, confirmed by computer simulations, is that in the upper region of CF a maximum-current phase will appear. Its local density profile satisfies the inequalities , which follow from the conditions , and the larger probability of gap formation near the end of the chain. In the lower region , the gaps are scarce and short-living, which is indicative of a high-density phase. Again, the left-hand side of the local density profile bends upward to .
Note that in the above consideration . In the case of , the critical values should decrease down to , as .
Iv Phase diagram and phase transitions
We performed Monte Carlo simulations of the generalized TASEP on open chains of mainly and sites. Each run started with relaxation updates and had not less than attempted updates per lattice site. The stationary properties were evaluated by averaging over 100 (quasi)independent runs. The estimated accuracy is for the local particle density and for the current.
First, we compare the behavior of current, , and the local density at the midpoint of the chain, , under modified hopping probabilities and , as a function of the input rate , at chain length , fixed and output rate , see Fig. 2.
We recall that in the standard backward-sequential TASEP with , the exact results in the thermodynamic limit are
in the maximum-current phase
The critical value , shown in Fig. 2 by a vertical green dashed line, corresponds to the transition of the asymptotic behavior of the current near from a parabolic one on the left-hand side of the segment , to a constant value in the MC phase on the right-hand side of it. Due to finite-size effects, the current slightly grows with up to at .
Similarly, at , we see that the phase transition across the segment is continuous too and estimate the critical values , see the vertical red dash-dotted line in Fig. 2. This indicates that the unusual phase transition, found in Bunzarova and Pesheva (2017) at across the boundary becomes a continuous one. Note that in the MC phase the current grows up to and the local density up to 0.7531.
To check the continuity of the phase transition across the segment , we consider in more detail both the - and -dependance of the current on larger lattice and closer to 1, namely and . The results are shown in Fig. 3 (a) and (b).
It is instructive to analytically check the above estimates for the critical values of the injection (ejection) () probabilities, commonly denoted by . To this end we use the continuity condition for the first derivative of the current: . Close to the critical value , we have a parabolic approximation for the current below ,
and a linear approximation above ,
Hence, in the case of a continuous second order phase transition, when , we obtain an estimate for :
In the case of the current as a function of , see Fig. 3(b), the best least-square fit yields:
which leads to the estimate . In spite of the large error interval, this estimate coincides with our former value of .
In the complementary case of the -dependent current, see Fig. 3(a), the best least-square fit yields:
which leads to the estimate . Again, the error bars are rather large, but this estimated value also coincides with our former assessment of . Finally, we assume that within error bars .
Now we focus on the phase transitions taking place by changing across the segment . In Ref. Bunzarova and Pesheva (2017), a mixed MP+CF phase was found at , see Fig. 1(b). This phase is characterized by the nonzero probability of clusters completely filling the chain of sites: changes with from zero at the left phase boundary to at the right boundary , with the CF phase. In Ref. Brankov et al. (2018) the MP+CF was interpreted as a boundary perturbed one. Here we show that as , exponentially decreases to zero, not only in the MP+CF phase but also in the subregion , which at belongs to the CF phase. The results of our computer simulations for the cluster-size distribution at , , and lattice sizes are shown in Fig. 4.
By using a larger series of chain sizes, from to , we obtain that decays exponentially fast with the unlimited increase of ,
The quality of the fit is illustrated in Fig. 5.
Therefore, by continuity arguments, we conclude that in the thermodynamic limit the regions between the left-hand boundary and the right-hand boundary at and , belong to the same phase. The fact that across the coexistence line there occurs a first-order phase transition is seen from the shape of the local density profiles, shown in Fig. 6 at different points in the plane.
By using continuity arguments, we generalize the above results to conjecture a generic phase diagram for with the same topology as in the case of the backward-sequential TASEP, see Fig. 1(a), but with -dependent triple point . In Fig. 7 we exemplify the phase diagram in the particular case of and .
Now we see the similarity in the behavior of the local density profiles in the cases and . In the low-density phase LD = LDILDII the bulk density is less than 1/2, the difference between the LDI and LDII regions is in the right-hand end of the local density profile: in LDI it bends upward, while in LDII it bends downward, similarly to the case of the standard backward-sequential TASEP. This can be readily explained by using the exact relationship , and assuming the same value of the current in the two regions. In the considered particular case of and , the bulk density in the high-density phase is very close to 1, the difference between the regions HDI and HDII being at left-hand side of the profile: in HDI it sharply bends downward, while in HDII it bends upward, again similarly to the case of the standard backward-sequential TASEP.
Additional information can be found in the different gaps evolution in regions LDI and HDI: in both cases , which implies , so that gaps can appear at the first site and evolve throughout the chain; however, in LDI the gaps are wide and long-living, while in HDI they are scarce and very short living, compare Figs. 8(a) and 8 (b). The typical gaps pattern in LDII (HDII) is similar to the one shown for LDI (LDII). These features may explain the large difference in the particle densities in the low-density and high-density phases.
Here we emphasize that the generalized TASEP does not satisfy the particle-hole symmetry inherent to the standard versions of TASEP. For example, the fundamental diagram of the generic model, obtained by Hrabák, see Fig. 6.2 in Hrabák (2014), where , is not symmetric under the replacement . The above mentioned phase diagram was obtained in the thermodynamic limit under periodic boundary conditions, which means that the very bulk dynamics at does not respect the particle whole symmetry. In addition, in our case of open boundaries, the left boundary condition Eq. (1) is not appropriate for introduction of holes from the right chain end. Therefore, it is rather unexpected that the currents in the phases LD and HD, calculated at symmetric points and , are equal: .
Another unexpected result is that as the critical point has a singular behavior at the point . The least-square allometric fit to as is shown in Fig. 9 for , , and .
The analytic expression for the fit reads:
This result supports our conjecture that the standard topology of the phase diagram appears suddenly from the topology, as soon as .
iv.1 Cluster-size probability distribution
To complete the properties of the different stationary phases, we present in Fig. 10 the simulated cluster-size probability distributions at different points of the phase diagram.
As expected, in the low-density subregions LDI and LDII the cluster-size distributions are almost identical and sharply peaked at single particles and clusters containing just a few particles. The least-square fit to that distribution yields the exponential decay
with and . This result is in conformity with the strongly gapped configurations shown in Fig. 8 (a).
More complicated is the behavior of the cluster-size distribution in subregion HDI of the high-density phase: it consists of a sum of two exponentially decaying functions with reciprocal rate parameters and which differ by two orders of magnitude:
where , , , and . Such a behavior could result from a finite-size mixture of typical low-density and high-density configurations. An indication of the finite-size effects is the residual probability of finding completely filled configurations.
In contrast, the behavior of the cluster-size distribution in subregion HDII is very smooth, being described by a single exponential:
with and . Remarkably, the reciprocal rate parameters and have very close values. The relatively slow decay of the probability with the increase of the cluster size is in conformity with the rare and short-living gaps in the typical high-density configurations. However, finite-size effects are also present in this region, as evidenced by the residual probability of completely filled configuration, .
Finally, the cluster-size probability distribution in the maximum-current phase takes somewhat intermediate position with the fit
where and , in the sense that .
We studied the generalized TASEP in the regime of particle attraction between hopping nearest-neighboring particles, when . A central problem to solve was to find how the topology of the phase diagram in the case of irreversible aggregation , see Fig. 1(b), transforms into the topology of the well-known phase diagram of the usual TASEP with backward-ordered sequential update, see Fig. 1(a), when decreases from down to . Based on an incomplete random-walk theory and on extensive Monte Carlo calculations, we conjectured that the above phenomenon takes place sharply, as long as becomes less than 1. The main difference between the phase diagrams for turned out to be the dependence of the critical probabilities on , at fixed . Apart of that, we have shown the similarity of the local density profiles and the current as a function of the injection and ejection probabilities, in the cases and . The main effect of increasing the modified hopping probability turns out to be increase in the values of critical point coordinates, the bulk density and the current. For example, on passing from to , these values grow from
On the ground of our random walk theory and the computer simulations, we have conjectured that the simple criteria , for growing gaps, and , for decreasing gaps, hold true on the average.
An interesting result is the exponential decay to zero of the probability of a complete cluster in the CF phase, when of and the chain length increases unboundedly, see Fig. 9.
In general, we believe that our analysis of the stationary cluster-size probability distribution under different boundary conditions and chain lengths could be helpful for the interpretation of experimental results on cluster growth and fragmentation processes in systems of particles confined to one-dimensional channels.
One of the significant findings of our study is the singular behavior of as a function of at the point , see Fig. 9. Our late colleague V.B. Priezzhev has conjectured that there is some Renormalization Group transformation at the background of that singularity. We will attempt to develop this idea in a future work.
There are still open problems, such as an elaboration of the random walk theory to the extend of yielding both qualitative and quantitative predictions, the analytical derivation of the local density at the chain ends and the value of the current in the different stationary phases, just to mention some.
Acknowledgements.The authors gratefully acknowledge the fruitful discussions with their late colleague and coauthor Professor V.B. Priezzhev, held at an early stage of preparation of the present study. NCP gratefully acknowledges the partial support of grant BG05M2OP001-1.001-0003.
- MacDonald et al. (1968) C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Biopolymers 6, 1 (1968).
- Zhang et al. (2006) F. Zhang, C. L. Liu, and B. R. Hu, J. Neurochemistry 4, 102 (2006).
- Tipathi and Chowdhury (2008) T. Tipathi and D. Chowdhury, Phys. Rev. E 77, 011921 (2008).
- Guelich et al. (2007) P. Guelich, A. Garai, K. Nashinari, A. Schadshneider, and D. Chowdhury, Phys. Rev. E 75, 041905 (2007).
- Kolomeisky (2007) A. B. Kolomeisky, Phys. Rev. Lett. 98, 048105 (2007).
- Zilman et al. (2009) A. Zilman, J. Pearson, and G. Bel, Phys. Rev. Lett. 103, 128103 (2009).
- Nagle (1996) K. Nagle, Phys. Rev. E 53, 4655 (1996).
- Chowdhury et al. (2000) D. Chowdhury, L. Santen, and A. Schadschneider, Phys. Rep. 329, 199 (2000).
- Helbing (2001) D. Helbing, Rev. Mod. Phys. 73, 1067 (2001).
- Derrida et al. (1992) B. Derrida, E. Domany, and D. Mukamel, J. Stat. Phys. 69, 667 (1992).
- Schütz and Domany (1993) G. M. Schütz and E. Domany, J. Stat. Phys. 72, 277 (1993).
- Derrida et al. (1993) B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, J. Phys. A 26, 1493 (1993).
- Hinrichsen (1996) H. Hinrichsen, J. Phys. A 29, 3659 (1996).
- Honecker and Peschel (1997) A. Honecker and I. Peschel, J. Stat. Phys 88, 319 (1997).
- Rajewski et al. (1996) N. Rajewski, A. Schadschneider, and M. Schreckenberg, J. Phys. A 29, L305 (1996).
- Rajewski and Schreckenberg (1997) N. Rajewski and M. Schreckenberg, Physica A 139, 245 (1997).
- Evans et al. (1998) M. R. Evans, N. Rajewsky, and E. R. Speer, J. Stat. Phys. 95, 45 (1998).
- de Gier and Nienhuis (1999) J. de Gier and B. Nienhuis, Phys. Rev. E 59, 4899 (1999).
- Wölki (2005) M. Wölki, Steady States of Discrete Mass Transport Models (Master thesis) (University of Duisburg-Essen, 2005).
- Derbyshev et al. (2012) A. E. Derbyshev, S. S. Poghosyan, A. M. Povolotsky, and V. B. Priezzhev, J. Stat. Mech. 2012, P05014 (2012).
- Derbyshev et al. (2015) A. E. Derbyshev, A. M. Povolotsky, and V. B. Priezzhev, Phys. Rev. E 91, 022125 (2015).
- Aneva and Brankov (2016) B. L. Aneva and J. G. Brankov, Phys. Rev. E 94, 022138 (2016).
- Bunzarova and Pesheva (2017) N. Z. Bunzarova and N. C. Pesheva, Phys. Rev. E 95, 052105 (2017).
- Bunzarova et al. (2017) N. Z. Bunzarova, N. C. Pesheva, V. Priezzhev, and J. G. Brankov, J. Phys: Conf. Series 936, 012026 (2017).
- Brankov et al. (2018) J. G. Brankov, N. Z. Bunzarova, N. C. Pesheva, and V. Priezzhev, Physica A 494, 340 (2018).
- Hrabák (2014) P. Hrabák, Ph.D. Thesis (Czech Technical University, Faculty of Nuclear Sciences and Physical Engineering, Praga, 2014).