The exact phase diagram for a class of open multispecies asymmetric exclusion processes.

The asymmetric exclusion process is an idealised stochastic model of transport, whose exact solution has given important insight into a general theory of nonequilibrium statistical physics. In this work, we consider a totally asymmetric exclusion process with multiple species of particles on a one-dimensional lattice in contact with reservoirs. We derive the exact nonequilibrium phase diagram for the system in the long time limit. We find two new phenomena in certain regions of the phase diagram: dynamical expulsion when the density of a species becomes zero throughout the system, and dynamical localisation when the density of a species is nonzero only within an interval far from the boundaries. We give a complete explanation of the macroscopic features of the phase diagram using what we call nested fat shocks.

Introduction

The one-dimensional asymmetric simple exclusion process (ASEP) with open boundaries has been of great importance as a model system towards understanding nonequilibrium phenomena. The matrix ansatz, which has since become an important tool[1], was developed to compute the nonequilibrium steady state (NESS) of the totally asymmetric version of the single-species ASEP[2] (TASEP). Using this ansatz, various measurable quantities such as the density and current were calculated in the NESS, from which the exact nonequilibrium phase diagram was derived. While many rigorous results are known in one dimension, there are none for higher dimensional exclusion processes in contact with reservoirs.

The exact calculation of various out-of-equilibrium quantities was also useful in formulating general principles. For example, the diffusion constant was first calculated exactly for the open TASEP[3]. Using the integrability of the ASEP, the spectrum of the transition matrix was computed[4]. The large deviation functional for the density profiles[5] and the large deviation function for the current[6,7] was derived using the matrix ansatz. These calculations helped in the formulation of two general principles for driven diffusive systems; the additivity principle[8] and the macroscopic fluctuation theory[9].

It is natural to extend the ASEP to several species of particles. Multispecies exclusion processes have applications to studies of traffic flow[10], cell motility[11,12], chemotaxis[13], chemical reactions[14] and biological transport in ion-channels[15-17]. The ASEP with second-class particles was first considered with periodic boundary conditions[18]. The first model with two species and open boundaries whose NESS was determined exactly was a model where positive and negative charges moved in the lattice under the influence of an electric field[19]. Much later, a two-species model called the semipermeable TASEP [20,21] was considered, where second-class particles were confined to the lattice. The computation of the NESS was later generalised to the semipermeable ASEP[22]. The NESS of a version with more general boundary conditions was also determined exactly[23]. Later on, the phase diagram for a large class of ASEPs with two species and open boundaries was obtained[24,25]. More recently, classes of integrable ASEPs with multiple species of particles and open boundaries were determined[26].

In this work, we study a multispecies exclusion process on a finite one-dimensional lattice with r species of charges called the mASEP introduced recently[27]. The hopping rates in the bulk are asymmetric and those in the boundary are defined in such a way that there are r + 1 conserved particle numbers. The main results are the following. We obtain the complete (r + 2)-dimensional phase diagram and present formulas for all densities and currents in the thermodynamic limit in all regions of the phase diagram. It will turn out that all the macroscopic features can be explained by a new structure which we call a nested fat shock. To make the presentation self-contained, we review the features for the semipermeable ASEP in Section I of the Supplementary material. We will prove these results in Section II of the Supplementary material by using projections to the semipermeable ASEP.

Results

The mASEP is defined on a one-dimensional lattice of size n, where each site is occupied by exactly one particle of type . The barred particles are negative charges, the unbarred ones are positive charges, and 0’s are vacancies. There are r species of charges, with the total number of particles of charge j being fixed to be for . As a consequence, the number n 0 of vacancies is also fixed, with . More precisely, fix an (r + 1)-tuple of positive integers . The state space consists of all words of length n in the alphabet such that the total number of j′s and ’s is equal to n for and the total number of 0’s is n 0.

The dynamics is the effect of a rightward-pointing electric field. In the bulk, we have the asymmetric hopping rulewhere we think of the barred particles as negative numbers, and set q < 1. On the left and right boundaries, positive charges can only replace and be replaced by their negatively charged partners with rates given by

