Multi species asymmetric simple exclusion process with impurity activated flips

We obtain an exact matrix product steady state for a class of multi species asymmetric simple exclusion process with impurities, under periodic boundary condition. Alongside the usual hopping dynamics, an additional flip dynamics is activated only in the presence of impurities. Although the microscopic dynamics renders the system to be non-ergodic, exact analytical results for observables are obtained in steady states for a specific class of initial configurations. Interesting physical features including negative differential mobility and transition of correlations from negative to positive with changing vacancy density, have been observed. We discuss plausible connections of this exactly solvable model with multi lane asymmetric simple exclusion processes as well as enzymatic chemical reactions.

Apart from its extensive success in modeling numerous real-world phenomena, ASEP and its variations have been instrumental in understanding the mathematical structures and physical characteristics of generic non-equilibrium steady states and dynamics [2,12,[17][18][19][20].In particular, the exact steady states of the totally asymmetric simple exclusion process (TASEP) [12,21] and the general ASEP [22] with open boundary conditions, have been obtained using matrix product ansatz.Except for the infinite dimensional representations [21,22], interesting and especially useful finite dimensional matrix representations have been achieved for corresponding quadratic algebra for certain conditions on the transition rates [23,24].The matrix product ansatz has been extremely effective in deriving the non-equilibrium steady states of several generalizations of ASEP including two species [25] and multi-species processes [26], see Ref. [27] for a detailed review.In fact, the stationary state for the multi species TASEP has been solved remarkably by a different method of multiline queuing process [28], which is explored further in terms of combinatorial R in crystal base theory [29].Several two point and three point correlations have been studied analytically in the multi species TASEP [30] and inhomogeneous multi species TASEP with species dependent rates have been analyzed [31,32].The multi species ASEP has also been investigated with integrable open boundary conditions [33] and matrix product solutions are found [34].Due to the connection of TASEP to integrable spin chains [35], the algebraic Bethe Ansatz has been applied to study the dynamics of TASEP [36] and ASEP [37] with open boundaries.Interestingly, ASEP belongs to Kardar-Parisi-Zhang universality class [38] with dynamic exponent 3  2 [39,40].With the aid of Monte Carlo simulations and several improved versions of mean-field theories, TASEP has also been generalized to non-conserved dynamics [41], two-lane [42][43][44] and multi-lane [45][46][47][48] models relating to traffic flow and complex networks [49].
It is quite natural to expect the presence of multiple species of particles with a variety of microscopic dynamics in a system in general.Often due to the distinction between the dynamics of different species, some species are referred to as impurities and give rise to fascinating physical and mathematical structures.For example, the presence of a single impurity which hops with a different rate and allows overtaking of ordinary particles in the TASEP on a periodic lattice, leads to a matrix product steady state with six distinct phases including the creation of a shock in one of the phases [50].This impurity model has been generalized by considering bidirectional asymmetric hopping of the ordinary particles [51] which allows for finite dimensional matrix representations in certain regions of the parameter space, in comparison to the infinite dimensional representation in Ref. [50].A phase transition arising from the motion of the single impurity in the direction opposite to the ordinary particles has also been observed [52].The long time limit behavior of the TASEP with a single impurity has been solved using the Bethe Ansatz [53] and the diffusion constant of the impurity has been calculated from both the Bethe Ansatz [53] and the matrix product ansatz [54].Remarkably, a disordered ASEP with species dependent hop rates, has been shown to exhibit Bose-Einstein condensation [55].A variation of the ASEP with ordinary particles and many impurities has been considered in Ref. [56], where the impurities are not allowed to hop to the vacant neighbors but they can exchange positions with ordinary particles.Interestingly, such an impurity model in [56] possesses a different scaling exponent 5  2 in comparison to the usual KPZ exponent 3  2 for the ASEP without impurities [39,40].Other than the disorders or impurities associated with the particles themselves, there are many exciting studies with position dependent or site-wise disorders for ASEP [57][58][59][60][61][62][63].
In this article, we study a class of the multi-species (I = 1, 2, . . ., µ) ASEP in the presence of impurities, under periodic boundary conditions.In addition to the usual hopping of particles to vacant sites in ASEP, we consider flips of different species among each other (e.g.species I transforming to species J and vice versa).Importantly, these flip processes are initiated only in the presence of a special type of particles (as nearest neighbors) that we denote impurities.These impurities activate the flip processes.Therefore, we name this non-equilibrium stochastic process to be multi species asymmetric simple exclusion process with impurity activated flips (µ-ASEP-IAF).Note that the total number of impurities, along with that of the vacancies, remain conserved in the µ-ASEP-IAF.Specifically, we emphasize that the flip processes between two non-conserved species (I, J) do not occur through the interaction with any non-conserved species K (= 1, 2, . . ., µ) at the nearest neighbors.Thus, the microscopic dynamics considered here is different from previously studied models like TASEP with internal degrees of freedom [64] and multi-species reaction-diffusion processes [65,66].Notably, the distinction in the microscopic dynamics also makes µ-ASEP-IAF non-ergodic in nature in contrast to the ergodic models [64][65][66].We should mention that the non-ergodicity of exactly solvable models is related to undecidability of thermalization in integrable models [67].
The motivations for studying the µ-ASEP-IAF are as follows.(i) We aim to obtain an exact non-equilibrium steady state of the µ-ASEP-IAF under periodic boundary condition, so that it would be an important addition to the category of exact solvable models in disordered sys-tems.(ii) The µ-ASEP-IAF being non-ergodic, it would be interesting to derive exact analytical expressions for partition function and observables for suitable choice of initial configurations and compare the corresponding steady state results with that of a random initial configuration.(iii) The µ-TASEP-IAF can be mapped to multi-lane TASEP which is a basic model for multi-lane traffic flow.Different species of particles in µ-ASEP-IAF play the roles of particles in different lanes of multi-lane TASEP and the impurities in µ-ASEP-IAF act as bridges between lanes that allow particles to exchange lanes in multi lane TASEP.See Appendix A for details.(iv) Considering the conserved impurities as enzymes (E) and different non-conserved species as substrates (S) and products (P), the flip process of µ-ASEP-IAF can be thought as an enzymatic chemical reaction like S + E → P + E, which is a crude approximation of the Michaelis-Menten reaction scheme S + E ⇋ S E → P + E [68][69][70].See Appendix C for details.Notably, both the mappings in (iii) and (iv) would not be possible if the impurities could also flip to other species.
Below we briefly summarize our main results.(i) We find that the steady states of µ-TASEP-IAF (totally asymmetric hopping) and µ-ASEP-IAF (bidirectional hopping) under periodic boundary conditions, can be obtained exactly as matrix product states, where distinct matrices represent different components (each species, impurity, vacancy) of the system.We provide explicit finite dimensional matrix representations for the totally asymmetric case, whereas the matrices for the general asymmetric case are found to be infinite dimensional .(ii) For a specific choice of initial configuration, we could analytically calculate the partition function in the sector of allowed configurations in the steady state and consequently the observables of interests (average densities of non-conserved species, currents and spatial correlations).The analytical results are in agreement with Monte Carlo simulations.For a fixed set of input parameters, we show considerable quantitative deviations between steady state observable values for different initial configurations, establishing the initial configuration dependence or non-ergodicity of the dynamics.(iii) Two-point nearest neighbor correlations exhibit interesting non-trivial behaviors.Particularly, with the variation of the vacancy density, we observe characteristics like certain correlations changing signs i.e. varying from negative to positive with some intermediate zero correlation point, and, non-monotonic behavior with both local maximum and local minimum.(iv) We find negative differential mobility in µ-ASEP-IAF.For special choices of hopping rates, both the drift current and flip current decrease with increasing bias giving rise to negative differential mobility.
The article is organized as follows.In Sec. 2 we describe the µ-TASEP-IAF in details and show that the steady state can be achieved using matrix product ansatz.The analytical calculation of partition function starting from a suitably chosen initial configuration is presented in Sec. 3. We discuss the behaviors of observables like species densities, drift current, flip current and spatial correlations from both analytical calculation and Monte Carlo simulations with variation of input parameters in Sec. 4. The µ-TASEP-IAF is generalized to µ-ASEP-IAF with bidirectional motions of the species in Sec. 5, where we show the corresponding matrix product states and discuss the negative differential mobility of particles.In Sec. 6, we summarize the results with future directions.We discuss the mapping between µ-TASEP-IAF and multi lane TASEP in Appendix A. A variation of µ-TASEP-IAF that comes up with better features in connection to traffic in multi-lane problems is discussed in Appendix B. The connections between µ-TASEP-IAF and enzymatic chemical reactions are briefly presented in Appendix C. In Appendix D, we provide explicit solutions for the fugacity, in the grand canonical ensemble, for some special choices of the input parameters.The block-diagonal structure of the transition rate matrix dictating the transitions between configurations in the configuration space, is presented in Appendix E.

Microscopic dynamics
Let us consider a system of µ different species of particles and impurities on a one dimensional periodic lattice with L sites i = 1, 2, . . ., L. Each site can either be vacant or it can be occupied by only one particle of any of the species I = 1, 2, . . ., µ or by an impurity.All the particles obey hardcore exclusion.The impurity and the vacancy are denoted by + and 0, respectively.The system evolves according to the microscopic dynamics given below, drift (species) : According to the dynamics of the 2-TASEP-IAF in Eq. (1), a particle of species I can hop to its right nearest neighbor with rate p I if the target site is vacant.The impurity (+) hopping rate is ε.If the right neighbor of a particle of species I is occupied by an impurity, then the species I can transform to species K with rate w I K and the reverse transformation from species K to I occurs with rate w K I .Clearly, this flip dynamics is activated by the presence of the impurities (+).The total number of impurities N + along with the total number of vacancies N 0 are conserved quantities, which can be readily seen from Eq. (1).The complete set of input parameters for the µ-TASEP-IAF is (p I , ε, w I K , ρ + , ρ 0 ), where ρ + = N + /L and ρ 0 = N 0 /L are the conserved densities for the impurities and the vacancies respectively.To illustrate the dynamics, we present a schematic figure of the allowed dynamical processes for the µ = 2 case in Fig. 1.
From the microscopic dynamics in Eq. ( 1), it is clear that starting from a specific initial configuration, the different species and the impurities cannot overtake each other.Indeed the flip dynamics changes the number of accessible configurations by transforming one species to another, but does not allow the dynamics to be ergodic.To discuss the non-ergodicity with an example, let us consider an initial configuration (for 2-TASEP-IAF) of the form {. . .0 + 102011 + 2 . . .}.In the rest of this section, we will denote the particle under consideration by italics e.g. 1 for the chosen particle 1, 2 for the chosen particle 2 etc.If we consider the 1 in {. . .0 + 102011 + 2 . . .}, it can transform into 2 by the + at its right neighbor, thereby changing the configuration to {. . .0 + 102012 + 2 . . .}.However, another 1 in the initial configuration {. . .0 + 102011 + 2 . . .} can never transform to 2 at any stage of the evolution because it can never come in contact with any + (due to the non-overtaking nature of the dynamics), so that the configuration {. . .0 + 102021 + 2 . . .} is never accessible.

Steady state: matrix product ansatz
Any configuration of the µ-TASEP-IAF can be represented by {s i } ≡ {s 1 , s 2 , . . ., s L }, where s i denotes the occupation at site i.Clearly, s i can be one of the species K = 1, 2, . . ., µ or it can be an impurity (+) or it can be a vacancy (0).We find that the steady state of the present model can be written in the following matrix product form In Eq. ( 2), any configuration {s i } is represented by a string of matrices {X i } where the matrices D K , A and E corresponds to a particle of species K, impurity and vacancy respectively.The time evolution of any configuration of the µ-TASEP-IAF is dictated by the Master equation which in steady state becomes M |P〉 = 0.Here |P〉 is a column vector containing all possible configurations and M is the rate matrix made up of the transition rates between configurations.Since the dynamics in Eq. ( 1) is a two-site microscopic dynamics, the transition rate matrix, under the periodic boundary condition, can be expressed as where M i,i+1 is a (µ+2) 2 ×(µ+2) 2 dimensional matrix acting on the pair of sites (i, i +1) and I is (µ + 2) × (µ + 2) dimensional identity matrix placed at every site except the pair (i, i + 1).Then the steady state M |P〉 = 0 of the µ-TASEP-IAF can be achieved through the following two-site flux (probability current) cancellation condition where and where (.) T denotes the transpose of the row vector (.) and Ẽ, Ã, DK are auxiliary matrices that are introduced to satisfy the steady state equation and these have to be found out consistently along with the matrix representations for E, D, D K (K = 1, 2, . . ., µ).We find that suitable choices for the auxiliary matrices for the µ-TASEP-IAF are Correspondingly, the matrices E, A and D K have to obey the matrix algebra consisting of the equations given below The last relation in Eq. ( 9) is reminiscent of the Kirchhoff's current law for each species K, in the sense that the total flip current from all other species to species K is equal to the total flip current from K to all other species.Note that the matrix algebra in Eq. ( 9) allows scalar solutions when the hopping rates for every species and the impurity become equal i.e. p K = ε for all K. Naturally for this special set of rates, since the matrices reduce to scalars, no spatial correlations exist between the constituents of the system.For any other choice of rates, we expect matrix solutions to the Eq. ( 9).Below we discuss the cases µ = 2, µ = 3 extensively with explicit matrix representations and then generalize them to get the matrix representations for general µ.µ = 2 : For the 2-TASEP-IAF (K = 1, 2), the matrix algebra [Eq.( 9)] simplifies to Clearly, the matrix relation for the flip process becomes trivial for the two-species case implying the absence of net flip current between the two species.More precisely, the flip process satisfies detailed balance condition for µ = 2.However, there are non-zero drift currents in the system.We find the following 3 × 3 matrix representations that satisfy the matrix algebra in Eq. ( 10), The matrix representation of the impurity i.e.A, in the projector form, resembles that of the defect of second class particles in case of TASEP with first and second class particles [25,71] except the fact that the matrices are infinite dimensional in Refs.[25,71].µ = 3 : The matrix algebra in Eq. ( 9) for the 3-TASEP-IAF process reads as In comparison to the last relation in Eq. ( 10), clearly the flip processes for the three-species case given by the last three relations in Eq. ( 12) do not require the detailed balance as a necessary condition.Rather, the general condition (without putting any constraint on the set of flip rates) that satisfies the flip processes in Eq. ( 12) is We obtain the following 4 × 4 representations of the matrices that satisfy Eq. ( 12) along with Eq. ( 13), where We note that the condition Eq. ( 13), with the explicit matrix representations from Eq. ( 14), becomes The parameter α in Eq. ( 16) quantifies the deviations of the flip processes from the detailed balance condition between any pair of species.As we would see later, the net flip current between any two species is proportional to α.In fact α = 0, which puts some constraints on the flip-rates, correspond to a straightforward generalization of the two-species process to three-species process with similar flip process matrix relations w K I D K A = w I K D I A .
general µ : For the general case of µ-TASEP-IAF (K = 1, 2, . . ., µ), the structures of the matrices D K , E, A would be similar to that of Eq. (11) and Eq. ( 14).More precisely, the matrices are (µ+1)×(µ+1) dimensional and the corresponding explicit representations of the matrices are given by In Eq. ( 17), the vector 〈I| = (0, . . ., 0, 1, 0, . . .0) where 1 is placed at the I-th element with all other elements being zero and |I〉 is the transpose of 〈I|.Notably, the values of the coefficients d K associated with matrices D K in Eq. ( 17), can be calculated by solving the set of µ homogeneous linear equations of the form The summations over the index I in Eq. ( 18) generally includes all I = 1, . . ., µ except K.This corresponds to the general dynamics in Eq. ( 1) where any two species can transform into one another in the presence of the impurity (at right neighbor).However, one might be interested in special cases where the flip processes are restricted between certain pairs of species only.For example, one particular situation can be where any species K can only transform to species numbers (K + 1) and (K − 1) .In that case, the parameter α [Eq.(16)] which dictates the flipcurrent between any two species, will be given by α this reduces to the general solution of µ = 3 as given in Eq. ( 16).However, unlike the cases of µ = 2 or µ = 3 or the special instance of the multi-species case stated above which allow α to be independent of the species I, in general α would depend on the particulars of the pairs (I, K) for µ > 3. Mathematically, Eq. ( 16) would be generalized as follows where α K I have different values for different pairs (I, K).As an example, one can consider µ = 4 with the allowed set of flip dynamics described in Fig. 2. In Fig. 2, among twelve total possible flip rates, only six are present.Specifically, in this example, the impurity cannot activate flips between species (2, 3), implying α 23 = 0 whereas for other pairs α K I ̸ = 0.

Partition function for special initial configuration
The non-ergodic nature of the microscopic dynamics in Eq. ( 1) ensures that we cannot express the partition function of the µ-TASEP-IAF in the usual form of Tr[T L ], even under periodic boundary conditions.Here the "transfer matrix" T refers to with z 0 and z + being the fugacities corresponding to the vacancies and impurities respectively, in the grand canonical ensemble.This is because T L generates all configurations in the configuration space irrespective of the initial ordering of the species, but the dynamics in Eq. ( 1) allows only those configurations which preserve certain orderings from the initial configuration.To illustrate this with an example, let us consider a string of (n + 1) number of species 1 particles followed by an impurity in the initial configuration i.e. {..11..11 + ..} ≡ ..1 n+1 + .. .At any step in the µ-TASEP-IAF, the ordering of the first n number of species 1 particles in this string cannot be broken, i.e. no species 2 particle or impurity can appear inside this string, except for the last 1 (which may transform to 2).On the other hand, T L generates configurations which do not preserve such orderings.Naturally, to calculate partition function for µ-TASEP-IAF analytically, it becomes essential to choose suitable initial configurations for which we can correctly identify the accessible set of configurations in the steady state.In the rest of this section, we discuss such special initial configurations, corresponding steady states and partition functions.

µ = 2
One special initial configuration C(t = 0) ≡ C(0) (represented by matrices) for the 2-TASEP-IAF is where X ...X represents an uninterrupted sequence of the matrix X .We consider the densities of the uninterrupted sequences of D 2 -s and D 1 -s to be equal, which is ρ = N /L.Further, we have taken the two sequences of D 1 A and D 2 A to be of equal density ρ + so that the density of impurity (A) in each of these sequences is ρ + /2.This ensures that the total density of impurities is ρ + .The initial configuration in Eq. ( 21) satisfies the relation ρ = 1 2 (1−ρ 0 −2ρ + ), where ρ 0 and ρ + are densities of the vacancies and impurities respectively.In the steady state we have ρ 0 + ρ + + ρ 1 + ρ 2 = 1, with ρ 1 and ρ 2 being the average densities of species 1 and species 2 particles in the steady state.We emphasize that ρ 0 and ρ + are input parameters while ρ 1 and ρ 2 are derived quantities.Starting from Eq. ( 21), any accessible configuration C ss in the steady state is of the following generic form subjected to the constraint The parameter τ in Eq. ( 22) can take value either 1 or 0. The partition function, obeying the above mentioned constraint, is given by It would be useful to get rid of the δ(.) constraint by associating a fugacity z 0 to the vacancy (represented by E) and considering the system in a grand canonical ensemble.The matrix strings (D 1 + D 2 )E m AE n , D 2 E r and D 1 E s can be evaluated by incorporating the matrix algebra in Eq. ( 10) along with the explicit representations from Eq. ( 11).It should be mentioned that the projector form of A [Eq. (11)] leads to factorization of the matrix strings, e.g.
which helps significantly in carrying out the analytical calculations.Consequently, the partition function in the grand canonical ensemble under the periodic boundary condition, finally becomes For the special initial configuration in Eq. ( 21), we have derived the partition function in Eq. (25).The fugacity z 0 can be obtained as a function of the vacancy density and other input parameters by inverting the density-fugacity relation In general, the solution of z 0 obtained from Eq. ( 26), using Mathematica, appears to be complicated and lengthy.However, in Appendix D, we would discuss two special cases (with specific choices of the input parameter) that provide closed form solutions for the fugacity.The other conserved quantity, the impurity density is already fixed at ρ + = N + /L.The expression Eq. ( 25) would be used for evaluating the observables of interest in the next Sec.4.

µ = 3
Similar to the initial configuration C(t = 0) for µ = 2, a suitable initial configuration for the three species case that enables us to perform analytical calculation of partition function and observables, is The density of the uninterrupted sequence of each species 1, 2, 3 in Eq. ( 27) is taken to be equal to ρ = N /L.Moreover, we have chosen the density of each of the sequences D 1 A, D 2 A and D 3 A to be 2ρ + /3 where the density of impurities in each of these sequences is ρ + /3.Consequently, the choice of initial configuration in Eq. ( 27) ensures that the total density of impurities remain ρ + .In the steady state, ρ 0 +ρ + + , where ρ I is the average density of the non-conserved species I. Starting from Eq. ( 27), the form of any accessible configuration C ss in steady state would be with the conservation of the total number of vacancies where δ τ,K is the Kronecker delta symbol with K = 1, 2, 3.As done in case of µ = 2, here also we associate a fugacity z 0 with the vacancy.Using the matrix algebra from Eq. ( 12) alongside the matrix representations in Eqs.( 14) and ( 15), the partition function in grand canonical ensemble becomes where the explicit expressions for d I -s have been presented earlier in Eq. ( 15).

General µ
For the µ-TASEP-IAF (K = 1, . . ., µ), the generalization of the initial configurations in Eq. ( 21) (µ = 2) and Eq. ( 27) (µ = 3) would be The above initial configuration is chosen in a way that the density of each sequence D I A (I = 1, . . ., µ) is 2ρ + /µ in which the density of impurities is equal to ρ + /µ, so that the total impurity density adds up to ρ + .We have where d I is the solution of Eq. ( 18) and z 0 is the fugacity associated with the vacancy in the grand canonical ensemble.
In this section, we have derived the partition functions of the µ-TASEP-IAF with µ = 2, µ = 3 and general µ in Eq. ( 25), Eq. ( 29) and Eq. ( 31), respectively, for specific initial configuration Eq. ( 21), Eq. ( 27) and Eq. ( 30), respectively, under periodic boundary conditions.These results would be useful to calculate the average values of observables in the next section for the same initial configurations discussed here.

Observables: comparisons of analytical results with Monte Carlo simulations
In this section, we analytically calculate the following observables in the steady state, (i) average density ρ I of the non-conserved species I, (ii) average drift currents J I0 and J +0 , for species I and impurities respectively, (iii) average flip current J I↔K between species pair (I, K) and (iv) two-point correlations C 0I between vacancies (0) and species I. Mostly we will restrict the calculations to the number of species µ = 2, except using µ = 3 for the case of average flip current (since there is no net flip current between pair of species for µ = 2).In particular, starting from the special initial configuration Eq. ( 21) for µ = 2 (or, Eq. ( 27) for µ = 3), we will show agreements between the analytical calculations and the Monte Carlo simulation results.

Species densities
First we consider the average densities (ρ I ) of the non-conserved species I = 1, 2. The formal expression for ρ I in the steady state under the periodic boundary condition, can be written as ..
To elaborate Eq. ( 32), the main point is to note the expression inside the trace (Tr[.]) that denotes configurations with at least one D 1 .This can be understood more clearly by comparing it with the expression for any possible configuration in Eq. ( 22).From Eq. ( 22), to arrive at the matrix string inside the trace in Eq. ( 32), one has to put τ = 1 for one k value to ensure the presence of at least one species 1 particle (D 1 ) in the configuration.Since this D 1 could have been placed for any k = 1, 2, . . ., N + , we have a combinatorial pre-factor ρ + = N + /L in Eq. (32).The summations over all the variables {m, n, r, s} from zero to infinity are performed as the system is considered in the grand canonical ensemble.The first factor (1 in Eq. ( 32) is to take care of the N number of D 1 -s present in the initial configuration Eq. ( 21) that cannot flip.Using the matrix algebra and matrix representations from Eqs. ( 10) and ( 11) respectively, we finally arrive at the following expressions for the average densities of the non-conserved species, The fugacity z 0 is obtained in terms of the input parameters by solving the density-fugacity relation Replacing the solution of z 0 (which is too lengthy to provide here) in Eq. ( 33), we finally get the average densities as functions of input parameters (p 1,2 , w 12,21 , ε, ρ 0 , ρ + ) only.
We compare the analytical results with those of Monte Carlo simulations, starting from the same initial configuration Eq. ( 21).In the simulation, we vary the vacancy density in the initial configuration by changing the lengths of the uninterrupted strings of D 2 and D 1 in Eq. ( 21) i.e. simply by tuning ρ which is related to vacancy density ρ 0 as ρ 0 + 2 ρ = 1 − 2ρ + .In Figs.3(a) and (b), we observe that the analytical and simulation results are in agreement with each other, where ρ I is plotted against ρ 0 and w 12 (flip rate of species 1 to species 2), respectively.The species densities decrease linearly with increasing ρ 0 [Fig.3 The parameters used are L = 10 3 , p 1 = 0.3, p 2 = 1.0, ε = 0.1, w 12 = 0.4, w 21 = 1.0, ρ + = 0.2.The ensemble average is done over 10 5 samples.shifts to Consequently in Fig. 3(b), for a particular set of chosen parameters, we observe that for w 21 < w 12 < w ⋆ 12 , one still has ρ 2 < ρ 1 .In other words, when w 12 ∈ (w 21 , w ⋆ 12 ), although the species 1 particles more often transform to species 2 particles, still the average density of species 2 particles is less than that of species 1 particles.For any value of µ, the general expression for the average density ρ I for the non-conserved species I (I = 1, 2, . . ., µ) is where d I is the solution of Eq. ( 18) (e.g. the solution for µ = 3 is explicitly given in Eq. ( 15)).

Drift current
Next we consider the average drift currents J I0 and J +0 for the non-conserved species I (I =1, 2) and the impurity respectively.We focus on I = 1 to explain the procedure for calculating the current J 10 , because the parallel procedure applies for any other species.The average drift current J 10 is equal to p 1 〈10〉, where 〈10〉 is the ensemble average of the pair 10.In terms of matrices the expression 〈10〉 simply translates to 〈D 1 E〉.The current J 10 can be calculated in two parts, 10 + J (2) where J (1) (1) is the contribution from the drift of species 1 particles that can flip and J (2) (2) is the corresponding contribution from species 1 particles that cannot The nonlinear monotonic decrease and increase in J 10 and J 20 respectively are similar in nature to the effect of w 12 on the species densities (Fig. 3).The impurity drift current increases slowly with increasing w 12 .The parameters used are The ensemble average is done over 10 5 samples.
flip (as they cannot have any impurity as right neighbor) according to the initial configuration in Eq. (21).Correspondingly, the term D 1 E in the averages 〈D 1 E〉 (1) and 〈D 1 E〉 (2) would come from the product sequences (τD 1 + (1 − τ)D 2 )E m A and D 1 E s respectively Eq. ( 22).The expression for J (1) 10 is given by ..
The construction of Eq. ( 37) follows similar arguments as of Eq. ( 32), except now we have to place D 1 E instead of D 1 .This also reflects in the summations, note that the lower limit of the index m 1 has been changed to 1 instead of 0 to ensure the presence of one D 1 E. Similarly, the formal expression for J (2) Obviously in Eq. ( 38), the lower limit of the index s 1 is shifted to 1 from 0. Using the matrix algebra and matrix representations from Eqs. ( 10) and ( 11), we get from Eqs. ( 37) and (38): , J (2) Substituting Eq. ( 39) into Eq.( 36), we obtain J 10 .Following the same procedures, one can calculate J 20 and J +0 .We finally arrive at the analytical expressions for the drift currents under the periodic boundary condition, given below: In Figs. 4 and 5 we present the variation of the drift currents as functions of vacancy density ρ 0 and flip rate w 12 , respectively.For both cases, the analytical results match with the Monte Carlo simulation results.The drift currents for the species 1 and 2 exhibit non-monotonic behaviors with increasing vacancy density whereas the impurity drift current increases monotonically (Fig. 4).At lower vacancy densities, as we increase ρ 0 , the chances for hopping increase, thereby increasing the drift current in Fig. 4.However, after a particular value of ρ 0 , if we increase it further, the densities of the non-conserved species fall considerably so that drift current ultimately decreases, although there are many vacancies in the system.Since the impurities do not flip, the density of impurities is fixed and the corresponding impurity drift current can only increase with increasing ρ 0 (Fig. 4).The maximum of the drift current for different species generally occur at distinct values of ρ 0 .With variation of the flip rate w 12 (Fig. 5), the drift currents of the non-conserved species show similar non-linear behaviors like their corresponding densities [Fig.3(b)].Notably, although the flip dynamics does not affect the drift of the impurities explicitly, still we observe that the impurity drift current increases with increasing w 12 , albeit weakly.For any positive integer value of µ, the general expression for the drift current of any non-conserved species I (I = 1, . . ., µ) is obtained to be where the d I is the solution of Eq. ( 18).

Two-point correlations
Besides the currents, we would like to calculate some other two point functions which have interesting features.We have basically calculated the two-point functions 〈10〉 and 〈20〉 in the process of determining the drift currents J 10 and J 20 respectively.Now we focus on the nearest neighbors two-point correlations involving 〈01〉 and 〈02〉.We find the exact expressions for the corresponding two-point correlations under the periodic boundary condition as We observe that this correlation initially starts increasing with increasing ρ 0 , reaches to a local maximum, then decreases.However, as ρ 0 is increased further, C 01 reaches to a local minimum and then starts increasing again.The parameters used are L = 10 3 , p 1 = 0.3, p 2 = 1.0, ε = 0.1, w 12 = 0.8, w 21 = 1.0, ρ + = 0.2.The ensemble average is done over 10 5 samples.
The correlation C 01 is plotted against the vacancy density ρ 0 in Fig. 6.As ρ 0 is increased starting from zero, the correlation also increases and reaches a local maximum, followed by a decrease and reaching a local minimum.After this point, if the vacancy density is increased further, C 01 increases again.So, instead of a single maximum or single minimum, the correlation C 01 interestingly exhibits both local maximum and local minimum with the variation of ρ 0 .In other words, C 01 increases with increasing vacancy density for both sufficiently high and sufficiently low values of ρ 0 , with intermediate non-monotonic character.In Fig. 7, it is interesting to see that the non-monotonic behavior of C 02 is such that it goes from negative correlation values to positive correlation values.Naturally, there exists some intermediate value of ρ 0 for which the special arrangements of the accessible steady state configurations makes the average correlation C 02 to be zero.

Flip current
All the observables we have discussed up to now (average species densities, drift currents, correlations), correspond to 2-TASEP-IAF.However, the net flip current is zero for µ = 2. Therefore, in order to have a net non-zero flip current between pairs of species, here we consider the case µ = 3 (I = 1, 2, 3).We denote the net flip current between species I and K as J I↔K .
For µ = 3, we find that the net flip current between any two species [(1, 2), (2, 3), (3, 1)] are equal to each other (i.e.independent of the indices I and K) and its exact form is given by The initial configuration that we have used to arrive at Eq. ( 44) is the one in Eq. ( 27).In Eq. ( 44), we have calculated the current in the cyclic ordering i.e. (I = 1, K = 2), (I = 2, K = 3), (I = 3, K = 1).We have presented the behavior of the flip current as functions of the vacancy density ρ 0 and flip rate w 12 in Figs.8(a) and 8(b) respectively.We observe that the analytical calculation are in agreement with the Monte Carlo simulation results.In Fig. 8(a), the flipcurrent decreases monotonically in a nonlinear manner with increasing vacancy density.The reason behind this, is the decrease in species densities with increasing ρ 0 [Fig.3(a)].On the other hand, flip current between any pair of species increases monotonically with increasing flip rate w 12 (Fig. 8(b)).The generalization of the formula Eq. ( 44) for any µ, under the periodic boundary condition, is obtained as where d I , d K are the solutions of Eq. (18).We should mention that, for µ > 3, the flip currents between different pairs of species would be in general distinct from one another.

Non-ergodicity: dependence on initial configuration
Here we establish the non-ergodicity of the µ-TASEP-IAF by showing explicitly the dependence of the average steady state values of observables on the choice of the initial configuration.
For simplicity, we restrict ourselves to the case of µ = 2.We choose two different initial configurations as follows.(i) The initial configuration given in Eq. ( 21) for which we know exactly which constituent (any species or impurity or vacancy) is placed at a given lattice site.
Since the initial arrangement is specified completely, we call this configuration specified initial configuration (SIC).(ii) A random initial configuration where the constituent at each site is selected randomly such that the densities of impurities and vacancies, and the initial densities of species I (I = 1, 2), are exactly the same as that of the SIC described in (i).Since the initial arrangement for this configuration is randomly carried out, we call it random initial configuration (RIC).To clarify, both SIC and RIC are characterized by the same set of rates (p 1 , p 2 , ε, w 12 , w 21 ) and the densities (ρ + , ρ 0 , ρ 1 (0), ρ 2 (0)), where ρ 1 (0), ρ 2 (0) represent the initial (t = 0) densities of the non-conserved species 1 and 2 respectively.Although the analysis of the steady state for the SIC can be performed exactly as discussed already, the same could not be done for the RIC.Therefore, in this section we use Monte Carlo simulations to compare the steady state observable values for SIC and RIC.We compare the steady state observable values for SIC and RIC with the same set of input parameters (p 1 , p 2 , ε, w 12 , w 21 , ρ + , ρ 0 ) and same initial densities of the species (ρ 1 (0), ρ 2 (0)) in Fig. 9.We denote the data points for SIC and RIC with different symbols, circles and rectangles respectively, in Fig. 9.The variations of the non-conserved species densities ρ I are presented in Figs.9(a) and 9(c) as functions of flip rate w 12 and hop rate p 2 , respectively.Both figures exhibit clear quantitative differences between the density values for SIC and RIC.We observe that the deviations between SIC and RIC decreases (increases) with increasing w 12 (p 2 ).The initial configuration dependence of the steady state values of nearest neighbor two-point correlations are shown in Figs.9(b) and (d).In Fig. 9(b), we observe that the correlation between species 1 particles C 11 = 〈11〉 − ρ 2  1 has distinct numerical values for SIC and RIC when the parameter w 12 is tuned.Interestingly, the correlation corresponding to RIC changes from negative to positive whereas the same for SIC remains positive with increasing w 12 .This implies the existence of some intermediate w 12 which corresponds to uncorrelated species 1 particles for RIC, whereas they are correlated for the SIC.Similar kind of interesting behavior is observed for the correlation between species 2 particles C 22 = 〈22〉 − ρ 2  2 when plotted against p 2 in Fig. 9(d).Thus we have illustrated the dependence of steady state values of species densities and correlations on the choice of the initial configuration.The same can also be investigated in other two point and higher point functions.We end this section with a general comment regarding the non-ergodicity in the present model.If we consider a sequence of the form {+s i s i+1 . . .s n s n+1 +} in an initial configuration, where s j = 0, 1, . . .µ but s j ̸ = + for i ⩽ j ⩽ (n+1), then the ordering of different species 1, . . ., µ (not vacancies) for i ⩽ j ⩽ n remains intact for the allowed subspace of configurations in the steady state.Naturally, number of such orderings increase with system size.Recent study [72] in context of classical reversible cellular automaton shows the number of local conservation laws increase exponentially with system size, leading to block diagonal form of the propagator with exponential scaling of the number of blocks with system size.In fact there are quantum systems like certain Lindbladian for quantum ASEP [73], dipole-conserving Hamiltonian [74] etc. for which the space of operators or states fragment into invariant subspaces whose number again scale exponentially with system size.It would be interesting to investigate in future how does the number of conserved orderings scale with system size in our non-ergodic model, the detailed block diagonal structure of the transition rate matrix [M in Eq. ( 3)] and the role of corresponding underlying symmetries.An explicit illustration of the block-diagonal structure of the rate matrix in the µ-ASEP-IAF, for small system sizes, is presented in Appendix E.

Partially asymmetric generalization: µ-ASEP-IAF
In this section, we consider the µ-ASEP-IAF under periodic boundary conditions, a generalization of the µ-TASEP-IAF in Eq. ( 1), by including partially asymmetric motions of different species of particles.A particle of species I can hop towards right with rate p I and it can hop towards left with rate q I (I = 1, . . ., µ), if the target site is empty.Notably, the impurities are not allowed to hop to left.This naturally adds another way to distinguish the conserved impurities from all non-conserved species.The microscopic dynamics is given by, drift (species) : The µ-ASEP-IAF remains non-ergodic in nature.Since we could obtain the steady state of the µ-TASEP-IAF using matrix product ansatz [Eq.( 2)], we assume the same can be done for the partially asymmetric motion also.

Matrix algebra, auxiliaries and matrix representations
The matrix algebra for the dynamics in Eq. ( 45) under the periodic boundary condition, is In comparison to the matrix algebra [Eq.( 9)] for the µ-TASEP-IAF, the only changes occurring in Eq. ( 46) correspond to the drifts of the non-conserving species.At this point, we should mention that the matrix equation p K D K E − q K E D K = D K has been studied in Ref. [55], in context of a conserved disordered ASEP model.Due to the presence of the impurities and the flip processes activated by them, the matrix algebra for µ-ASEP-IAF in Eq. ( 46), can be considered as a generalization of the matrix algebra in Ref. [55].To arrive at the matrix algebra in Eq. ( 46) from the dynamics (45), ansatz (2) and flux cancellation condition (5), the choice of the auxiliary matrices are the same as the totally asymmetric case, i.e.
However, unlike the totally asymmetric case, we find the matrix representations for the µ-ASEP-IAF to be infinite dimensional.Notably, this does not necessarily eliminate the possibility of getting alternate finite dimensional representations of the matrices.Below we present the matrix representations for µ = 3 case explicitly (as we will stick to µ = 3 for the discussion of observable in this section) and mention the changes required to construct the matrices for any µ > 0.
In the absence of the impurities (A) and the flip processes, the term d 1,1 µ for every µ becomes unity, and corresponding the matrix representations for D I and E in Eq. ( 48) (and their generalizations for general µ) are the same as that of the conserved disordered ASEP [27,55].In Eq. ( 48 is the solution of Eq. ( 18).

Partition function for special initial configuration
In the totally asymmetric case, for general µ, we have considered the specific initial configuration Eq. ( 30) which leads us to acquire analytical expressions for observables of interest.We choose a particular case of Eq. ( 30), namely the ρ = 0 case, as the special initial configuration for analytical calculation in µ-ASEP-IAF.More precisely, the choice of our special initial configuration for µ-ASEP-IAF is, The initial configuration Eq. ( 49) is chosen in a way that fixes the total impurity density to ρ + .In comparison to Eq. ( 30), the initial configuration in Eq. ( 49) is simpler and does not contain consecutive D I -s.We would see shortly that this specific choice is sufficient to show negative differential mobility in µ-ASEP-IAF.We find the partition function in the steady state under the periodic boundary condition corresponding to the initial configuration Eq. ( 49), to be We have used short hand notations which are essentially the solutions of Eq. ( 18).When q I = 0 for all species, it is straightforward to check that the partition function in Eq. ( 50) reduces to the partition function Eq. ( 31) of the totally asymmetric case, under the condition ρ = 0.

Species densities, drift current and flip current
Just like we did in the µ-TASEP-IAF, we can analytically calculate several observables of interest in the partially asymmetric case also, using the matrix algebra (46) and matrix representations (48) following the same procedures as before.Starting from the initial configuration stated in Eq. ( 49), the average density ρ I of any non-conserved species I (I = 1, . . ., µ) for the µ-ASEP-IAF is obtained as where ρ 0 and ρ + are the conserved densities for the vacancies and the impurities, respectively.If we put q I = 0 for all I in Eq. ( 51), the expression of ρ I for the totally asymmetric case [Eq.(35)] is correctly recovered.The drift currents J I0 and the flip currents J I↔K for the µ-ASEP-IAF (I, K = 1, . . ., µ) are given by

Negative differential mobility
In what follows, we will show that the species in the µ-ASEP-IAF under the periodic boundary condition, exhibit negative differential mobility [75,76].More precisely, we would see that both the drift currents and the flip current can decrease with increasing bias (which we define later), giving rise to the phenomena of negative differential mobility (NDM).NDM has been observed for driven tracer particles in the presence of static obstacles [78,79] or in crowded medium [77] and for many particle systems in presence of kinetic constraints [80] or obstacles [81].There have been many studies to understand the mechanism of NDM in driven systems and it appears that some kind of trapping that leads to decrease in dynamical activity, acts as a main cause of NDM [79,82,83].In connection to asymmetric simple exclusion process, a two dimensional variant of ASEP where the kinetic constraint is implemented by restricting the motion of the particles depending on the number of its occupied neighbors, has been shown to exhibit NDM at high density and high bias values [80].In one dimension, a single driven tracer hopping asymmetrically in the environment of bath particles executing symmetric exclusion process, exhibits negative differential mobility as well as absolute negative mobility (current flowing in a direction opposite to the bias direction), where the kinetic constraint is imposed by an additional exchange dynamics of the tracer with a distant bath particle depending on the vacant nearest neighbors [84].Another way to incorporate the effect of the kinetic constraint leading to NDM, is to consider the escape rate from a configuration as a decreasing function of the bias, shown elaborately for a biased random walker in Ref. [79].
Recently in Ref. [85], the authors have proposed that slowing down of non-driven degrees of freedom (modes) through the biasing of the driven mode, can give rise to negative differential mobility for both the driven and non-driven degrees of freedom in an interacting many particle system.Here, we apply this mechanism to show that indeed the µ-ASEP-IAF can exhibit NDM for particular choices of the rates in the microscopic dynamics.
To illustrate NDM in µ-ASEP-IAF, we will focus on the µ = 3 case.We have three species of particles (I = 1, 2, 3), impurities (+) and vacancies in the system following the microscopic dynamics Eq. ( 45).We choose the drift rates p I and q I of the species I = 1, 2, 3 to be The choices of the hopping rates in Eq. ( 53) are inspired by similar choices in Ref. [85] in context of NDM for different models.The special choices of the hopping rates in Eq. ( 53) allow us to identify the parameter b as the hopping bias in the system.This is because ln(p 1 /q 1 ) = b and the unbiased case p 1 = q 1 = 1 corresponds to b = 0.Then, the species 1 particle is a driven mode for any b > 0.Even when b > 0, Eq. ( 53) clearly states that species 2 and species 3 particles are non-driven modes in the system because the right and left hopping rates are equal for both of them.However, there is a key difference between the hopping rates of species 2 and species 3 particles.The hopping rates of species 2 depend explicitly on the bias b of the driven mode (species 1).More precisely, the hopping rates p 2 and q 2 decrease with increasing bias b.This corresponds to the slowing down of non-driven mode and turns out to be the key for negative differential mobility.The hopping rates p 3 and q 3 of the other non-driven mode species 3, does not depend on b.Here we should mention the presence of another driven mode in the system, which is the impurity.Since, the impurity motion is only unidirectional, it is a driven mode by construction.Consequently, even at b = 0 the system is in a non-equilibrium steady state for ε > 0 where the impurity acts as the lone driven mode.In the present context, we choose ε to be a constant independent of the bias b.We focus on the behavior of the currents with the variation of b.To summarize, the driven modes in the 3-ASEP-IAF are (i) species 1 (driven by b) and (ii) impurity (driven due to unidirectional motion with constant rate ε, independent of b), whereas the non-driven modes are (i) species 2 (hopping rates are decreasing function of b) and (ii) species 3 (hopping rates independent of b).All the flip rates w I K (I, K = 1, 2, 3) are kept constants independent of the drift bias b.With this set up, we now investigate the variation of the drift currents and flip current as functions of the bias b, both from analytical formulae Eq. ( 52) and Monte Carlo simulations.In Fig. 10(a), we present the behaviors of the drift currents of the driven modes with variation of the bias b, under the periodic boundary condition.Interestingly, although the bias b is directly applied to species 1 to increase its current, the drift current for species 1 decreases monotonically with increasing bias giving rise to the phenomena of negative differential mobility.The current of the other driven mode, the impurities, initially increase with increasing bias, reaches to a maximum, but then decreases as the bias is further increased, thereby leading to NDM.The drift currents of both the non-driven modes exhibit non-monotonic behaviors with increasing bias as shown in Fig. 10(b).Both of them decrease with increasing bias for sufficiently large values of b, showing negative differential mobility.Notably, the flip dynamics is not directly affected by the drift bias since all the flip rates are kept to be constants independent of b.Therefore, it is intriguing to observe that the net flip current still decreases with increasing bias (for large b) and therefore exhibits NDM, as presented in Fig. 10(c).The mechanism behind the negative differential mobility in drift current is related to the decreasing dynamical activity (number of hops per unit time) of the species 2 particles (one of the non-driven modes) with increasing forward bias b for the species 1 particles.Since with increasing b, the hop rate of species 2 [Eq.(53)] decreases, its waiting time at the residing site increases i.e. it becomes more and more prone to stay at the residing site rather to leave the site as b increases.That is why, although the increasing bias tries to push particles forward, their ways are blocked by the slowed down species 2 particles.The exclusion interaction and the non-overtaking dynamics facilitates the NDM even better by not allowing other species or impurities to overtake the slowed down species 2 particles.The reason behind the negative differential mobility occurring in the flip current requires further investigation.
We end this section with mentioning the possibility of further nontrivial transport properties in the steady state of the µ-ASEP-IAF when one considers the counter flow scenario.The counter flow in the system arises when the net bias of some species of particles are opposite to that of the others.For example, in the 2-ASEP-IAF, species 1 can have net bias to right i.e. p 1 > q 1 whereas the species 2 can have net bias in the opposite direction i.e. q 2 > p 2 .Counter flow can give rise to interesting physical features e.g.phase transitions [48,52,86].This urges for detailed investigation of the counter flow situation in µ-ASEP-IAF in future works.

Summary and future directions
In this article, we have obtained an exact steady state probability distribution of the µ-ASEP-IAF on a one dimensional lattice under periodic boundary conditions, using the matrix product ansatz.The µ-ASEP-IAF consists of (i) drift of the species (I = 1, 2, . . ., µ) and impurities, and (ii) flip between different species initiated by the impurities.In steady state, we provide the explicit finite dimensional [(µ+1)×(µ+1)] matrix representations for any µ > 0 for the totally asymmetric case i.e. µ-TASEP-IAF.For the partially asymmetric scenario i.e. µ-ASEP-IAF, we obtain the corresponding matrices with infinite dimensional representations.Importantly, due to the non-ergodicity of the µ-ASEP-IAF dynamics, the partition function and observables in the steady state depend on the specific choice of the initial configuration.However, for a special class of initial configurations, we could indeed analytically calculate the partition function for both the totally asymmetric and partially asymmetric cases with any µ > 0, under periodic boundary conditions.We present exact analytical expressions for steady state observables like the average densities of the non-conserved species, drift current, flip current and some other two-point correlations.We show that our analytical calculations are in agreement with the Monte Carlo simulations for the analytically tractable specific initial configuration.In this connection, the non-ergodicity of the model has been established extensively (Monte Carlo simulations) by showing the deviations of the steady state observable values for a random initial configuration from that of the specific initial configuration mentioned above.Along with the important exactly solvable analytical structure, the µ-ASEP-IAF also has interesting physical features.Notably, with the variation of vacancy density, several two-point correlations exhibit interesting behaviors e.g.transiting between negative and positive correlations, showing both local maximum and local minimum etc.The effect of the drift on the flip processes are evident from the functional dependence of the species densities on the flip rates.Interestingly, both the drift current and flip current in the µ-ASEP-IAF are shown (analytically and numerically) to display negative differential mobility (decreasing current with increasing bias) for certain choices of the drift rates.The mechanism behind the negative differential mobility relies on slowing down a non-driven mode in the system through the biasing of a driven mode, which eventually leads to decreased dynamical activity of all the modes in the steady state.
Apart from its own intriguing mathematical and physical characteristics, the µ-ASEP-IAF studied here is relevant in two other important contexts.The µ-ASEP-IAF has interesting connections to (i) multi lane asymmetric simple exclusion proces (m-ASEP) which serves as a simple yet remarkable model for multi lane traffic flow, and (ii) enzymatic chemical reactions.For the totally asymmetric model µ-TASEP-IAF, these connections are discussed in details in Appendix A and Appendix C respectively.Importantly, the exact solution of µ-ASEP-IAF suggests possible exact solutions in corresponding multi lane ASEP and traffic models with correlations between particles in different lanes and non-zero net current between lanes.The detailed and rigorous analysis to develop these connections and incorporate them for studying multi-lane traffic flow, constitutes one of the main future directions.We also propose a variation of the µ-ASEP-IAF that exhibits better prospects for being a model for multi lane traffic flow (see Appendix B).In future we plan to investigate the exact steady state and observables of this varied µ-ASEP-IAF model using matrix product ansatz.Just like the multi lane traffic flow, the connections between µ-ASEP-IAF and enzymatic chemical reactions both in steady state as well as dynamics, should be analyzed in more details by considering observables relevant for the chemical reactions.In the present article, we have considered impurities with fixed finite density and constant hopping rate.It would be useful to explore the effects of the variations of impurity density and impurity hopping rate on the observables in the µ-ASEP-IAF.Another important future direction would be to analyze the exact steady state of the µ-ASEP-IAF with open boundary conditions which is more pertinent in context of transport processes, and might also lead to rich phase transitions.It would be interesting to look into the effect of counter flow (i.e.some species having net bias in opposite direction relative to the other species) on the transport properties and possibility of phase separations in the µ-ASEP-IAF.We would also like to investigate the dynamics of the µ-ASEP-IAF in detail, in particular if the product form of the steady state also prevails in the dynamics (using dynamical matrix product ansatz), the dynamical activity in terms of large deviations and possibility of dynamical phase transitions in related models.We can generalize the approaches described above to establish connections between the multi-lane TASEP and the µ-TASEP-IAF for any µ ⩾ 3. It is noteworthy that, for the multispecies case (µ ⩾ 3), we have shown the existence of non-zero net flip current in the µ-TASEP-IAF.It implies the existence of net non-zero lane change current between neighboring lanes in the multi lane TASEP.Also, the correlations between different species of particles and vacancies in the µ-TASEP-IAF suggest non-zero correlations between particles in the different lanes in the m-TASEP.The mapping can also be extended for the partially asymmetric motion of particles.We must mention that the connections between the muti-lane TASEP and µ-TASEP-IAF described here, are approximate.To establish more accurate relations between the two processes, one has to perform rigorous calculations for observables in the multi-lane TASEP and compare the corresponding results with that of the µ-TASEP-IAF.

B A variation of µ-TASEP-IAF, connection to multi-lane traffic flow
The multi lane TASEP has been widely regarded as a simplistic yet important model for multi lane traffic flow [2].Due to the connections between the µ-TASEP-IAF and the multi-lane TASEP discussed in Appendix A, it is natural to ask about the applicability of µ-TASEP-IAF [Eq.(1)] as a suitable model for multi lane traffic flow.Before addressing this question, we should mention that the lane change dynamics in realistic traffic flow must facilitate the traffic as a whole.More precisely, the change of lanes should increase the total flow or total current along the lanes.To investigate this for the µ-TASEP-IAF with µ = 2, we plot the total drift current of species 1 and species 2, J total = J 10 + J 20 as a function of the flip rate w 12 in Fig. 13.Of course, for the two lane TASEP, this amounts to investigating the variation of the total drift current of lane 1 and lane 2 by changing the lane change rate w 12 .
In Fig. 13 we observe that the total current, although increases with the flip rate, the amplitude of the increment is quite small keeping in mind the wide range of variation in the tuning parameter w 12 ∈ (0, 1).The reason behind this, as revealed by a careful observation, is the approximate mapping between the µ-TASEP-IAF and multi lane TASEP described in Fig. 12.In the µ-TASEP-IAF dynamics, when a species encounters an impurity, it flips but does not change its position.On the other hand, in the multi lane TASEP, when a particle in any lane comes in contact with an active bridge, it actually changes the lane i.e. not only changes its characteristics (lane 1 particle to lane 2 particle or vice versa) but also changes its position.To incorporate this in our present model, we propose a variation of the µ-TASEP-IAF dynamics in Eq. (1).Specifically, the change is made only in the flip dynamics.Earlier in Eq. ( 1), the flip dynamics has been with I, K = 1, ..., µ and I ̸ = K.Whereas we propose the new flip dynamics to be Note that, in comparison to Eq. (B.1 ), the flip process in Eq. (B.2 ) accompanies the flip of the species with a hop towards right.The other hopping processes in Eq. ( 1) remain the same for this varied µ-TASEP-IAF.The effect of the dynamics can be immediately observed in Fig. 14 where we present the variation of the total drift current as a function of the flip rate (lane change rate) for both the dynamics in Eq. (B.1 ) and Eq.(B.2 ) (using Monte Carlo simulations).Indeed, the Fig. 14 shows that J total increases considerably with increasing lane change rate for Eq.(B.2 ) whereas it grows weakly for Eq.(B.1 ) (see Fig. 13).This observation implies that the proposed variation of the µ-TASEP-IAF acts as a better model for multi lane traffic flow in comparison to the original model.It would be interesting to study this variation of the µ-TASEP-IAF both analytically and numerically and to build connections with the multi lane traffic flow.

C Connection between µ-TASEP-IAF and enzymatic chemical reactions
In this appendix, we briefly discuss some connections between the µ-TASEP-IAF and enzymatic chemical reactions.One of the simplest form of the enzymatic chemical reaction is, where E, S, P denotes enzyme, substrate, product respectively and ES corresponds to the intermediate complex.The parameters k f and k b are the rate constants for forward and backward reaction for the intermediate complex formation, while kf and kb are the rate constants for forward and backward reactions between the intermediate complex and the product (along with enzyme).Clearly, the initially present enzyme in the reaction remains intact after the reaction is completed.Here, to discuss some connections to the µ-ASEP-IAF, we would rather consider a much simplified version of the chemical reaction (C.1 ) as where we ignore the intermediate complex formation.k S (k P ) is the rate constant for S transforming to P (P transforming to S).
Let us consider a spatially extended narrow channel with many units of substrates, enzymes and products all of which drift through the channel at different rates.This system of chemical reagents can be approximately mapped to an equivalent µ-TASEP-IAF.As explained in Fig. 15 for µ = 2, the impurity in 2-ASEP-IAF plays the role of the enzyme, as it can transform one species of particle to another species.One species, e.g.species 1 can be considered as S for chemical reaction Eq.(C.2 ), whereas the species 2 particle acts as P. The flip rates w 12 and w 21 mimic the rate constants k S and k P .With this set up, our study of the 2-TASEP-IAF reveals the effect of drift on the resultant concentrations of substrates and products in the steady state.Notably, the multi-species (µ > 2) case of the µ-TASEP-IAF can be mapped to generalized version of the chemical reaction Eq.(C.where i = 1, . . ., µ 1 are identified as µ 1 number of substrates and j = 1, . . ., µ 2 are identified as µ 2 number of products in the enzymatic chemical reaction system (with E acting as the enzyme for each chemical reaction) where µ 1 + µ 2 = µ (µ being the total number of species in the µ-TASEP-IAF).To study the time evolution of the enzymatic chemical reaction, one has to study the dynamics of the µ-TASEP-IAF.The connection between the µ-TASEP-IAF and the enzymatic chemical reactions should be studied more thoroughly with proper attention to the observables of interest for the chemical reactions.
D Solution of Eq. ( 26) for fugacity z 0 : special cases We have calculated the partition function (Sec.3) and observables (Sec.4) in the grand canonical ensemble, by associating a fugacity z 0 with the vacancies.Since the observables have to be finally expressed in terms of the input parameters (p 1 , p 2 , ε, w 12 , w 21 , ρ 0 , ρ + ) only, an important step in the calculation is to solve Eq. ( 26) to obtain the fugacity as a function of these input parameters i.e. z 0 (p 1 , p 2 , ε, w 12 , w 21 , ρ 0 , ρ + ).However, in most cases the solutions of z 0 from Eq. ( 26), cannot be obtained explicitly.In this appendix, we would provide two simple cases for particular choices of the hop-rates and the initial configuration where the solutions for z 0 get simplified significantly.Specifically, we would consider ρ = 0 for the initial configuration (see Sec. Below we discuss two special cases.

D.1 Case I:
A particularly simple solution can be acquired for the choice p 1 = p 2 = 1 ̸ = ε.As a result, Eq. (D.1 ) is reduced to a quadratic equation which leads to the following solution Note that in this case, the fugacity does not depend explicitly on the flip rates w 12 and w 21 .

D.2 Case II:
A comparatively cumbersome yet closed form solution is attained for the case p 1 = 1/2, p 2 = ε = w 21 = 1, ρ + = 1/4.Here the fugacity would be a function of w 12 and ρ 0 i.e. z 0 (w 12 , ρ 0 ).The reason behind keeping w 12 and ρ 0 as free parameters is that, the observables in the main text have been mostly analyzed as functions of these two parameters.The

E Block-diagonal structure of the transition rate matrix
In this appendix, we show the block-diagonal structure of the transition rate matrix M [Eq.(3)], reflecting the non-ergodicity of µ-ASEP-IAF.To illustrate this with an example for µ = 2, we consider a small system of size L = 4 where the number of impurity and vacancy are given by N + = 1 and N 0 = 1, respectively, and the total number of species 1 and species 2 particles is N 1 + N 2 = 2.Total number of configurations in the configuration space, in this case, is 48.However, since there is no spatial disorder in the transition rates, we take into account the translational invariance of the model on a periodic lattice.Consequently, there are 12 independent configurations of the system, which we denote as follows (depending on the sequence of species 1 and 2) Note that, in absence of the flip dynamics (i.e.w 12 = w 21 = 0), sectors I become disconnected from I I, similarly I I I gets disconnected from I V , resulting in four blocks in the transition matrix.On the other hand, in the special limit when N 1 + N 2 = 1, we would have a single block with the system becoming ergodic.
Next we explore the variation in the number of blocks as the system size is increased.We keep N 0 = 1 throughout, because it appears that the number of blocks depends on the arrangements of 1, 2 and +, but not on the location of vacancies.This might be better understood in a box-particle representation of the model where 1, 2, + denote boxes and 0-s are particles.
For L = 5, the special case N 1 + N 2 = 1 (N + = 3) keeps the system ergodic with a single block only.But, as we increase N 1 + N 2 , e.g.N + = 2 and N 1 + N 2 = 2, one can check that the rate matrix is block-diagonal with 3 blocks.With further increase in N 1 + N 2 (=3) which also corresponds to N + = 1, we have 4 blocks in the transition rate matrix.Below we present N blocks in a tabular form, explicitly for a few sets of (L, N + ), with N 0 = 1 and In fact, for fixed system size L, with N 0 = 1, the general formulae for number of blocks N blocks in the transition rate matrix, for cases N + = 1 and N + = 2 turn out to be

Figure 1 :
Figure 1: The figure illustrates all possible microscopic dynamical processes for the 2-TASEP-IAF with µ = 2.The species 1 and species 2 particles can hop to right (if vacant) with rates p 1 and p 2 whereas the corresponding hopping rate for the impurity (+) is ε.The flip process between the species 1 and 2 can occur only in the presence of an impurity at the right neighbor.The corresponding flip rates are w 12 (for 1 transforming to 2) and w 21 (for 2 transforming to 1).

Figure 2 :
Figure 2: The figure illustrates an example of µ = 4 case where certain flips are activated between pairs of species by the impurity, whereas some of the flips are absent e.g.w 24 = 0.In this scenario, α K I would be different for different pairs (I, K).
1 in the steady state, with ρ I being the average density of the non-conserved species I. Proceeding in the same way as shown in cases of µ = 2 and µ = 3, we obtain the partition function for the general µ-TASEP-IAF to be

Figure 5 :
Figure5:The figure shows comparison of drift currents J I0 between theory (solid lines) and Monte Carlo simulations (points) by changing flip rate w 12 .The nonlinear monotonic decrease and increase in J 10 and J 20 respectively are similar in nature to the effect of w 12 on the species densities (Fig.3).The impurity drift current increases slowly with increasing w 12 .The parameters used are L = 10 3 , p 1 = 0.3, p 2 = 1.0, ε = 0.1, w 21 = 0.6, ρ + = 0.2, ρ 0 = 0.2.The ensemble average is done over 10 5 samples.

Figure 6 :
Figure6: The figure illustrates the variation of two-point correlation C 01 between the vacancy and species 1 by tuning the vacancy density ρ 0 .We observe that this correlation initially starts increasing with increasing ρ 0 , reaches to a local maximum, then decreases.However, as ρ 0 is increased further, C 01 reaches to a local minimum and then starts increasing again.The parameters used are L = 10 3 , p 1 = 0.3, p 2 = 1.0, ε = 0.1, w 12 = 0.8, w 21 = 1.0, ρ + = 0.2.The ensemble average is done over 10 5 samples.

1 1 = w 21 and d 1 , 1 2
), we observe that the matrices (D I ) corresponding to the species (I) are upper triangular.The subscript I in matrix element d i, j I denotes the species I whereas the superscript (i, j) refers to the i-th row and j-th column of the matrix.The notation (m) r used in the expression of d m,m+r I corresponds to the Pochhammer symbol for rising factorials, (m) r = m(m + 1)(m + 2) . . .(m + r − 1) with (m) 0 = 1.The matrix E representing vacancy is a lower shift matrix and the matrix A representing impurity has non-zero terms in a single row only.For the simpler case µ = 2 the only changes in comparison to Eq. (48) will be in the values of d 1,1 I (I = 1, 2), which would be simply d 1,= w 12 .In fact, the matrix representations [Eq.(48)] for the 3-ASEP-IAF can be generalized for any µ > 0 in a straightforward manner.The representations will remain the same, only the values of d 1,1 I (I = 1, 2, . . ., µ) would change where d 1,1 I

Figure 10 :
Figure 10: The figure illustrates the negative differential mobility of the currents in the 3-ASEP-IAF.In (a), the drift current of the species 1 (driven by bias b) monotonically decreases with increasing bias.The current of the other driven mode (driven due to unidirectional motion with rate ε) impurity (+) shows non-monotonic behavior, it increases initially but ultimately decreases with increasing bias for large values of b.In (b), the drift current of both non-driven modes (species 2 and species 3) decrease with increasing bias b for large values of the bias.The figure (c) shows that the flip current between any pair of species also exhibits negative differential mobility with increasing b.Although the flip current is not related to the drift bias b directly, still it decreases with increasing bias for large values of b.The parameters used here are L = 10 3 , ε = 0.1 , p 3 = q 3 = 1.0 , w 31 = w 12 = 0.8 , w 13 = w 32 = 0.2 , w 21 = w 23 = 0.5 , ρ + = 0.3 , ρ 0 = 0.4 .The ensemble average is done over 10 6 samples.

Figure 11 :
Figure 11: The figure illustrates a two lane totally asymmetric simple exclusion process and the identification of its components to the equivalent constituents of the 2-TASEP-IAF.

Figure 12 :
Figure 12: The figure shows the connection between each microscopic dynamics of the two lane TASEP with the equivalent microscopic dynamics in the 2-TASEP-IAF.

Figure 13 :
Figure 13: The figure shows the variation of the total drift current of species 1 and 2 with the flip rate w 12 .While the flip rate (equivalent lane change rate in two lane TASEP) is varied over a large range w 12 ∈ (0, 1), the corresponding increase in the total drift current (total flow along the two lanes in TASEP) is reasonably small.The parameters used are L = 10 4 , p 1 = 0.3 , p 2 = 1.0 , ε = 0.1 , w 21 = 1.0 , ρ + = 0.2 , ρ 0 = 0.2 .The ensemble average is done over 10 7 samples.

Figure 15 :
Figure 15: The figure illustrates the connection between the simplified form of the enzymatic chemical reaction Eq.(C.2 ) (in a narrow channel with many units of drifting enzymes, substrates and products) and the 2-TASEP-IAF.The enzyme, substrate and product in the chemical system can be identified as impurity, species 1 and species 2 particle in the 2-TASEP-IAF.

11 + 0 ⃝ is 6 ×
≡ I 1 , 110+ ≡ I 2 , 101+ ≡ I 3 , 12 + 0 ≡ I I 1 , 120+ ≡ I I 2 , 102+ ≡ I I 3 , 21 + 0 ≡ I I I 1 , 210+ ≡ I I I 2 , 201+ ≡ I I I 3 , 22+ 0 ≡ I V 1 , 220+ ≡ I V 2 , 202+ ≡ I V 3 .(E.1)We have divided the 12 configurations in Eq. (E.1 ) into 4 sectors I, I I, I I I, I V where the three configurations within a given sector are connected through the drift dynamics.To investigate the connectivity between these sectors through the flip dynamics, below we provide the full transition rate matrix for these 12 configurations (enumerated consecutively from I 1 to I V 3 ), 6 null matrix.In Eqs.(E.2 ) and (E.3 ), we clearly observe that the transition rate matrix in in block-diagonal form with two blocks.We observe that sector I is connected to sector I I through flip dynamics, whereas sector I I I and I V are also connected to each other via flip dynamics.However, sectors (I, I I) are disconnected from sectors (I I I, I V ), thereby creating two separate blocks in the rate matrix.