The mASEP possesses charge-conjugation symmetry in the following sense: interchanging positively and negatively charged particles as well as the rates 1 and q, α and β, and γ and , and changing the direction of motion leaves the mASEP invariant. The model with (r = 1) and is the single-species open ASEP;[2] for arbitrary n 0, this is the semipermeable ASEP[22]. Furthermore, if , this is the semipermeable TASEP[20,21] (see Section I of the Supplementary material). We thus denote the mASEP with as the mTASEP. Results of simulations for the mTASEP with r = 2 are given in Fig. 1.

Figure 1

Plots of the densities of particles 2 (black crosses), 1 (green squares), 0 (red diamonds), (blue circles), and (pink triangles), versus the scaled position for the mTASEP with r = 2 where , and in the regions (a) (), (b) (), (c) (), (d) (), (e) (), (f) the shock line (), and (g) the shock line ().

Phase Diagram

For each integer j between 0 and r, let be the total density of (both positively and negatively charged) particles of species j. We consider the behaviour of the mASEP in the limit and for each j such that the total density of species j particles converges to θ  > 0. Define the quantitiesand the function . Set and for . Then the exact phase diagram is given in Fig. 2.

Figure 2

The nonequilibrium phase diagram of the mASEP with r species of charges. There are 2r + 1 different regions, which are labelled . Each region is characterised by different bulk densities of all particles. The explanation for the nomenclature of the phases in given below. See Tables 1 and 2 for the densities and currents respectively in the NESS in these 2r + 1 regions.

Table 1

The densities of all species of particles in phase , as well as phases and , and the boundary for . Piecewise constant densities correspond to phase separation and piecewise linear profiles correspond to averaging over shocks. The exact formulas can be calculated from the schematic plots in the top row of Fig. 5.

Phase ↓	Densities in the bulk
	Species k	Values of ρ_k,ρ_k¯
0	k = 0	ρ 0 = θ 0
	1 ≤ k ≤ r	ρ k = ρ_k¯=θ_k/2
j	k = 0	ρ 0 piecewise constant
	1 ≤ k ≤ j − 1	ρ_k¯=0 ρ k piecewise constant
	k = j	ρ_j¯=f(b)−Θ_j+1 ρ j piecewise constant
	j + 1 ≤ k ≤ r	ρ k = ρ_k¯=θ_k/2
j¯	k = 0	piecewise constant
	1≤k≤j−1	ρ k = 0ρ_k¯ piecewise constant
	k = j	ρ_j=f(a)−Θ_j+1 ρ_j¯ piecewise constant
	j+1≤k≤r	ρ_k=ρ_k¯=θ_k/2
j−j¯ boundary	k = 0	ρ 0 piecewise linear
	1≤k≤j−1	ρ_k,ρ_k¯ piecewise linear
	j+1≤k≤r	ρ_k=ρ_k¯=θ_k/2

Table 2

The currents of all species of particles in phase , as well as phases and for . All currents are seen to be continuous across the boundary. For the special cases of and , take . Note that J 0 = 0 and .

Phase ↓	Currents
	Species k	Value of J k
0	1≤k<r	(1−q)(Θ_k−Θ_k+1)(1−Θ_k−Θ_k+1)
	k = r	(1−q)Θ_r(1−Θ_r)
j	1≤k≤j−1	0
	k = j	(1−q)(f(b)(1−f(b))−Θ_j+1(1−Θ_j+1))
	j+1≤k<r	(1−q)(Θ_k−Θ_k+1)(1−Θ_k−Θ_k+1)
	k = r	(1−q)Θ_r(1−Θ_r)
j¯	1≤k≤j−1	0
	k = j	(1−q)(f(a)(1−f(a))−Θ_j+1(1−Θ_j+1))
	j+1≤k<r	(1−q)(Θ_k−Θ_k+1)(1−Θ_k−Θ_k+1)
	k = r	(1−q)Θ_r(1−Θ_r)

The nonequilibrium phase diagram of the r-species mASEP in Fig. 2 comprises 2r + 1 phases, , each one of which is characterised by the bulk densities ρ for , as well as the currents J for of all types of particles; we tabulate these in Tables 1 and 2 respectively. Note that there is no current of 0’s and the current of ′s to the left is the same as that of j’s to the right. The mean densities of the j and jump discontinuously across the boundary. By contrast, the mean densities vary continuously along the (and ) boundary. All currents J change continuously across all phase boundaries in Fig. 2. In all phases except , the system shows phase coexistence with a sharp interface separating intervals of different density for some particle type, as we show by the illustrative density profiles in Fig. 1 for the mTASEP with r = 2 and Supplementary Fig. S2 for the mASEP with r = 1. The proofs of the density profiles in Table 1 and currents in Table 2, which explain the phase diagram in Fig. 2, are given in Section II of the Supplementary material.

Figure 5

The top row shows schematic plots of the densities ρ and , for all j, versus the normalised position x illustrating the nested fat shock (a) pinned to the left in region , (b) pinned to the right in region and (c) unpinned in the boundary in (c). The densities are plotted against the rescaled location x. The value of is the height of the region containing particle a at x. The bottom row shows simulation results in multiples of 2000 steps as spatiotemporal plots for the mTASEP with r = 2 and in (a) region (), (b) region (), and (c) the boundary (). The blue, red and green colours represent and 1 particles respectively. The particles of type and are shown in white. See Supplementary Videos 1, 2 and 3 showing the spatiotemporal evolution of the nested fat shocks in simulation runs for (a), (b) and (c) respectively.

We illustrate the density profiles for various regions in the phase diagram for the example of the 2-species mTASEP (i.e. ) in Fig. 1. Therefore, and from (3), and , . There are five phases for this system, corresponding to r = 2, and two relevant phase-boundaries. The currents can be calculated by mean-field type calculations from the densities in each of the regions. More precisely, and . When there is phase-separation, the densities of all species conspire to ensure that the currents are constant across the system (because of particle conservation in the bulk). The value of the currents in each phase can be compared with Table 2.

In phase , all densities are constant, and the densities of oppositely charged particles are equal. Therefore, , and for k = 1, 2. This is seen in Fig. 1(e) and matches with the first row of Table 1.

In phase , the densities of 2’s, and ’s are constant, whereas those of 1’s and 0’s undergo phase separation. As in phase , . Moreover, is also constant throughout the system. In the phase-separated regions, either (forcing ), or from which ρ 0 can be calculated. The density plots can be seen in Fig. 1(b) and match the calculation of the densities in the second row of Table 1 with . The density profiles in can be calculated analogously using charge-conjugation symmetry.

In phase , the only constant densities are given by and . It is not immediately obvious why ’s are excluded from the system, and we give an explanation for this phenomenon of dynamical expulsion in the next section. Particles of species 0, 1 and 2 are phase segregated in three distinct parts. The density of 2’s, , is given by the 1 − β in the rightmost part and in the other two parts. Particles of species 0 exist only in the leftmost part with density , and those of species 1 exist only in the middle part with the same density. In the thermodynamic limit, the middle part is infinitely far away from both boundaries and it is not immediately clear how 1’s can be localised in the bulk. We call this phenomenon dynamical localisation and explain how this occurs in the next section. The densities can be seen in Fig. 1(d) and match the calculation of the densities in the second row of Table 1 with . Again, the profiles in can be calculated using charge-conjugation symmetry.

The nomenclature for the phases can now be explained. Each phase is denoted by the phase-segregated species with largest absolute value. For example, 0’s and 1’s are segregated in phase , 0’s and ’s are segregated in phase , 0’s, 1’s and 2’s are segregated in phase , and 0’s, ’s and ’s are segregated in phase . The sole exception is phase , where all species have constant density.

To understand the density profiles in the and boundaries, we appeal to the nested fat shock construction, which we explain below. Recall that the phase diagram is calculated in the limit where the system size, . We rescale the system by a factor of so that the locations lie in the interval . In the boundary, as shown in Fig. 3, the densities of particles 2 and are constant and equal to . All the particles of type 0 form a ‘bound state’ of fixed width. We call this the nested fat shock (the nesting is of order 1 here). Both ends of the bound state execute a synchronised symmetric random walk with reflecting boundary conditions. As a consequence, , ρ 0 and ρ 1 are piecewise linear after averaging. In particular, is constant towards the left, since the right end of the bound state cannot move all the way to the left, and similarly for . This is shown in Fig. 1(f).

Figure 3

Instantaneous picture of the nested fat shock in the boundary in the rescaled mTASEP with r = 2. Each connected region is labelled with the species of a particle and the height of a region at a given location represents the density of that species at that point.

In the boundary as shown in Fig. 4, none of the densities are constant, and the picture is more complicated. The nested fat shock here consists of the regions containing ’s, 0’s and 's, in that order from left to right. The nesting is of order 2 here. There are four boundaries between the regions , , and , and all of them perform synchronised symmetric random walks in the bulk so that the widths of the regions containing , 0 and 1 is fixed. When one of them touches the boundary the widths of either or 1 can decrease, causing the opposites charged region to increase in size so that the sum of the widths of these two remains constant. The width of the region containing 0 never changes. This behaviour results in the piecewise linear profile shown in Fig. 1(g).

Figure 4

Instantaneous picture of the nested fat shock in the boundary in the rescaled mTASEP with r = 2. Each connected region is labelled with the species of a particle and the height of a region at a given location represents the density of that species at that point.

One can now derive the density profiles in regions , , and from these nested fat shocks. For example, in , one has the same nested fat shock structure as in Fig. 3, but the ends of the bound state containing 0 execute a random walk with negative drift, which ensures that the nested fat shock is pinned to the left. Similarly, the nested fat shock is pinned to the right in . Similarly, the density profiles in and can be calculated by forcing the nested fat shock in Fig. 4 to be pinned to the left and right respectively.

The general structure of the nested fat shock is explained in the next section.

Nested fat shock

All coarse features of the phase diagram in Fig. 2 are explained by the nested fat shock construction. This is a generalisation of the fat shock construction, which explains the phase diagram for the semipermeable TASEP[21] (i.e. the mTASEP with r = 1). Roughly speaking, the fat shock consists of a macroscopic interval of the system, where all the ’s are localised. The ’s form two simultaneous shocks with the ’s and ’s, with a constant macroscopic width. For more details on the fat shock, see Section I of the Supplementary material.

The nested fat shock is a macroscopic interval of the system where, for some j ≥ 0, particles of species are localised in a very specific way.

Particles of species 0 have a nonzero constant density in a subinterval of fixed width inside this interval. Particles of species 1 (resp. ) have a nonzero constant density in a subinterval to the right (resp. left) of the 0’s. Although the widths of the 1 and subintervals may vary, the sum of their widths is fixed. This pattern continues until species j on the right and species on the left. The boundary between any two adjacent subintervals is a shock-front. Depending on which part of the phase diagram the system finds itself in, the nested fat shock can have either negative, positive or zero drift. If the drift is negative, the negatively charged subintervals containing will not exist, and similarly if there is positive drift, the positively charged subintervals will vanish. If there is zero drift, all subintervals will exist and move in a synchronised fashion.

The top row of Fig. 5 shows the structure of the nested fat shock in these three cases in the most general scenario. We give a concrete example of a simulation run of the mTASEP with r = 2 in the bottom row of the figure, which shows the results for phases (a) (negative drift), (b) (positive drift), and (c) the boundary (zero drift). The simulations show the densities of particles of species and 1 only. If the nested fat shock has nonzero drift, it gets pinned to one of the boundaries; it is pinned on the left in (a) and on the right in (b). When the nested fat shock is pinned to the left, it consists of species 0 and 1, and particles of species exit the system from the left. Similarly, when pinned to the right, it consists of species and 0, and particles of species 1 vacate the system from the right. When the nested fat shock has zero drift as in (c), blocks of , 0 and 1 are present, with all shock fronts performing lockstep symmetric random walks.

The precise details for each phase are given below.

Region

Here, only the particles of species 0 participate in the nested fat shock and the width of the shock is larger than the system size. As a result, all densities are constant throughout the system. This explains the densities and currents in region . See the density plot in Fig. 1(e) for the result of simulations.

Regionsand

We give details only for region , since those of can be obtained by analogous arguments.

In region , particles of species participate in the nested fat shock. The velocity of this shock is negative and it gets pinned to the left boundary. However, this is not a stable situation. Initially, particles of species will be replaced by those of species j − 1, which will then move rightwards in the bulk, until they join the subinterval of the nested fat shock occupied by the particles of species j − 1. Once that process is completed, a similar phenomenon will happen with particles of species . This process will continue until all negatively charged particles in the nested fat shock have been replaced by their positive counterparts. In the steady state, we will only see particles of species participating in the shock, which will be pinned to the left.

We point out two new nonequilibrium features of this phase which can be seen in the top row of Fig. 5(a). First, note that species 0 through j − 1 are localised in the interior of the system. Each of these species has undergone phase separation, with one region of non-zero density and the remaining of zero density. The precise locations of the region with non-zero density can be calculated from the values of . What is more interesting is that species 1 through j − 1 are localised away from the boundary. This is somewhat counterintuitive since we have taken the thermodynamic limit and these particles are infinitely far away from the boundary. We call this phenomenon dynamical localisation. Such a phenomenon cannot occur in an equilibrium system. The second new feature is the complete absence of particles of species in the system, i.e. . This is related to the previous phenomenon since these particles can only enter at the expense of the dynamically localised particles. We call this phenomenon dynamical expulsion. In the extreme case of the phase, all the barred particles except are expelled.

The picture in the phase can be derived analogously. See the density plots in Regions and of Fig. 1(a,b,c and d) respectively for the results of simulations. See Fig. 5(a) (resp. (b)) for an illustration of the nested fat shock pinned on the left (resp. right) in the top row and the result of a simulation for r = 2 and (resp. ) in the bottom row. Supplementary Video 1 shows the spatiotemporal evolution of the same fat shock.

Boundary of the region

In the boundary, particles of species participate in the nested fat shock. The velocity of the shock fronts are now zero. Therefore, all these fronts perform a lockstep symmetric random walk in the bulk of the system. All subinterval widths will remain constant until the nested fat shock hits the boundary. When one of the extreme fronts gets pinned to the boundary, the widths of the subinterval containing ’s and ’s can change, but the sum of their widths will remain constant. The other subinterval widths will remain the same. While this front is pinned, the other fronts continue to move synchronously until either another one gets pinned or the one stuck to the boundary gets unpinned. If one more (either j − 2 or ) gets pinned, the same phenomenon will repeat for that species. Note that, for instance, if the nested fat shock gets (temporarily) pinned to the right and then the next shock front also gets pinned, the density of j − 2’s becomes zero. More and more shock fronts can get pinned until the fat shock containing 0’s touches the boundary, at which point the latest front to get pinned can only get unpinned. We thus end up with an instantaneous profile whose schematic is given in the top row of Fig. 5(c). A simulation of the movement of the shock fronts can be seen in the spatiotemporal plot in the bottom row of Fig. 5(c). Supplementary Video 3 shows the spatiotemporal evolution of the same fat shock. The steady state density profile can be obtained by averaging over the uniform shock locations and gives rise to piecewise linear profiles for species . The calculation of these profiles is not difficult, but is tedious and is skipped. The currents have the same values as those in regions and with a = b.

Discussion

In this work, we have found the complete phase diagram for a very general multispecies exclusion process, called the mASEP, in contact with reservoirs. We have found two nonequilibrium phenomena, namely dynamical localisation and dynamical exclusion. Just like the fat shock construction explained the gross features of the semipermeable TASEP in all phases, we find a new object called a nested fat shock which explains the features here. Since the widths of the subintervals in the nested shock are fixed because of the conservation of the total number of j and species, the system is extremely constrained. These constraints play a crucial role in establishing the structure in the various phases. It is an interesting open question as to how the phase diagram will look like with more general boundary rates.

The proofs of our results have appealed in a fundamental way to the colouring argument. It is natural to ask how general this argument is. We are working on a large class of multispecies exclusion processes by systematically exploiting this argument, and we expect new kinds of dynamical structures to appear.

In recent times, multispecies exclusion processes have found applications in physical, chemical and biological systems, as mentioned in the introduction. It would be interesting to see whether experimental realisations can be found for the various phases that we have shown in this work.

Methods

Two different methods were used to obtain the results in this paper. The first method is a mathematically rigorous analysis, leading to the picture of the phase diagram. The second method uses simulations to justify the physical picture of the nested fat shock explained here.

The phase diagram was rigorously obtained by applying a suitable projection on the stochastic process. The details are explained in Section II of the Supplementary material.

The simulations were performed using random sequential updates (RSUs) according to the rules of the dynamics. After reaching steady state, we calculate densities by recording fractional occupation numbers at every site, and currents by noting incoming and outgoing species on the left boundary during subsequent sequences of RSUs. All density profiles show densities for each species averaged over such large sequences. Supplementary Videos showing the evolution of the nested fat shock were obtained by keeping track of individual particles during the RSUs.

Supplementary Information: The exact phase diagram for a class of open multispecies asymmetric exclusion processes

Supplementary Video 1

Supplementary Video 2

Supplementary Video 3