Multi species asymmetric simple exclusion process with impurity activated
flips

Amit Kumar Chatterjee1*, Hisao Hayakawa1,2†,

1 Yukawa Institute for Theoretical Physics, Kyoto University, Sakyo-ku, Kyoto
606-8502, Japan

2 Center for Gravitational Physics and Quantum Information, Yukawa Institute
for Theoretical Physics, Kyoto University Sakyo-ku, Kyoto 606-8502, Japan

*<EMAIL_ADDRESS>${\dagger}$<EMAIL_ADDRESS>
## Abstract

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.

###### Contents

  1. 1 Introduction
  2. 2 Model: $\mu$-TASEP-IAF
    1. 2.1 Microscopic dynamics
    2. 2.2 Steady state: matrix product ansatz
  3. 3 Partition function for special initial configuration
    1. 3.1 $\mu=2$
    2. 3.2 $\mu=3$
    3. 3.3 general $\mu$
  4. 4 Observables: comparisons of analytical results with Monte Carlo simulations
    1. 4.1 Species densities
    2. 4.2 Drift current
    3. 4.3 Two-point correlations
    4. 4.4 Flip current
    5. 4.5 Non-ergodicity: dependence on initial configuration
  5. 5 Partially asymmetric generalization: $\mu$-ASEP-IAF
    1. 5.1 Matrix algebra, auxiliaries and matrix representations
    2. 5.2 Partition function for special initial configuration
    3. 5.3 Species densities, drift current and flip current
    4. 5.4 Negative differential mobility
  6. 6 Summary and future directions
  7. A Connection between $\mu$-TASEP-IAF and multi lane TASEP
  8. B A variation of $\mu$-TASEP-IAF, connection to multi-lane traffic flow
  9. C Connection between $\mu$-TASEP-IAF and enzymatic chemical reactions
  10. D Solution of Eq. (75) for fugacity $z_{0}$: special cases
    1. D.1 Case I:
    2. D.2 Case II:
  11. E Block-diagonal structure of the transition rate matrix

##  1 Introduction

Non-equilibrium stochastic processes are ubiquitous in nature, with wide range
of applicability in physics [1, 2, 3], chemistry [1], biology [4] and
interdisciplinary areas [5, 6, 7]. In fact, even in one dimension, several
models of non-equilibrium statistical mechanics exhibit surprisingly rich
physical phenomena including phase transitions [8] along with the important
feature of analytical tractability [9]. The asymmetric simple exclusion
process (ASEP) [10, 11, 2, 12] is broadly regarded as a paradigmatic model for
non-equilibrium transport processes as diverse as traffic and pedestrian flow
[13], mRNA translation by ribosomes [14] and motor protein transport through
single filaments [15, 16] etc.

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) [21, 12] 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
$\frac{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 $\frac{5}{2}$ in
comparison to the usual KPZ exponent $\frac{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 [58, 57, 59, 60, 61, 62, 63].

In this article, we study a class of the multi-species $(I=1,2,\dots,\mu)$
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 ($\mu$-ASEP-IAF). Note that the total number of impurities, along with
that of the vacancies, remain conserved in the $\mu$-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,\dots,\mu)$ 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 $\mu$-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 $\mu$-ASEP-IAF are as follows. (i) We aim to
obtain an exact non-equilibrium steady state of the $\mu$-ASEP-IAF under
periodic boundary condition, so that it would be an important addition to the
category of exact solvable models in disordered systems. (ii) The $\mu$-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 $\mu$-TASEP-IAF can be
mapped to multi-lane TASEP which is a basic model for multi-lane traffic flow.
Different species of particles in $\mu$-ASEP-IAF play the roles of particles
in different lanes of multi-lane TASEP and the impurities in $\mu$-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 $\mu$-ASEP-IAF can be thought as
an enzymatic chemical reaction like $S+E\rightarrow P+E$, which is a crude
approximation of the Michaelis-Menten reaction scheme $S+E\rightleftharpoons
SE\rightarrow 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 $\mu$-TASEP-IAF (totally asymmetric
hopping) and $\mu$-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 $\mu$-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 $\mu$-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 $\mu$-TASEP-IAF
is generalized to $\mu$-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 $\mu$-TASEP-IAF
and multi lane TASEP in Appendix A. A variation of $\mu$-TASEP-IAF that comes
up with better features in connection to traffic in multi-lane problems is
discussed in Appendix B. The connections between $\mu$-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.

##  2 Model: $\mu$-TASEP-IAF

###  2.1 Microscopic dynamics

Let us consider a system of $\mu$ different species of particles and
impurities on a one dimensional periodic lattice with $L$ sites
$i=1,2,\dots,L$. Each site can either be vacant or it can be occupied by only
one particle of any of the species $I=1,2,\dots,\mu$ 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,

$\displaystyle\mathrm{drift\,(species):}\hskip 14.22636ptI0\,\,$
$\displaystyle\stackrel{{\scriptstyle p_{I}}}{{\longrightarrow}}$
$\displaystyle\,\,0I,\hskip 22.76228ptI=1,2,...,\mu,$
$\displaystyle\mathrm{drift\,(impurity):}\hskip 8.5359pt+0,\,\,$
$\displaystyle\stackrel{{\scriptstyle\epsilon}}{{\longrightarrow}}$
$\displaystyle\,\,0+$ $\displaystyle\mathrm{flip:}\hskip 11.38092ptI+\,\,$
$\displaystyle\xrightleftharpoons[w_{KI}]{w_{IK}}$ $\displaystyle\,\,K+,\hskip
17.07182ptI,K=1,...,\mu.$ (1)

Figure 1: The figure illustrates all possible microscopic dynamical processes
for the $2$-TASEP-IAF with $\mu=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 $\epsilon$. 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$).

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 $\epsilon$. 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_{IK}$ and the
reverse transformation from species $K$ to $I$ occurs with rate $w_{KI}$.
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 $\mu$-TASEP-IAF is
$(p_{I},\epsilon,w_{IK},\rho_{+},\rho_{0})$, where $\rho_{+}=N_{+}/L$ and
$\rho_{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 $\mu=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 $\left\\{\dots 0+102011+2\dots\right\\}$. 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
$\left\\{\dots 0+10201{\textit{1}}+2\dots\right\\}$, it can transform into 2
by the $+$ at its right neighbor, thereby changing the configuration to
$\left\\{\dots 0+10201{\textit{2}}+2\dots\right\\}$. However, another 1 in the
initial configuration $\left\\{\dots 0+1020{\textit{1}}1+2\dots\right\\}$ 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 $\left\\{\dots 0+1020{\textit{2}}1+2\dots\right\\}$ is
never accessible.

###  2.2 Steady state: matrix product ansatz

Any configuration of the $\mu$-TASEP-IAF can be represented by
$\left\\{s_{i}\right\\}\equiv\left\\{s_{1},s_{2},\dots,s_{L}\right\\}$, where
$s_{i}$ denotes the occupation at site $i$. Clearly, $s_{i}$ can be one of the
species $K=1,2,\dots,\mu$ 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

$\displaystyle
P(\left\\{s_{i}\right\\})\propto\mathrm{Tr}\left[\prod_{i=1}^{L}X_{i}\right],$
(2) $\displaystyle
X_{i}=E\,\delta_{s_{i},0}+A\,\delta_{s_{i},+}+\sum_{K=1}^{\mu}D_{K}\,\delta_{s_{i},K}.$
(3)

In Eq. (3), any configuration $\left\\{s_{i}\right\\}$ is represented by a
string of matrices $\left\\{X_{i}\right\\}$ 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 $\mu$-TASEP-IAF
is dictated by the Master equation

$\frac{d}{dt}|P(t)\rangle=M|P(t)\rangle,$ (4)

which in steady state becomes $M|P\rangle=0$. Here $|P\rangle$ 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

$M=\sum_{i=1}^{L}\left(I\otimes\dots I\otimes\mathcal{M}_{i,i+1}\otimes
I\dots\otimes I\right),$ (5)

where $\mathcal{M}_{i,i+1}$ is a $(\mu+2)^{2}\times(\mu+2)^{2}$ dimensional
matrix and $I$ is $(\mu+2)\times(\mu+2)$ dimensional identity matrix placed at
every site except the pair $(i,i+1)$. Then the steady state $M|P\rangle=0$ of
the $\mu$-TASEP-IAF can be achieved through the following two-site flux
(probability current) cancellation condition

$\displaystyle\mathcal{M}_{i,i+1}\mathbf{X}_{i}\otimes\mathbf{X}_{i+1}$
$\displaystyle=$
$\displaystyle\tilde{\mathbf{X}}_{i}\otimes\mathbf{X}_{i+1}-\mathbf{X}_{i}\otimes\tilde{\mathbf{X}}_{i+1},$
(6)

where

$\displaystyle\mathbf{X}=\left(E,A,D_{1},D_{2}\dots,D_{\mu}\right)^{T},$ (7)

and

$\displaystyle\tilde{\mathbf{X}}=\left(\tilde{E},\tilde{A},\tilde{D}_{1},\tilde{D}_{2}\dots,\tilde{D}_{\mu}\right)^{T},$
(8)

where $(.)^{T}$ denotes the transpose of the row vector $(.)$ and
$\tilde{E},\tilde{A},\tilde{D}_{K}$ 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,\dots,\mu$). We find that suitable choices for the auxiliary matrices
for the $\mu$-TASEP-IAF are

$\displaystyle\tilde{E}=1,\,\,\tilde{A}=0,\,\,\tilde{D}_{K}=0\hskip
14.22636ptK=1,2,\dots,\mu.$ (9)

Correspondingly, the matrices $E,A$ and $D_{K}$ have to obey the matrix
algebra consisting of the equations given below

$\displaystyle p_{K}D_{K}E$ $\displaystyle=$ $\displaystyle D_{K},\hskip
64.01869ptK=1,\dots,\mu$ (10) $\displaystyle\epsilon AE$ $\displaystyle=$
$\displaystyle A,$ (11) $\displaystyle\sum_{\begin{subarray}{c}I=1\\\ I\neq
K\end{subarray}}^{\mu}w_{IK}D_{I}A$ $\displaystyle=$ $\displaystyle
D_{K}A\sum_{\begin{subarray}{c}I=1\\\ I\neq
K\end{subarray}}^{\mu}w_{KI},\hskip 14.22636ptK=1,\dots,\mu.$ (12)

The last relation in Eq. (12) 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. (12) allows scalar
solutions when the hopping rates for every species and the impurity become
equal i.e. $p_{K}=\epsilon$ 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. (12). Below we discuss the cases $\mu=2$,
$\mu=3$ extensively with explicit matrix representations and then generalize
them to get the matrix representations for general $\mu$.

$\mu=2:$ For the $2$-TASEP-IAF ($K=1,2$), the matrix algebra [Eq. (12)]
simplifies to

$\displaystyle p_{1}D_{1}E$ $\displaystyle=$ $\displaystyle D_{1},\hskip
14.22636ptp_{2}D_{2}E=D_{2},$ (13) $\displaystyle\epsilon AE$ $\displaystyle=$
$\displaystyle A,$ (14) $\displaystyle w_{12}D_{1}A$ $\displaystyle=$
$\displaystyle w_{21}D_{2}A.$ (15)

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
$\mu=2$. However, there are non-zero drift currents in the system. We find the
following $3\times 3$ matrix representations that satisfy the matrix algebra
in Eq. (15),

$\displaystyle D_{1}=w_{21}\left(\begin{array}[]{ccc}1&1&0\\\ 0&0&0\\\
0&0&0\\\ \end{array}\right),\hskip
14.22636ptD_{2}=w_{12}\left(\begin{array}[]{ccc}1&0&1\\\ 0&0&0\\\ 0&0&0\\\
\end{array}\right),$ (22) $\displaystyle
E=\left(\begin{array}[]{ccc}\frac{1}{\epsilon}&0&0\\\
\frac{1}{p_{1}}-\frac{1}{\epsilon}&\frac{1}{p_{1}}&0\\\
\frac{1}{p_{2}}-\frac{1}{\epsilon}&0&\frac{1}{p_{2}}\\\
\end{array}\right),\,\,\,A=\left(\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\
0&0&0\\\ \end{array}\right).$ (29)

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].

$\mu=3:$ The matrix algebra in Eq. (12) for the $3$-TASEP-IAF process reads as

$\displaystyle p_{1}D_{1}E=D_{1},\hskip 8.5359ptp_{2}D_{2}E$ $\displaystyle=$
$\displaystyle D_{2},\hskip 8.5359ptp_{3}D_{3}E=D_{3},$ (30)
$\displaystyle\epsilon AE$ $\displaystyle=$ $\displaystyle A,$ (31)
$\displaystyle w_{21}D_{2}A+w_{31}D_{3}A$ $\displaystyle=$
$\displaystyle(w_{12}+w_{13})D_{1}A,$ (32) $\displaystyle
w_{12}D_{1}A+w_{32}D_{3}A$ $\displaystyle=$
$\displaystyle(w_{21}+w_{23})D_{2}A,$ (33) $\displaystyle
w_{13}D_{1}A+w_{23}D_{2}A$ $\displaystyle=$
$\displaystyle(w_{31}+w_{32})D_{3}A.$ (34)

In comparison to the last relation in Eq. (15), clearly the flip processes for
the three-species case given by the last three relations in Eq. (34) 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. (34) is

$\displaystyle
w_{12}D_{1}A-w_{21}D_{2}A=w_{23}D_{2}A-w_{32}D_{3}A=w_{31}D_{3}A-w_{13}D_{1}A.$
(35)

We obtain the following $4\times 4$ representations of the matrices that
satisfy Eq. (34) along with Eq. (35),

$\displaystyle D_{1}=d_{1}\left(\begin{array}[]{cccc}1&1&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\
\end{array}\right),D_{2}=d_{2}\left(\begin{array}[]{cccc}1&0&1&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ \end{array}\right),$ (44) $\displaystyle
D_{3}=d_{3}\left(\begin{array}[]{cccc}1&0&0&1\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ \end{array}\right),$ (49) $\displaystyle
E=\left(\begin{array}[]{cccc}\frac{1}{\epsilon}&0&0&0\\\
\frac{1}{p_{1}}-\frac{1}{\epsilon}&\frac{1}{p_{1}}&0&0\\\
\frac{1}{p_{2}}-\frac{1}{\epsilon}&0&\frac{1}{p_{2}}&0\\\
\frac{1}{p_{3}}-\frac{1}{\epsilon}&0&0&\frac{1}{p_{3}}\\\
\end{array}\right),A=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ \end{array}\right),$ (58)

where

$\displaystyle d_{1}$ $\displaystyle=$ $\displaystyle
w_{21}w_{31}+w_{23}w_{31}+w_{32}w_{21},$ (59) $\displaystyle d_{2}$
$\displaystyle=$ $\displaystyle w_{12}w_{32}+w_{13}w_{32}+w_{31}w_{12},$ (60)
$\displaystyle d_{3}$ $\displaystyle=$ $\displaystyle
w_{13}w_{23}+w_{12}w_{23}+w_{21}w_{13}.$ (61)

We note that the condition Eq. (35), with the explicit matrix representations
from Eq. (58), becomes

$\displaystyle w_{KI}D_{K}A-w_{IK}D_{I}A=\alpha A,$
$\displaystyle\alpha=w_{12}w_{23}w_{31}-w_{21}w_{13}w_{32}.$ (62)

The parameter $\alpha$ in Eq. (62) 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 $\alpha$. In fact $\alpha=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_{KI}D_{K}A=w_{IK}D_{I}A$.

general $\mu:$ For the general case of $\mu$-TASEP-IAF ($K=1,2,\dots,\mu$),
the structures of the matrices $D_{K},E,A$ would be similar to that of Eq.
(29) and Eq. (58). More precisely, the matrices are $(\mu+1)\times(\mu+1)$
dimensional and the corresponding explicit representations of the matrices are
given by

$\displaystyle D_{K}$ $\displaystyle=$ $\displaystyle
d_{K}|1\rangle\left(\langle 1|+\langle K|\right),\,\,\,\,\,\,K=1,\dots,\mu,$
(63) $\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon}|1\rangle\langle
1|+\sum_{K=2}^{\mu+1}\frac{1}{p_{K-1}}|K\rangle\langle
K|+\left(\frac{1}{p_{K-1}}-\frac{1}{\epsilon}\right)|K\rangle\langle 1|,$ (64)
$\displaystyle A$ $\displaystyle=$ $\displaystyle|1\rangle\langle 1|.$ (65)

In Eq. (65), the vector $\langle I|=(0,\dots,0,1,0,\dots 0)$ where $1$ is
placed at the $I$-th element with all other elements being zero and
$|I\rangle$ is the transpose of $\langle I|$. Notably, the values of the
coefficients $d_{K}$ associated with matrices $D_{K}$ in Eq. (65), can be
calculated by solving the set of $\mu$ homogeneous linear equations of the
form

$d_{K}\sum_{I\neq K}w_{KI}+\sum_{I\neq K}w_{IK}d_{I}=0,\,\,\,K=1,2,\dots,\mu.$
(66)

The summations over the index $I$ in Eq. (66) generally includes all
$I=1,\dots,\mu$ 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 $\alpha$ [Eq. (62)] which dictates the flip-current between any two
species, will be given by
$\alpha=(\prod\limits_{K=1}^{\mu}w_{KK+1}-\prod\limits_{K=1}^{\mu}w_{K+1K})$.
Notably, this reduces to the general solution of $\mu=3$ as given in Eq. (62).
However, unlike the cases of $\mu=2$ or $\mu=3$ or the special instance of the
multi-species case stated above which allow $\alpha$ to be independent of the
species $I$, in general $\alpha$ would depend on the particulars of the pairs
$(I,K)$ for $\mu>3$. Mathematically, Eq. (62) would be generalized as follows

$w_{KI}D_{K}A-w_{IK}D_{I}A=\alpha_{KI}A,\hskip
14.22636ptI,K=1,2,3,4,\dots,\mu,$ (67)

where $\alpha_{KI}$ have different values for different pairs $(I,K)$. As an
example, one can consider $\mu=4$ with the allowed set of flip dynamics
described in Fig. 2.

Figure 2: The figure illustrates an example of $\mu=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, $\alpha_{KI}$ would be
different for different pairs $(I,K)$.

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 $\alpha_{23}=0$ whereas for other pairs
$\alpha_{KI}\neq 0$.

##  3 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 $\mu$-TASEP-IAF in the usual form
of $\mathrm{Tr}[T^{L}]$, even under periodic boundary conditions. Here the
“transfer matrix” $T$ refers to

$T=z_{0}E+z_{+}A+\sum_{K}D_{K},$ (68)

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.
$\left\\{..11..11+..\right\\}\equiv\left\\{..1^{n+1}+..\right\\}$. At any step
in the $\mu$-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
$\mu$-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.

###  3.1 $\mu=2$

One special initial configuration $C(t=0)\equiv C(0)$ (represented by
matrices) for the $2$-TASEP-IAF is

$\displaystyle C(0)=\underbrace{D_{1}A\dots D_{1}A}\,\,\underbrace{D_{2}A\dots
D_{2}A}\,\,\underbrace{D_{2}\dots D_{2}}\,\,\underbrace{D_{1}\dots
D_{1}}\,\,\underbrace{E\dots E},$ (69)

where $\underbrace{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 $\bar{\rho}=\bar{N}/L$. Further, we have taken
the two sequences of $D_{1}A$ and $D_{2}A$ to be of equal density $\rho_{+}$
so that the density of impurity ($A$) in each of these sequences is
$\rho_{+}/2$. This ensures that the total density of impurities is $\rho_{+}$.
The initial configuration in Eq. (69) satisfies the relation
$\bar{\rho}=\frac{1}{2}(1-\rho_{0}-2\rho_{+})$, where $\rho_{0}$ and
$\rho_{+}$ are densities of the vacancies and impurities respectively. In the
steady state we have $\rho_{0}+\rho_{+}+\rho_{1}+\rho_{2}=1$, with $\rho_{1}$
and $\rho_{2}$ being the average densities of species $1$ and species $2$
particles in the steady state. We emphasize that $\rho_{0}$ and $\rho_{+}$ are
input parameters while $\rho_{1}$ and $\rho_{2}$ are derived quantities.
Starting from Eq. (69), any accessible configuration $C_{ss}$ in the steady
state is of the following generic form

$\displaystyle C_{ss}=\prod_{k=1}^{N_{+}}\left(\tau
D_{1}+(1-\tau)D_{2}\right)E^{m_{k}}AE^{n_{k}}\,\prod_{i=1}^{\bar{N}}D_{2}E^{r_{i}}\,\prod_{j=1}^{\bar{N}}D_{1}E^{s_{j}},$
(70)

subjected to the constraint
$\sum\limits_{i=1}^{N_{+}}(m_{i}+n_{i})+\sum\limits_{j=1}^{\bar{N}}(r_{j}+s_{j})=N_{0}$.
The parameter $\tau$ in Eq. (70) can take value either $1$ or $0$. The
partition function, obeying the above mentioned constraint, is given by

$\displaystyle
Q_{N_{0},N_{+}}=\sum_{\left\\{n_{i}\right\\}}\,\sum_{\left\\{m_{i}\right\\}}\,\sum_{\left\\{r_{j}\right\\}}\,\sum_{\left\\{s_{j}\right\\}}\,\,\mathrm{Tr}\left[\prod_{k=1}^{N_{+}}(D_{1}+D_{2})E^{m_{k}}AE^{n_{k}}\prod_{i=1}^{\bar{N}}D_{2}E^{r_{i}}\prod_{j=1}^{\bar{N}}D_{1}E^{s_{j}}\right]$
(71) $\displaystyle\hskip
133.72786pt\times\,\delta(\sum_{i=1}^{N_{+}}(n_{i}+m_{i})+\sum_{j=1}^{\bar{N}}(r_{j}+s_{j})-N_{0}).$
(72)

It would be useful to get rid of the $\delta(.)$ 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. (15) along with the explicit
representations from Eq. (29). It should be mentioned that the projector form
of $A$ [Eq. (29)] leads to factorization of the matrix strings, e.g.

$\dots AD_{1}ED_{2}EAD_{2}EEA\dots=\dots|1\rangle\,\,\langle
1|D_{1}ED_{2}E|1\rangle\,\,\langle 1|D_{2}EE|1\rangle\,\,\langle 1|\dots,$
(73)

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

$\displaystyle
Q_{N_{+}}(z_{0})=\left(\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]\left(\frac{1}{1-\frac{z_{0}}{\epsilon}}\right)\right)^{N_{+}}\left(\frac{w_{21}w_{12}}{(1-\frac{z_{0}}{p_{1}})(1-\frac{z_{0}}{p_{2}})}\right)^{\bar{N}}.$
(74)

For the special initial configuration in Eq. (69), we have derived the
partition function in Eq. (74). The fugacity $z_{0}$ can be obtained as a
function of the vacancy density and other input parameters by inverting the
density-fugacity relation

$\rho_{0}=\frac{z_{0}}{L}\frac{d}{dz_{0}}\mathrm{ln}(Q_{N_{+}}(z_{0})).$ (75)

In general, the solution of $z_{0}$ obtained from Eq. (75), 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 $\rho_{+}=N_{+}/L$. The expression
Eq. (74) would be used for evaluating the observables of interest in the next
Sec. 4.

###  3.2 $\mu=3$

Similar to the initial configuration $C(t=0)$ for $\mu=2$, a suitable initial
configuration for the three species case that enables us to perform analytical
calculation of partition function and observables, is

$\displaystyle C(0)\equiv\underbrace{D_{1}A\dots
D_{1}A}\,\,\underbrace{D_{2}A\dots D_{2}A}\,\,\underbrace{D_{3}A\dots
D_{3}A}\underbrace{D_{3}\dots D_{3}}\,\,\underbrace{D_{2}\dots
D_{2}}\,\,\underbrace{D_{1}\dots D_{1}}\,\,\underbrace{E\dots E}.$ (76)

The density of the uninterrupted sequence of each species $1,2,3$ in Eq. (76)
is taken to be equal to $\bar{\rho}=\bar{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\rho_{+}/3$ where the density of impurities in each of these sequences is
$\rho_{+}/3$. Consequently, the choice of initial configuration in Eq. (76)
ensures that the total density of impurities remain $\rho_{+}$. In the steady
state, $\rho_{0}+\rho_{+}+\sum\limits_{I=1}^{3}\rho_{I}=1$, where $\rho_{I}$
is the average density of the non-conserved species $I$. Starting from Eq.
(76), the form of any accessible configuration $C_{ss}$ in steady state would
be

$\displaystyle
C_{ss}\equiv\prod_{k=1}^{N_{+}}\left(D_{1}\delta_{\tau,1}+D_{2}\delta_{\tau,2}+D_{3}\delta_{\tau,3}\right)E^{m_{k}}AE^{n_{k}}\prod_{i=1}^{\bar{N}}D_{3}E^{l_{i}}\,\prod_{i=1}^{\bar{N}}D_{2}E^{r_{i}}\,\prod_{i=1}^{\bar{N}}D_{1}E^{s_{i}},$
(77)

with the conservation of the total number of vacancies
$N_{0}=\sum\limits_{k=1}^{N_{+}}(m_{k}+n_{k})+\sum\limits_{i=1}^{\bar{N}}(l_{i}+r_{i}+s_{i})$,
where $\delta_{\tau,K}$ is the Kronecker delta symbol with $K=1,2,3$. As done
in case of $\mu=2,$ here also we associate a fugacity $z_{0}$ with the
vacancy. Using the matrix algebra from Eq. (34) alongside the matrix
representations in Eqs.(58) and (61), the partition function in grand
canonical ensemble becomes

$\displaystyle
Q_{N_{+}}(z_{0})=\left(\left[\sum_{I=1}^{3}\frac{d_{I}}{1-\frac{z_{0}}{p_{I}}}\right]\left(\frac{1}{1-\frac{z_{0}}{\epsilon}}\right)\right)^{N_{+}}\left(\prod_{I=1}^{3}\frac{d_{I}}{1-\frac{z_{0}}{p_{I}}}\right)^{\bar{N}},$
(78) (79)

where the explicit expressions for $d_{I}$-s have been presented earlier in
Eq. (61).

###  3.3 general $\mu$

For the $\mu$-TASEP-IAF ($K=1,\dots,\mu$), the generalization of the initial
configurations in Eq. (69) ($\mu=2$) and Eq. (76) ($\mu=3$) would be

$\displaystyle
C(0)\equiv\left(\prod_{i=1}^{N_{+}/\mu}D_{1}A\prod_{i=1}^{N_{+}/\mu}D_{2}A\dots\prod_{i=1}^{N_{+}/\mu}D_{\mu}A\right)\,\left(\prod_{i=1}^{\bar{N}}D_{1}\prod_{i=1}^{\bar{N}}D_{2}\dots\prod_{i=1}^{\bar{N}}D_{\mu}\right)\,\prod_{i=1}^{N_{0}}E.$
(80)

The above initial configuration is chosen in a way that the density of each
sequence $D_{I}A$ ($I=1,\dots,\mu$) is $2\rho_{+}/\mu$ in which the density of
impurities is equal to $\rho_{+}/\mu$, so that the total impurity density adds
up to $\rho_{+}$. We have
$\rho_{0}+\rho_{+}+\sum\limits_{I=1}^{\mu}\rho_{I}=1$ in the steady state,
with $\rho_{I}$ being the average density of the non-conserved species $I$.
Proceeding in the same way as shown in cases of $\mu=2$ and $\mu=3$, we obtain
the partition function for the general $\mu$-TASEP-IAF to be

$\displaystyle
Q_{N_{+}}(z_{0})=\left(\left[\sum_{I=1}^{\mu}\frac{d_{I}}{1-\frac{z_{0}}{p_{I}}}\right]\left(\frac{1}{1-\frac{z_{0}}{\epsilon}}\right)\right)^{N_{+}}\left(\prod_{I=1}^{\mu}\frac{d_{I}}{1-\frac{z_{0}}{p_{I}}}\right)^{\bar{N}},$
(81) (82)

where $d_{I}$ is the solution of Eq. (66) 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 $\mu$-TASEP-
IAF with $\mu=2$, $\mu=3$ and general $\mu$ in Eq. (74), Eq. (79) and Eq.
(82), respectively, for specific initial configuration Eq. (69), Eq. (76) and
Eq. (80), 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.

##  4 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 $\rho_{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\leftrightarrow 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 $\mu=2,$ except using $\mu=3$ for the case of average
flip current (since there is no net flip current between pair of species for
$\mu=2$). In particular, starting from the special initial configuration Eq.
(69) for $\mu=2$ (or, Eq. (76) for $\mu=3$), we will show agreements between
the analytical calculations and the Monte Carlo simulation results.

(a) (b)

Figure 3: The figures (a) and (b) show the variations of the average species
densities $\rho_{I}$ ($I=1,2$) against vacancy density $\rho_{0}$ and flip
rate $w_{12}$ respectively in the steady state. The species densities decrease
linearly with increasing $\rho_{0}$ in (a). In (b), with increasing $w_{12}$,
$\rho_{1}$ and $\rho_{2}$ decrease and increase respectively, both in
nonlinear manners. Notably, in the parameter range
$w_{21}<w_{12}<w_{12}^{\star},$ we observe $\rho_{2}<\rho_{1}$ in spite of the
higher flip rate of transformation from species $1$ to species $2$. The common
parameters for both figures (a) and (b) are
$L=10^{3},p_{1}=0.3,p_{2}=1.0,\epsilon=0.1,\rho_{+}=0.2$. For (a),
$w_{12}=0.4$ and $w_{21}=1.0$. For (b), $\rho_{0}=0.2$ and $w_{21}=0.6$. The
ensemble average is done over $10^{5}$ samples.

###  4.1 Species densities

First we consider the average densities $(\rho_{I})$ of the non-conserved
species $I=1,2$. The formal expression for $\rho_{I}$ in the steady state
under the periodic boundary condition, can be written as

$\displaystyle\rho_{1}=\frac{1}{2}(1-\rho_{0}-2\rho_{+})+\frac{\rho_{+}}{Q_{N_{+}}}\sum_{n_{1}=0}^{\infty}..\sum_{n_{N_{+}}=0}^{\infty}\sum_{m_{1}=0}^{\infty}..\sum_{m_{N_{+}}=0}^{\infty}\sum_{r_{1}=0}^{\infty}..\sum_{r_{\bar{N}}=0}^{\infty}\sum_{s_{1}=0}^{\infty}..\sum_{s_{\bar{N}}=0}^{\infty}$
$\displaystyle\mathrm{Tr}\left[D_{1}(z_{0}E)^{m_{1}}A(z_{0}E)^{n_{1}}\prod_{k=2}^{N_{+}}(D_{1}+D_{2})(z_{0}E)^{m_{k}}A(z_{0}E)^{n_{k}}\prod_{k=1}^{\bar{N}}D_{2}(z_{0}E)^{r_{k}}\prod_{k=1}^{\bar{N}}D_{1}(z_{0}E)^{s_{k}}\right].$
(84)

To elaborate Eq. (84), the main point is to note the expression inside the
trace ($\mathrm{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. (70). From Eq. (70), to
arrive at the matrix string inside the trace in Eq. (84), one has to put
$\tau=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,\dots,N_{+}$, we have a combinatorial pre-factor
$\rho_{+}=N_{+}/L$ in Eq. (84). The summations over all the variables
$\left\\{m,n,r,s\right\\}$ from zero to infinity are performed as the system
is considered in the grand canonical ensemble. The first factor
$(1-\rho_{0}-2\rho_{+})/2$ in Eq. (84) is to take care of the $\bar{N}$ number
of $D_{1}$-s present in the initial configuration Eq. (69) that cannot flip.
Using the matrix algebra and matrix representations from Eqs. (15) and (29)
respectively, we finally arrive at the following expressions for the average
densities of the non-conserved species,

$\displaystyle\rho_{1}$ $\displaystyle=$
$\displaystyle\rho_{+}\,\frac{w_{21}}{\left(1-\frac{z_{0}}{p_{1}}\right)}\,\frac{1}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]}\,+\,\frac{1}{2}(1-2\rho_{+}-\rho_{0})$
(85) $\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle\rho_{+}\,\frac{w_{12}}{\left(1-\frac{z_{0}}{p_{2}}\right)}\,\frac{1}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]}\,+\,\frac{1}{2}(1-2\rho_{+}-\rho_{0}).$
(86)

The fugacity $z_{0}$ is obtained in terms of the input parameters by solving
the density-fugacity relation
$\rho_{0}=z_{0}\frac{d}{dz_{0}}\mathrm{ln}(Q_{N_{+}}(z_{0}))$ using
Mathematica. Replacing the solution of $z_{0}$ (which is too lengthy to
provide here) in Eq. (86), we finally get the average densities as functions
of input parameters $(p_{1,2},w_{12,21},\epsilon,\rho_{0},\rho_{+})$ only.

Figure 4: The figure illustrates the non-monotonic behaviors of the drift
currents $J_{I0}$ ($I=1,2$) with increasing vacancy density $\rho_{0}$, unlike
the drift current for the impurities $J_{+0}$ which increases monotonically
with increasing $\rho_{0}$. The analytical results (solid lines) are in
agreement with the Monte Carlo simulation results (points). The parameters
used are
$L=10^{3},p_{1}=0.3,p_{2}=1.0,\epsilon=0.1,w_{12}=0.4,w_{21}=1.0,\rho_{+}=0.2$.
The ensemble average is done over $10^{5}$ samples. Figure 5: 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,\epsilon=0.1,w_{21}=0.6,\rho_{+}=0.2,\rho_{0}=0.2$.
The ensemble average is done over $10^{5}$ samples.

We compare the analytical results with those of Monte Carlo simulations,
starting from the same initial configuration Eq. (69). 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. (69) i.e. simply by
tuning $\bar{\rho}$ which is related to vacancy density $\rho_{0}$ as
$\rho_{0}+2\bar{\rho}=1-2\rho_{+}$. In Figs. 3(a) and (b), we observe that the
analytical and simulation results are in agreement with each other, where
$\rho_{I}$ is plotted against $\rho_{0}$ and $w_{12}$ (flip rate of species
$1$ to species $2$), respectively. The species densities decrease linearly
with increasing $\rho_{0}$ [Fig. 3(a)] whereas they decrease in non-linear
fashion with increasing $w_{12}$ [Fig. 3(b)]. Notably, in the absence of any
drift, we would have $\rho_{2}=\rho_{1}$ exactly at $w_{12}^{\star}=w_{21}$.
However, due to the hopping process, this point shifts to

$w_{12}^{\star}=w_{21}\frac{(1-\frac{z_{0}}{p_{2}})}{(1-\frac{z_{0}}{p_{1}})}.$
(87)

Consequently in Fig. 3(b), for a particular set of chosen parameters, we
observe that for $w_{21}<w_{12}<w_{12}^{\star}$, one still has
$\rho_{2}<\rho_{1}$. In other words, when $w_{12}\in(w_{21},w_{12}^{\star}),$
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 $\mu$, the general expression
for the average density $\rho_{I}$ for the non-conserved species $I$
($I=1,2,\dots,\mu$) is

$\displaystyle\rho_{I}=\rho_{+}\,\frac{d_{I}}{\left(1-\frac{z_{0}}{p_{I}}\right)}\,\frac{1}{\sum\limits_{K=1}^{\mu}\frac{d_{K}}{\left(1-\frac{z_{0}}{p_{K}}\right)}}\,+\,\frac{1}{\mu}(1-2\rho_{+}-\rho_{0}),$
(88)

where $d_{I}$ is the solution of Eq. (66) (e.g. the solution for $\mu=3$ is
explicitly given in Eq. (61)).

###  4.2 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}\langle 10\rangle$, where $\langle
10\rangle$ is the ensemble average of the pair $10$. In terms of matrices the
expression $\langle 10\rangle$ simply translates to $\langle D_{1}E\rangle$.
The current $J_{10}$ can be calculated in two parts,

$\displaystyle J_{10}=p_{1}\langle 10\rangle=p_{1}\langle
D_{1}E\rangle=J_{10}^{(1)}+J_{10}^{(2)}=p_{1}\langle
D_{1}E\rangle^{(1)}+p_{1}\langle D_{1}E\rangle^{(2)},$ (89)

where $J_{10}^{(1)}=p_{1}\langle D_{1}E\rangle^{(1)}$ is the contribution from
the drift of species $1$ particles that can flip and
$J_{10}^{(2)}=p_{1}\langle D_{1}E\rangle^{(2)}$ is the corresponding
contribution from species $1$ particles that cannot flip (as they cannot have
any impurity as right neighbor) according to the initial configuration in Eq.
(69). Correspondingly, the term $D_{1}E$ in the averages $\langle
D_{1}E\rangle^{(1)}$ and $\langle D_{1}E\rangle^{(2)}$ would come from the
product sequences $(\tau D_{1}+(1-\tau)D_{2})E^{m}A$ and $D_{1}E^{s}$
respectively Eq. (70). The expression for $J_{10}^{(1)}$ is given by

$\displaystyle J_{10}^{(1)}=p_{1}\langle
D_{1}E\rangle^{(1)}=\frac{\rho_{+}}{Q_{N_{+}}}\sum_{n_{1}=0}^{\infty}..\sum_{n_{N_{+}}=0}^{\infty}\sum_{m_{1}=1}^{\infty}..\sum_{m_{N_{+}}=0}^{\infty}\sum_{r_{1}=0}^{\infty}..\sum_{r_{\bar{N}}=0}^{\infty}..\sum_{s_{\bar{N}}=0}^{\infty}$
$\displaystyle\mathrm{Tr}\left[z_{0}p_{1}D_{1}E(z_{0}E)^{m_{1}-1}A(z_{0}E)^{n_{1}}\prod_{k=2}^{N_{+}}(D_{1}+D_{2})(z_{0}E)^{m_{k}}A(z_{0}E)^{n_{k}}\prod_{k=1}^{\bar{N}}D_{2}(z_{0}E)^{r_{k}}\prod_{k=1}^{\bar{N}}D_{1}(z_{0}E)^{s_{k}}\right].$
(91)

The construction of Eq. (91) follows similar arguments as of Eq. (84), 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_{10}^{(2)}$ is

$\displaystyle J_{10}^{(2)}=p_{1}\langle
D_{1}E\rangle^{(2)}=\frac{1}{Q_{N_{+}}}\sum_{n_{1}=0}^{\infty}..\sum_{r_{\bar{N}}=0}^{\infty}\sum_{s_{1}=1}^{\infty}..\sum_{s_{\bar{N}}=0}^{\infty}$
(92)
$\displaystyle\mathrm{Tr}\left[\prod_{k=1}^{N_{+}}(D_{1}+D_{2})(z_{0}E)^{m_{k}}A\prod_{k=1}^{\bar{N}}D_{2}(z_{0}E)^{r_{k}}z_{0}p_{1}D_{1}E(z_{0}E)^{s_{1}-1}\prod_{k=2}^{\bar{N}}D_{1}(z_{0}E)^{s_{k}}\right].$
(93)

Obviously in Eq. (93), the lower limit of the index $s_{1}$ is shifted to $1$
from $0$. Using the matrix algebra and matrix representations from Eqs. (15)
and (29), we get from Eqs. (91) and (93):

$\displaystyle
J_{10}^{(1)}=z_{0}\rho_{+}\,\frac{w_{21}}{\left(1-\frac{z_{0}}{p_{1}}\right)}\,\frac{1}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]},$
$\displaystyle J_{10}^{(2)}=\frac{z_{0}}{2}(1-2\rho_{+}-\rho_{0}).$ (94)

Substituting Eq. (94) into Eq. (89), 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:

$\displaystyle J_{10}$ $\displaystyle=$ $\displaystyle
z_{0}\rho_{+}\,\frac{w_{21}}{\left(1-\frac{z_{0}}{p_{1}}\right)}\,\frac{1}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]}\,+\,\frac{z_{0}}{2}(1-2\rho_{+}-\rho_{0}),$
(95) $\displaystyle J_{20}$ $\displaystyle=$ $\displaystyle
z_{0}\rho_{+}\,\frac{w_{12}}{\left(1-\frac{z_{0}}{p_{2}}\right)}\,\frac{1}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]}\,+\,\frac{z_{0}}{2}(1-2\rho_{+}-\rho_{0}),$
(96) $\displaystyle J_{+0}$ $\displaystyle=$ $\displaystyle\rho_{+}z_{0}.$
(97)

In Figs. 4 and 5 we present the variation of the drift currents as functions
of vacancy density $\rho_{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 $\rho_{0}$,
the chances for hopping increase, thereby increasing the drift current in Fig.
4. However, after a particular value of $\rho_{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
$\rho_{0}$ (Fig. 4). The maximum of the drift current for different species
generally occur at distinct values of $\rho_{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 $\mu,$ the general expression for the drift
current of any non-conserved species $I$ ($I=1,\dots,\mu$) is obtained to be

$\displaystyle
J_{I0}=z_{0}\rho_{+}\,\frac{d_{I}}{\left(1-\frac{z_{0}}{p_{I}}\right)}\,\frac{1}{\sum\limits_{K=1}^{\mu}\frac{d_{K}}{\left(1-\frac{z_{0}}{p_{K}}\right)}}\,+\,\frac{z_{0}}{\mu}(1-2\rho_{+}-\rho_{0}),$
(98)

where the $d_{I}$ is the solution of Eq. (66).

Figure 6: The figure illustrates the variation of two-point correlation
$C_{01}$ between the vacancy and species $1$ by tuning the vacancy density
$\rho_{0}$. We observe that this correlation initially starts increasing with
increasing $\rho_{0}$, reaches to a local maximum, then decreases. However, as
$\rho_{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,\epsilon=0.1,w_{12}=0.8,w_{21}=1.0,\rho_{+}=0.2$.
The ensemble average is done over $10^{5}$ samples. Figure 7: The figure shows
a crossover from negative correlation to positive correlation for $C_{02}$
(between vacancy and species $2$) with increasing vacancy density $\rho_{0}$.
There is some special intermediate density for which $C_{02}$ becomes zero.
The parameters used are
$L=10^{3},p_{1}=0.3,p_{2}=1.0,\epsilon=0.1,w_{12}=0.8,w_{21}=1.0,\rho_{+}=0.2$.
The ensemble average is done over $10^{5}$ samples.

###  4.3 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 $\langle 10\rangle$ and $\langle 20\rangle$ 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 $\langle
01\rangle$ and $\langle 02\rangle$. We find the exact expressions for the
corresponding two-point correlations under the periodic boundary condition as

$\displaystyle C_{01}=\langle 0\,1\rangle-\langle 0\rangle\langle 1\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{z_{0}}{\epsilon}-\rho_{0}\right)\,\frac{w_{21}}{\left(1-\frac{z_{0}}{p_{1}}\right)}\,\frac{\rho_{+}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]}$
(100)
$\displaystyle+\left(\frac{z_{0}}{2p_{1}}-\frac{\rho_{0}}{2}\right)(1-2\rho_{+}-\rho_{0}),$
$\displaystyle C_{02}=\langle 0\,2\rangle-\langle 0\rangle\langle 2\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{z_{0}}{\epsilon}-\rho_{0}\right)\,\frac{w_{12}}{\left(1-\frac{z_{0}}{p_{2}}\right)}\,\frac{\rho_{+}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}}\right]}$
(102)
$\displaystyle+\left(\frac{z_{0}}{2p_{2}}-\frac{\rho_{0}}{2}\right)(1-2\rho_{+}-\rho_{0}).$

The correlation $C_{01}$ is plotted against the vacancy density $\rho_{0}$ in
Fig. 6. As $\rho_{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 $\rho_{0}$. In other words, $C_{01}$ increases
with increasing vacancy density for both sufficiently high and sufficiently
low values of $\rho_{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 $\rho_{0}$ for
which the special arrangements of the accessible steady state configurations
makes the average correlation $C_{02}$ to be zero.

(a) (b)

Figure 8: The figures exhibit the variation of the flip current
$J_{I\leftrightarrow K}$ with vacancy density $\rho_{0}$ and flip rate
$w_{12}$ in (a) and (b), respectively. The analytical and Monte Carlo
simulation results show good agreement. The flip current monotonically
decreases with increasing $\rho_{0}$ whereas it increases monotonically with
increasing $w_{12}$. The common parameters for both figures (a) and (b) are
$L=10^{3},p_{1}=0.3,p_{2}=1.0,p_{3}=1.0,\epsilon=0.1,w_{21}=0.5,w_{23}=0.5,w_{32}=0.2,w_{31}=0.8,w_{13}=0.2,\rho_{+}=0.15$.
For (a), $w_{12}=0.4$ and for (b), $\rho_{0}=0.22$. The ensemble average is
done over $10^{5}$ samples.

(a) (b)

(c) (d)

Figure 9:  The set of figures exhibits the dependence of steady state values
of observables on the choice of initial configuration. We observe clear
deviations between SIC (denoted by circles) and RIC (denoted by rectangles) in
(a)-(d) and the amount of deviation in each figure changes with the variation
of input parameter. The common set of parameters used for (a)-(d) are
$L=10^{3},p_{1}=0.3,\epsilon=0.1,\rho_{+}=0.216,\rho_{0}=0.203,\rho_{1}(0)=0.321,\rho_{2}(0)=0.260$.
The other parameters for (a), (b) are $p_{2}=1.0$, $w_{21}=0.6$ and for (c),
(d) are $w_{12}=0.5,w_{21}=0.4$. The ensemble averages are done over $10^{5}$
samples.

###  4.4 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 $\mu=2$. Therefore, in order to have a net non-zero
flip current between pairs of species, here we consider the case $\mu=3$
($I=1,2,3$). We denote the net flip current between species $I$ and $K$ as
$J_{I\leftrightarrow K}$. For $\mu=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

$\displaystyle J_{I\leftrightarrow K}=w_{IK}\langle I+\rangle-w_{KI}\langle
K+\rangle=\rho_{+}\,\frac{(w_{12}w_{23}w_{31}-w_{21}w_{13}w_{32})}{\left[\frac{d_{1}}{1-\frac{z_{0}}{p_{1}}}+\frac{d_{2}}{1-\frac{z_{0}}{p_{2}}}+\frac{d_{3}}{1-\frac{z_{0}}{p_{3}}}\right]}.$
(103)

The initial configuration that we have used to arrive at Eq. (104) is the one
in Eq. (76). In Eq. (104), 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 $\rho_{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 flip-current decreases monotonically in a nonlinear
manner with increasing vacancy density. The reason behind this, is the
decrease in species densities with increasing $\rho_{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. (104) for any $\mu$, under the periodic
boundary condition, is obtained as

$\displaystyle J_{I\leftrightarrow
K}=\rho_{+}\,\frac{(d_{I}w_{IK}-d_{K}w_{KI})}{\sum\limits_{K=1}^{\mu}\frac{d_{K}}{\left(1-\frac{z_{0}}{p_{K}}\right)}},$
(104)

where $d_{I},d_{K}$ are the solutions of Eq. (66). We should mention that, for
$\mu>3$, the flip currents between different pairs of species would be in
general distinct from one another.

###  4.5 Non-ergodicity: dependence on initial configuration

Here we establish the non-ergodicity of the $\mu$-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 $\mu=2$. We choose two different initial configurations as
follows. (i) The initial configuration given in Eq. (69) 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},\epsilon,w_{12},w_{21}$) and the densities
($\rho_{+},\rho_{0},\rho_{1}(0),\rho_{2}(0)$), where $\rho_{1}(0),\rho_{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},\epsilon,w_{12},w_{21},\rho_{+},\rho_{0}$) and same initial
densities of the species ($\rho_{1}(0),\rho_{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
$\rho_{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}=\langle 11\rangle-\rho_{1}^{2}$ 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}=\langle 22\rangle-\rho_{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
$\left\\{+s_{i}s_{i+1}\dots s_{n}s_{n+1}+\right\\}$ in an initial
configuration, where $s_{j}=0,1,\dots\mu$ but $s_{j}\neq+$ for $i\leqslant
j\leqslant(n+1)$, then the ordering of different species $1,\dots,\mu$ (not
vacancies) for $i\leqslant j\leqslant 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. (4)] and the role of corresponding underlying
symmetries. An explicit illustration of the block-diagonal structure of the
rate matrix in the $\mu$-ASEP-IAF, for small system sizes, is presented in
Appendix E.

##  5 Partially asymmetric generalization: $\mu$-ASEP-IAF

In this section, we consider the $\mu$-ASEP-IAF under periodic boundary
conditions, a generalization of the $\mu$-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,\dots,\mu$), 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,

$\displaystyle\mathrm{drift\,(species):}\hskip 14.22636ptI0\,\,$
$\displaystyle\xrightleftharpoons[q_{I}]{p_{I}}$ $\displaystyle\,\,0I\hskip
14.22636ptI=1,2,...,\mu$ $\displaystyle\mathrm{drift\,(impurity):}\hskip
8.5359pt+0\,\,$
$\displaystyle\stackrel{{\scriptstyle\epsilon}}{{\longrightarrow}}$
$\displaystyle\,\,0+$ $\displaystyle\mathrm{flip:}\hskip 11.38092ptI+\,\,$
$\displaystyle\xrightleftharpoons[w_{KI}]{w_{IK}}$ $\displaystyle\,\,K+\hskip
4.26773ptI,K=1,...,\mu.$ (105)

The $\mu$-ASEP-IAF remains non-ergodic in nature. Since we could obtain the
steady state of the $\mu$-TASEP-IAF using matrix product ansatz [Eq. (3)], we
assume the same can be done for the partially asymmetric motion also.

###  5.1 Matrix algebra, auxiliaries and matrix representations

The matrix algebra for the dynamics in Eq. (105) under the periodic boundary
condition, is

$\displaystyle p_{K}D_{K}E-q_{K}ED_{K}$ $\displaystyle=$ $\displaystyle
D_{K},\hskip 65.44142ptK=1,\dots,\mu$ (106) $\displaystyle\epsilon AE$
$\displaystyle=$ $\displaystyle A,$ (107)
$\displaystyle\sum_{\begin{subarray}{c}I=1\\\ I\neq
K\end{subarray}}^{\mu}w_{IK}D_{I}A$ $\displaystyle=$ $\displaystyle
D_{K}A\sum_{\begin{subarray}{c}I=1\\\ I\neq
K\end{subarray}}^{\mu}w_{KI},\hskip 14.22636ptK=1,\dots,\mu.$ (108)

In comparison to the matrix algebra [Eq. (12)] for the $\mu$-TASEP-IAF, the
only changes occurring in Eq. (108) 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}ED_{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 $\mu$-ASEP-IAF in Eq.
(108), can be considered as a generalization of the matrix algebra in Ref.
[55]. To arrive at the matrix algebra in Eq. (108) from the dynamics (105),
ansatz (3) and flux cancellation condition (6), the choice of the auxiliary
matrices are the same as the totally asymmetric case, i.e.

$\displaystyle\tilde{E}=1,\,\,\tilde{A}=0,\,\,\tilde{D}_{K}=0\hskip
14.22636ptK=1,2,\dots,\mu.$ (109)

However, unlike the totally asymmetric case, we find the matrix
representations for the $\mu$-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 $\mu=3$ case explicitly (as we will stick to
$\mu=3$ for the discussion of observable in this section) and mention the
changes required to construct the matrices for any $\mu>0$.

$\mu=3:$ A possible set of representations of the matrices for the $3$-ASEP-
IAF $(I=1,2,3)$ is the following

$\displaystyle E$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}0&0&0&0&.&.\\\ 1&0&0&0&.&.\\\
0&1&0&0&.&.\\\ 0&0&1&0&.&.\\\ 0&0&0&1&&\\\ .&.&&&.&\\\ .&.&&&&.\\\
\end{array}\right),\hskip
9.95863ptA=\left(\begin{array}[]{cccccc}1&\frac{1}{\epsilon}&\frac{1}{\epsilon^{2}}&\frac{1}{\epsilon^{3}}&.&.\\\
0&0&0&0&.&.\\\ 0&0&0&0&.&.\\\ .&.&.&.&.&.\\\ .&.&.&.&.&.\\\
\end{array}\right)$ (122) $\displaystyle D_{I}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}d_{I}^{1,1}&d_{I}^{1,2}&d_{I}^{1,3}&d_{I}^{1,4}&.&.\\\
0&d_{I}^{2,2}&d_{I}^{2,3}&d_{I}^{2,4}&.&.\\\
0&0&d_{I}^{3,3}&d_{I}^{3,4}&.&.\\\ 0&0&0&d_{I}^{4,4}&.&.\\\ .&.&&&.&\\\
.&.&&&&.\\\ \end{array}\right),\hskip 8.5359ptI=1,2,3$ (133) $\displaystyle
d_{I}^{m,m+r}=\frac{(m)_{r}}{r!\,p_{I}^{r}}\left(\frac{q_{I}}{p_{I}}\right)^{m-1}d_{I}^{1,1},\hskip
22.76228pt\forall r\geqslant 0$ $\displaystyle
d_{1}^{1,1}=w_{21}w_{31}+w_{23}w_{31}+w_{32}w_{21},$ $\displaystyle
d_{2}^{1,1}=w_{12}w_{32}+w_{13}w_{32}+w_{31}w_{12},$ $\displaystyle
d_{3}^{1,1}=w_{13}w_{23}+w_{12}w_{23}+w_{21}w_{13}.$

In the absence of the impurities ($A$) and the flip processes, the term
$d_{\mu}^{1,1}$ for every $\mu$ becomes unity, and corresponding the matrix
representations for $D_{I}$ and $E$ in Eq. (133) (and their generalizations
for general $\mu$) are the same as that of the conserved disordered ASEP [55,
27]. In Eq. (133), we observe that the matrices $(D_{I})$ corresponding to the
species $(I)$ are upper triangular. The subscript $I$ in matrix element
$d_{I}^{i,j}$ 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_{I}^{m,m+r}$ corresponds to the Pochhammer symbol for
rising factorials, $(m)_{r}=m(m+1)(m+2)\dots(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 $\mu=2$ the only changes in comparison to Eq. (133) will be in the values
of $d_{I}^{1,1}$ ($I=1,2$), which would be simply $d_{1}^{1,1}=w_{21}$ and
$d_{2}^{1,1}=w_{12}$. In fact, the matrix representations [Eq. (133)] for the
$3$-ASEP-IAF can be generalized for any $\mu>0$ in a straightforward manner.
The representations will remain the same, only the values of $d_{I}^{1,1}$
($I=1,2,\dots,\mu$) would change where $d_{I}^{1,1}$ is the solution of Eq.
(66).

###  5.2 Partition function for special initial configuration

In the totally asymmetric case, for general $\mu$, we have considered the
specific initial configuration Eq. (80) which leads us to acquire analytical
expressions for observables of interest. We choose a particular case of Eq.
(80), namely the $\bar{\rho}=0$ case, as the special initial configuration for
analytical calculation in $\mu$-ASEP-IAF. More precisely, the choice of our
special initial configuration for $\mu$-ASEP-IAF is,

$\displaystyle
C(0)\equiv\prod_{i=1}^{N_{+}/\mu}D_{1}A\prod_{i=1}^{N_{+}/\mu}D_{2}A\dots\prod_{i=1}^{N_{+}/\mu}D_{\mu}A\prod_{i=1}^{N_{0}}E.$
(134)

The initial configuration Eq. (134) is chosen in a way that fixes the total
impurity density to $\rho_{+}$. In comparison to Eq. (80), the initial
configuration in Eq. (134) 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 $\mu$-ASEP-IAF. We find the partition
function in the steady state under the periodic boundary condition
corresponding to the initial configuration Eq. (134), to be

$\displaystyle
Q_{N_{+}}(z_{0})=\left(\left[\sum_{I=1}^{\mu}\frac{d_{I}}{1-\frac{z_{0}}{p_{I}}-\frac{q_{I}}{p_{I}}\frac{z_{0}}{\epsilon}}\right]\left(\frac{1}{1-\frac{z_{0}}{\epsilon}}\right)\right)^{N_{+}}.$
(135)

We have used short hand notations $d_{I}\equiv d_{I}^{1,1}$ which are
essentially the solutions of Eq. (66). When $q_{I}=0$ for all species, it is
straightforward to check that the partition function in Eq. (135) reduces to
the partition function Eq. (82) of the totally asymmetric case, under the
condition $\bar{\rho}=0$.

###  5.3 Species densities, drift current and flip current

Just like we did in the $\mu$-TASEP-IAF, we can analytically calculate several
observables of interest in the partially asymmetric case also, using the
matrix algebra (108) and matrix representations (133) following the same
procedures as before. Starting from the initial configuration stated in Eq.
(134), the average density $\rho_{I}$ of any non-conserved species $I$
($I=1,\dots,\mu$) for the $\mu$-ASEP-IAF is obtained as

$\displaystyle\rho_{I}=\rho_{+}\,\frac{d_{I}}{\left(1-\frac{z_{0}}{p_{I}}-\frac{q_{I}}{p_{I}}\frac{z_{0}}{\epsilon}\right)}\,\frac{1}{\sum\limits_{K=1}^{\mu}\frac{d_{K}}{\left(1-\frac{z_{0}}{p_{K}}-\frac{q_{K}}{p_{K}}\frac{z_{0}}{\epsilon}\right)}},$
(136)

where $\rho_{0}$ and $\rho_{+}$ are the conserved densities for the vacancies
and the impurities, respectively. If we put $q_{I}=0$ for all $I$ in Eq.
(136), the expression of $\rho_{I}$ for the totally asymmetric case [Eq. (88)]
is correctly recovered. The drift currents $J_{I0}$ and the flip currents
$J_{I\leftrightarrow K}$ for the $\mu$-ASEP-IAF ($I,K=1,\dots,\mu$) are given
by

$\displaystyle J_{I0}$ $\displaystyle=$ $\displaystyle
z_{0}\,\rho_{+}\,\frac{d_{I}}{\left(1-\frac{z_{0}}{p_{I}}-\frac{q_{I}}{p_{I}}\frac{z_{0}}{\epsilon}\right)}\,\frac{1}{\sum\limits_{K=1}^{\mu}\frac{d_{K}}{\left(1-\frac{z_{0}}{p_{K}}-\frac{q_{K}}{p_{K}}\frac{z_{0}}{\epsilon}\right)}},$
(137) $\displaystyle J_{I\leftrightarrow K}$ $\displaystyle=$
$\displaystyle\rho_{+}\,\frac{(d_{I}w_{IK}-d_{K}w_{KI})}{\sum\limits_{K=1}^{\mu}\frac{d_{K}}{\left(1-\frac{z_{0}}{p_{K}}-\frac{q_{K}}{p_{K}}\frac{z_{0}}{\epsilon}\right)}}.$
(138)

###  5.4 Negative differential mobility

In what follows, we will show that the species in the $\mu$-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 [82, 79, 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 $\mu$-ASEP-IAF can exhibit NDM for
particular choices of the rates in the microscopic dynamics.

(a) (b)

(c)

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
$\epsilon$) 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},\epsilon=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,\rho_{+}=0.3,\rho_{0}=0.4$.
The ensemble average is done over $10^{6}$ samples.

To illustrate NDM in $\mu$-ASEP-IAF, we will focus on the $\mu=3$ case. We
have three species of particles $(I=1,2,3)$, impurities $(+)$ and vacancies in
the system following the microscopic dynamics Eq. (105). We choose the drift
rates $p_{I}$ and $q_{I}$ of the species $I=1,2,3$ to be

$\displaystyle\begin{array}[]{ccc}p_{1}=1,q_{1}=e^{-b},&p_{2}=\frac{1}{1+b^{2}}=q_{2},&p_{3}=1=q_{3}.\end{array}$
(140)

The choices of the hopping rates in Eq. (140) are inspired by similar choices
in Ref. [85] in context of NDM for different models. The special choices of
the hopping rates in Eq. (140) allow us to identify the parameter $b$ as the
hopping bias in the system. This is because $\mathrm{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. (140) 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 $\epsilon>0$ where the impurity acts as the lone
driven mode. In the present context, we choose $\epsilon$ 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 $\epsilon$, 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_{IK}$ ($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. (138) 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. (140)] 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 $\mu$-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 [86, 52, 48]. This urges
for detailed investigation of the counter flow situation in $\mu$-ASEP-IAF in
future works.

##  6 Summary and future directions

In this article, we have obtained an exact steady state probability
distribution of the $\mu$-ASEP-IAF on a one dimensional lattice under periodic
boundary conditions, using the matrix product ansatz. The $\mu$-ASEP-IAF
consists of (i) drift of the species ($I=1,2,\dots,\mu$) and impurities, and
(ii) flip between different species initiated by the impurities. In steady
state, we provide the explicit finite dimensional [$(\mu+1)\times(\mu+1)$]
matrix representations for any $\mu>0$ for the totally asymmetric case i.e.
$\mu$-TASEP-IAF. For the partially asymmetric scenario i.e. $\mu$-ASEP-IAF, we
obtain the corresponding matrices with infinite dimensional representations.
Importantly, due to the non-ergodicity of the $\mu$-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
$\mu>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 $\mu$-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 $\mu$-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
$\mu$-ASEP-IAF studied here is relevant in two other important contexts. The
$\mu$-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 $\mu$-TASEP-IAF, these connections are discussed in
details in Appendix A and Appendix C respectively. Importantly, the exact
solution of $\mu$-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 $\mu$-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 $\mu$-ASEP-IAF model using matrix product ansatz. Just like the
multi lane traffic flow, the connections between $\mu$-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 $\mu$-ASEP-IAF. Another important future direction would be
to analyze the exact steady state of the $\mu$-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 $\mu$-ASEP-IAF. We would also like to
investigate the dynamics of the $\mu$-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.

## Acknowledgements

We thank Kazuaki Takasan for fruitful discussions. We gratefully acknowledge
Arvind Ayyer for pointing out important references. This work is partially
supported by the Grants-in-Aid for Scientific Research (Grant No. 21H01006).

##  Appendix A Connection between $\mu$-TASEP-IAF and multi lane TASEP

Here we explore the connections between the multi lane totally asymmetric
simple exclusion process ($m$-TASEP) and the one dimensional $\mu$-ASEP-IAF.
For simplicity, we consider $\mu=2$ i.e. the $2$-TASEP-IAF (and
correspondingly $2$-TASEP or two lane TASEP).

In Fig. 11, we present a two lane TASEP where particles can hop in forward
directions in lane $1$ and lane $2$ with rates $p_{1}$ and $p_{2}$
respectively. The particles in lane $1$ and lane $2$ can be interpreted as two
types i.e. species $1$ and species $2$ particles in the $2$-TASEP-IAF process.

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.

Except for the hopping of particles in the two lanes in the $2$-TASEP, the
particles can change lanes through bridges connecting the lanes. There are two
types of bridges, active (green vertical lines in Fig. 11) that allows
vertical hopping i.e. lane change of particles and inactive (red vertical
lines in Fig. 11) that does not allow lane change of particles. The active
bridges in the two lane TASEP mimic the impurities ($+$) in the $2$-TASEP-IAF.
However, the inactive bridges are not counted in the equivalent $2$-TASEP-IAF.
Notably, a neighboring pair of (active, inactive) bridges can change to
(inactive, active). This inactive-active transformation of neighboring lanes
can be interpreted as a resultant drift of the active bridges through the
system. Consequently, this accounts for the forward hopping of impurity in the
$2$-TASEP-IAF.

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.

The equivalence of the microscopic dynamics of the two lane TASEP and the
$2$-TASEP-IAF is shown in Fig. 12. The last (bottom) panel in Fig. 12 exhibits
the connection between the inactive-active lane transformations in $2$-TASEP
and the impurity hopping in the $2$-TASEP-IAF. The two panels above the bottom
panel in Fig. 12 describe the equivalence of the lane change of particles in
$2$-TASEP with the impurity activated flip in the corresponding $2$-TASEP-IAF.
To elaborate, when the lane $1$ ($2$) particle in $2$-TASEP comes in contact
with an active bridge, it can go to lane $2$ ($1$) with rate $w_{12}$
($w_{21}$). Similarly, when a species $1$ ($2$) particle in the $2$-TASEP-IAF
encounters an impurity as a right neighbor, it can flip to a species $2$ ($1$)
particle with rate $w_{12}$ ($w_{21}$). The first two panels in Fig. 12
present the relations between usual hopping dynamics in the two processes.

We can generalize the approaches described above to establish connections
between the multi-lane TASEP and the $\mu$-TASEP-IAF for any $\mu\geqslant 3$.
It is noteworthy that, for the multi-species case $(\mu\geqslant 3)$, we have
shown the existence of non-zero net flip current in the $\mu$-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 $\mu$-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
$\mu$-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 $\mu$-TASEP-IAF.

##  Appendix B A variation of $\mu$-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
$\mu$-TASEP-IAF and the multi-lane TASEP discussed in Appendix A, it is
natural to ask about the applicability of $\mu$-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 $\mu$-TASEP-IAF with $\mu=2,$ we plot the total drift current of
species $1$ and species $2$, $J_{\mathrm{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}$.

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}\in(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,\epsilon=0.1,w_{21}=1.0,\rho_{+}=0.2,\rho_{0}=0.2$.
The ensemble average is done over $10^{7}$ samples.

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}\in(0,1)$. The reason behind
this, as revealed by a careful observation, is the approximate mapping between
the $\mu$-TASEP-IAF and multi lane TASEP described in Fig. 12.

Figure 14: The figure illustrates the comparison of the total drift current
$J_{\mathrm{total}}$ between the flip process in Eq. (141) and the proposed
variation in Eq. (142). The current $J_{\mathrm{total}}$ for the variation Eq.
(142) increases considerably with increasing flip rate (equivalent lane change
rate in two lane TASEP) $w_{12}$, whereas the same for the original flip
process Eq. (141) increase much slowly (Fig. 13). The parameters used are
$L=10^{4},p_{1}=0.3,p_{2}=1.0,\epsilon=0.1,w_{21}=1.0,\rho_{+}=0.2,\rho_{0}=0.2$.
The ensemble average is done over $10^{7}$ samples.

In the $\mu$-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
$\mu$-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

$\displaystyle I+\,\,$ $\displaystyle\stackrel{{\scriptstyle
w_{IK}}}{{\longrightarrow}}$ $\displaystyle\,\,K+,$ (141)

with $I,K=1,...,\mu$ and $I\neq K$. Whereas we propose the new flip dynamics
to be

$\displaystyle I+\,\,$ $\displaystyle\stackrel{{\scriptstyle
w_{IK}}}{{\longrightarrow}}$ $\displaystyle\,\,+K.$ (142)

Note that, in comparison to Eq. (141), the flip process in Eq. (142)
accompanies the flip of the species with a hop towards right. The other
hopping processes in Eq. (1) remain the same for this varied $\mu$-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. (141) and Eq. (142)
(using Monte Carlo simulations). Indeed, the Fig. 14 shows that
$J_{\mathrm{total}}$ increases considerably with increasing lane change rate
for Eq. (142) whereas it grows weakly for Eq. (141) (see Fig. 13). This
observation implies that the proposed variation of the $\mu$-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 $\mu$-TASEP-IAF
both analytically and numerically and to build connections with the multi lane
traffic flow.

##  Appendix C Connection between $\mu$-TASEP-IAF and enzymatic chemical
reactions

In this appendix, we briefly discuss some connections between the $\mu$-TASEP-
IAF and enzymatic chemical reactions. One of the simplest form of the
enzymatic chemical reaction is,

$E+S\xrightleftharpoons[k_{b}]{k_{f}}ES\,\xrightleftharpoons[\bar{k}_{b}]{\bar{k}_{f}}\,E+P,$
(143)

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 $\bar{k}_{f}$ and $\bar{k}_{b}$ 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 $\mu$-ASEP-IAF, we would rather consider a much
simplified version of the chemical reaction (143) as

$S+E\,\xrightleftharpoons[k_{P}]{k_{S}}\,P+E,$ (144)

where we ignore the intermediate complex formation. $k_{S}$ ($k_{P}$) is the
rate constant for $S$ transforming to $P$ ($P$ transforming to $S$).

Figure 15: The figure illustrates the connection between the simplified form
of the enzymatic chemical reaction Eq. (144) (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.

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 $\mu$-TASEP-IAF. As explained in Fig. 15 for $\mu=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. (144), 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 ($\mu>2$) case of
the $\mu$-TASEP-IAF can be mapped to generalized version of the chemical
reaction Eq. (144) as

$S_{i}+E\,\xrightleftharpoons[k_{P_{j}}]{k_{S_{j}}}\,P_{j}+E,$ (145)

where $i=1,\dots,\mu_{1}$ are identified as $\mu_{1}$ number of substrates and
$j=1,\dots,\mu_{2}$ are identified as $\mu_{2}$ number of products in the
enzymatic chemical reaction system (with $E$ acting as the enzyme for each
chemical reaction) where $\mu_{1}+\mu_{2}=\mu$ ($\mu$ being the total number
of species in the $\mu$-TASEP-IAF). To study the time evolution of the
enzymatic chemical reaction, one has to study the dynamics of the $\mu$-TASEP-
IAF. The connection between the $\mu$-TASEP-IAF and the enzymatic chemical
reactions should be studied more thoroughly with proper attention to the
observables of interest for the chemical reactions.

##  Appendix D Solution of Eq. (75) 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},\epsilon,w_{12},w_{21},\rho_{0},\rho_{+})$
only, an important step in the calculation is to solve Eq. (75) to obtain the
fugacity as a function of these input parameters i.e.
$z_{0}(p_{1},p_{2},\epsilon,w_{12},w_{21},\rho_{0},\rho_{+})$. However, in
most cases the solutions of $z_{0}$ from Eq. (75), 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
$\bar{\rho}=0$ for the initial configuration (see Sec. 3). Consequently, the
density-fugacity relation Eq. (75) becomes a quartic equation in the variable
$z_{0}$, emerging from

$\displaystyle\rho_{0}=\rho_{+}z_{0}\left[\frac{1}{\epsilon-
z_{0}}+\frac{1}{p_{1}w_{21}(p_{2}-z_{0})+p_{2}w_{12}(p_{1}-z_{0})}\left(p_{1}w_{21}\frac{p_{2}-z_{0}}{p_{1}-z_{0}}+p_{2}w_{12}\frac{p_{1}-z_{0}}{p_{2}-z_{0}}\right)\right].$
(146)

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\neq\epsilon$. As a result, Eq. (146) is reduced to a quadratic
equation which leads to the following solution

$z_{0}=\frac{(1+\epsilon)(\rho_{0}+\rho_{+})-\sqrt{(1+\epsilon)^{2}(\rho_{0}+\rho_{+})^{2}-4\epsilon\rho_{0}(\rho_{0}+2\rho_{+})}}{2(\rho_{0}+2\rho_{+})}.$
(147)

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}=\epsilon=w_{21}=1,\rho_{+}=1/4$. Here the fugacity would be a
function of $w_{12}$ and $\rho_{0}$ i.e. $z_{0}(w_{12},\rho_{0})$. The reason
behind keeping $w_{12}$ and $\rho_{0}$ as free parameters is that, the
observables in the main text have been mostly analyzed as functions of these
two parameters. The corresponding solution (Eq. (146) essentially reduces to a
cubic equation) for the fugacity is given below,

$z_{0}=-\frac{a_{2}}{3a_{3}}+\frac{2^{1/3}(-a_{2}^{2}+3a_{3}a_{1})}{3a_{3}|\nu|}-\frac{2^{2/3}|\nu|}{6a_{3}},$
(148)

where $|\nu|$ denotes the absolute value of $\nu$ and its functional form is

$\displaystyle\nu(a_{0},a_{1},a_{2},a_{3})=$ (149)
$\displaystyle(-2a_{2}^{3}+9a_{3}a_{2}a_{1}-27a_{3}^{2}a0+\sqrt{-4(a_{2}^{2}-3a_{3}a_{1})^{3}+(2a_{2}^{3}-9a_{3}a_{2}a_{1}+27a_{3}^{2}a_{0})^{2}})^{\frac{1}{3}}.\hskip
14.22636pt$ (150)

In Eqs. (148) and (150), the parameters $a_{0},a_{1},a_{2},a_{3}$ are explicit
functions of $w_{12}$ and $\rho_{0}$, as follows

$\displaystyle a_{0}$ $\displaystyle=$ $\displaystyle-4\rho_{0}(1+w_{12}),$
(151) $\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle
3+2w_{12}+4\rho_{0}(4+5w_{12}),$ (152) $\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle-7-8w_{12}-4\rho_{0}(5+8w_{12}),$ (153) $\displaystyle a_{3}$
$\displaystyle=$ $\displaystyle 4(1+2\rho_{0})(1+2w_{12}).$ (154)

Although there seems to be no definite rule for obtaining closed-form
solutions of $z_{0}$ like Eqs. (147) and (148), one might achieve other
convenient solutions by searching for suitable subspace of the transition
rates.

##  Appendix 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. (4)], reflecting the non-ergodicity of $\mu$-ASEP-IAF. To
illustrate this with an example for $\mu=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$)

$\begin{array}[]{ccc}11+0\equiv I_{1},&110+\equiv I_{2},&101+\equiv I_{3},\\\
12+0\equiv II_{1},&120+\equiv II_{2},&102+\equiv II_{3},\\\ 21+0\equiv
III_{1},&210+\equiv III_{2},&201+\equiv III_{3},\\\ 22+0\equiv
IV_{1},&220+\equiv IV_{2},&202+\equiv IV_{3}.\\\ \end{array}\\\ $ (155)

We have divided the 12 configurations in Eq. (155) into 4 sectors
$I,II,III,IV$ 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 $IV_{3}$),

$M=\left(\begin{array}[]{cc}M_{I,II}&\bigcirc\\\ \bigcirc&M_{III,IV}\\\
\end{array}\right),$ (156)

where

$\displaystyle M_{I,II}=\left(\begin{array}[]{cccccc}-\epsilon-
w_{12}&0&p_{1}&w_{21}&0&0\\\ \epsilon&-p_{1}&0&0&0&0\\\
0&p_{1}&-p_{1}-w_{12}&0&0&w_{21}\\\ w_{12}&0&0&-\epsilon-w_{21}&0&p_{1}\\\
0&0&0&\epsilon&-p_{2}&0\\\ 0&0&w_{12}&0&p_{2}&-p_{1}-w_{21}\\\
\end{array}\right),$ (163) $\displaystyle
M_{III,IV}=\left(\begin{array}[]{cccccc}-\epsilon-w_{12}&0&p_{2}&w_{21}&0&0\\\
\epsilon&-p_{1}&0&0&0&0\\\ 0&p_{1}&-p_{2}-w_{12}&0&0&w_{21}\\\
w_{12}&0&0&-\epsilon-w_{21}&0&p_{2}\\\ 0&0&0&\epsilon&-p_{2}&0\\\
0&0&w_{12}&0&p_{2}&-p_{2}-w_{21}\\\ \end{array}\right),$ (170)

and $\bigcirc$ is $6\times 6$ null matrix. In Eqs. (156) and (170), 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 $II$ through flip
dynamics, whereas sector $III$ and $IV$ are also connected to each other via
flip dynamics. However, sectors $(I,II)$ are disconnected from sectors
$(III,IV)$, thereby creating two separate blocks in the rate matrix.

Note that, in absence of the flip dynamics (i.e. $w_{12}=w_{21}=0$), sectors
$I$ become disconnected from $II$, similarly $III$ gets disconnected from
$IV$, 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 $N_{1}+N_{2}=L-N_{0}-N_{+}$,

$L$ | $N_{+}$ | $N_{blocks}$
---|---|---
4 | 1 | 2
4 | 2 | 1
5 | 1 | 4
5 | 2 | 3
5 | 3 | 1
6 | 1 | 8
6 | 2 | 6
6 | 3 | 3
6 | 4 | 1
7 | 1 | 16
7 | 2 | 15
8 | 1 | 32
8 | 2 | 32
9 | 1 | 64
9 | 2 | 74
10 | 1 | 128
10 | 2 | 160

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

$\displaystyle N_{+}=1:$ $\displaystyle\hskip 28.45274ptN_{blocks}=2^{L-3},$
(171) $\displaystyle N_{+}=2:$ $\displaystyle\hskip
28.45274ptN_{blocks}=2^{L-6}L,\hskip 8.5359ptL\hskip 8.5359pt\mathrm{even}$
(172) $\displaystyle\hskip
64.01869pt=2^{L-6}(L-1)+2^{\frac{L-7}{2}}\left(2^{\frac{L-5}{2}}+1\right),\hskip
8.5359ptL\hskip 8.5359pt\mathrm{odd,}$ $\displaystyle N_{+}=L-2:$
$\displaystyle\hskip 28.45274ptN_{blocks}=1.$ (173)

It would be interesting to find out the analytical formula for the number of
blocks in the transition rate matrix for any general $N_{+}$, which would
contain the formulae in Eq. (173) as special cases.

## References

  * [1] N. G. van Kampen, 1992, Stochastic Processes in Physics and Chemistry (Elsevier, North-Holland).
  * [2] A. Schadschneider, D. Chowdury and K. Nishinari, 2011, Stochastic transport in complex system: From molecules to vehicles (Elsevier, Amstrdam).
  * [3] P. Krapivsky, S. Redner and E. Ben-Naim , 2010, A Kinetic View of Statistical Physics (Cambridge: Cambridge University Press).
  * [4] E. S. Allman and J. A. Rhodes , 2004, Mathematical models in biology : an introduction (Cambridge: Cambridge University Press).
  * [5] R. A. Blythe and A. J. McKane, Stochastic models of evolution in genetics, ecology and linguistics, J. Stat. Mech. P07018 (2007), 10.1088/1742-5468/2007/07/P07018.
  * [6] C. Castellano, S. Fortunato and V. Loreto, Statistical physics of social dynamics, Rev. Mod. Phys. 81, 591 (2009), 10.1103/RevModPhys.81.591.
  * [7] H. Qian and L. M. Bishop, The Chemical Master Equation Approach to Nonequilibrium Steady-State of Open Biochemical Systems: Linear Single-Molecule Enzyme Kinetics and Nonlinear Biochemical Reaction Networks, Int. J. Mol. Sci. 11, 3472 (2010), 10.3390/ijms11093472.
  * [8] B. Schmittmann and R. K. P. Zia, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. Lebowitz (Academic, London, 1994), Vol. 17.
  * [9] V. Privman, 1997, Nonequilibrium Statistical Mechanics in One Dimension (Cambridge University Press).
  * [10] F. Spitzer, Interaction of Markov processes, Adv. in Math. 5, 246 (1970), 10.1016/0001-8708(70)90034-4.
  * [11] T. M. Liggett, 1999, Stochastic Models of Interacting Systems:Contact, Voter and Exclusion Processes, (Springer-Verlag, New-York).
  * [12] G.M. Schütz, in C.Domb and J.Lebowitz (eds.) Phase Transitions and Critical Phenomena, Vol.19 (Academic, London, 2001).
  * [13] D. Helbing, Traffic and related self-driven many-particle systems, Rev. Mod. Phys. 73, 1067 (2001), 10.1103/RevModPhys.73.1067.
  * [14] C. T. MacDonald and J. H. Gibbs, Concerning the kinetics of polypeptide synthesis on polyribosomes, Biopolymers 7, 707 (1969), 10.1002/bip.1969.360070508.
  * [15] K. Nishinari, Y. Okada, A. Schadschneider and D. Chowdury, Intracellular Transport of Single-Headed Molecular Motors KIF1A, Phys. Rev. Lett. 95, 118101 (2005), 10.1103/PhysRevLett.95.118101.
  * [16] C. Leduc, K. Padberg-Gehle, V. Varga, D. Helbing, S. Diez and J. Howard, Molecular crowding creates traffic jams of kinesin motors on microtubules, Proc. Natl Acad. Sci. USA 109, 6100 (2012), 10.1073/pnas.1107281109.
  * [17] B. Derrida, An exactly soluble non-equilibrium system: The asymmetric simple exclusion process, Phys. Rep. 301, 65 (1998), 10.1016/S0370-1573(98)00006-4.
  * [18] T. Sasamoto, One-dimensional partially asymmetric simple exclusion process with open boundaries: orthogonal polynomials approach, J. Phys. A: Math. Gen. 32, 7109 (1999), 10.1088/0305-4470/32/41/306.
  * [19] O. Golinelli and K. Mallick, The asymmetric simple exclusion process: an integrable model for non-equilibrium statistical mechanics, J. Phys. A: Math. Gen. 39, 12679 (2006), 10.1088/0305-4470/39/41/S03.
  * [20] K. Mallick, The exclusion process: A paradigm for non-equilibrium behaviour, Physica A 418, 17 (2015), 10.1016/j.physa.2014.07.046.
  * [21] B. Derrida, M. R. Evans, V. Hakim and V. Pasquier, Exact solution of a 1D asymmetric exclusion model using a matrix formulation, J. Phys. A: Math. Gen. 26, 1493 (1993), 10.1088/0305-4470/26/7/011.
  * [22] M. Uchiyama, T. Sasamoto and M. Wadati, Asymmetric simple exclusion process with open boundaries and Askey–Wilson polynomials, J. Phys. A: Math. Gen. 37, 4985 (2004), 10.1088/0305-4470/37/18/006.
  * [23] F. H. L. Essler and V. Rittenberg, Representations of the quadratic algebra and partially asymmetric diffusion with open boundaries, J. Phys. A: Math. Gen. 29, 3375 (1996), 10.1088/0305-4470/29/13/013.
  * [24] K. Mallick and S. Sandow, Finite-dimensional representations of the quadratic algebra: Applications to the exclusion process, J. Phys. A: Math. Gen. 30, 4513 (1997), 10.1088/0305-4470/30/13/008.
  * [25] B. Derrida, S. A. Janowsky, J. L. Lebowitz and E. R. Speer, Exact Solution of the Totally Asymmetric Simple Exclusion Process: Shock Profiles, J. Stat. Phys. 73, 813 (1993), 10.1007/BF01052811.
  * [26] M. R. Evans, P. A. Ferrari and K. Mallick, Matrix Representation of the Stationary Measure for the Multispecies TASEP, J. Stat. Phys. 135, 217 (2009), 10.1007/s10955-009-9696-2.
  * [27] R. A. Blythe and M. R. Evans, Nonequilibrium steady states of matrix-product form: a solver’s guide, J. Phys. A: Math. Theor. 40, R333 (2007), 10.1088/1751-8113/40/46/R01.
  * [28] P. A. Ferrari and J. B. Martin, Stationary distributions of multi-type totally asymmetric exclusion processes, Ann. Probab. 35(3), 807 (2007), 10.1214/009117906000000944.
  * [29] A. Kuniba, S.Maruyama and M. Okado, Multispecies TASEP and combinatorial $R^{\ast}$, J. Phys. A: Math. Theor. 48, 34FT02 (2015), 10.1088/1751-8113/48/34/34FT02.
  * [30] A. Ayyer and S. Linusson, Correlations in the Multispecies TASEP and a Conjecture by Lam, Trans. Amer. Math. Soc. 369, 1097 (2017), 10.1090/tran/6806.
  * [31] A. Ayyer and S. Linusson, An Inhomogeneous Multispecies TASEP on a Ring, Advances in Applied Mathematics 57, 21 (2014), https://doi.org/10.1016/j.aam.2014.02.001.
  * [32] E. Lee, Integrability of the Multi-Species TASEP with Species-Dependent Rates, Symmetry 13(9), 1578 (2021), 10.3390/sym13091578.
  * [33] N. Crampe, C. Finn, E. Ragoucy and M. Vanicat, Integrable boundary conditions for multi-species ASEP, . Phys. A: Math. Theor. 49, 375201 (2016), 10.1088/1751-8113/49/37/375201.
  * [34] C. Finn, E. Ragoucy and M. Vanicat, Matrix product solution to multi-species ASEP with open boundaries, J. Stat. Mech. 043201 (2018), 10.1088/1742-5468/aab1b5.
  * [35] S. Sandow, Partially asymmetric exclusion process with open boundaries, Phys. Rev. E 50, 2660 (1994), 10.1103/PhysRevE.50.2660.
  * [36] N. Crampé, Algebraic Bethe ansatz for the totally asymmetric simple exclusion process with boundaries, J. Phys. A: Math. Theor. 48, 08FT01 (2015), 10.1088/1751-8113/48/8/08FT01.
  * [37] J. de Gier and F. Essler, Bethe Ansatz Solution of the Asymmetric Exclusion Process with Open Boundaries, Phys. Rev. Lett. 95, 240601 (2005), 10.1103/PhysRevLett.95.240601.
  * [38] M. Kardar, G. Parisi and Y-C. Zhang, Dynamic Scaling of Growing Interfaces, Phys. Rev. Lett. 56, 889 (1986), /10.1103/PhysRevLett.56.889.
  * [39] L. H. Gwa and H. Spohn, Bethe solution for the dynamical-scaling exponent of the noisy Burgers equation, Phys. Rev. A 46, 844 (1992), 10.1103/PhysRevA.46.844.
  * [40] T. Sasamoto and H. Sphon, One-Dimensional Kardar-Parisi-Zhang Equation: An Exact Solution and its Universality, Phys. Rev. Lett. 104, 230602 (2010), 10.1103/PhysRevLett.104.230602.
  * [41] A. Parmeggiani, T. Franosch and E. Frey, Phase Coexistence in Driven One-Dimensional Transport, Phys. Rev. Lett. 90, 086601 (2003), 10.1103/PhysRevLett.90.086601.
  * [42] E. Pronina and A. B. Kolomeisky, Two-channel totally asymmetric simple exclusion processes, J. Phys. A: Math. Gen. 37, 9907 (2004), 10.1088/0305-4470/37/42/005.
  * [43] T. Mitsudo and H. Hayakawa, Synchronization of kinks in the two-lane totally asymmetric simple exclusion process with open boundary conditions, J. Phys. A: Math. Gen. 38, 3087 (2005), 10.1088/0305-4470/38/14/002.
  * [44] R. Jiang, M. Hu, Y. Wu and Q. Wu, Weak and strong coupling in a two-lane asymmetric exclusion process, Phys. Rev. E 77, 041128 (2008), 10.1103/PhysRevE.77.041128.
  * [45] H. J. Hilhorst and C. Appert-Rolland, A multi-lane TASEP model for crossing pedestrian traffic flows, J. Stat. Mech. P06009 (2012), 10.1088/1742-5468/2012/06/P06009.
  * [46] Y. Wang, R. Jiang and Q. Wu, Dynamics in phase transitions of TASEP coupled with multi-lane SEPs, Nonlinear Dyn 88, 1631 (2017), 10.1007/s11071-017-3335-2.
  * [47] P. Helms , U. Ray and G. K.-L. Chan, Dynamical phase behavior of the single- and multi-lane asymmetric simple exclusion process via matrix product states, Phys. Rev. E 100, 022101 (2019), 10.1103/PhysRevE.100.022101.
  * [48] Q.-Y. Hao ,R. Jiang, M.-B. Hu, Y. Zhang, C.-Y. Wu and N. Guo, Theoretical analysis and simulation of phase separation in a driven bidirectional two-lane system, Phys. Rev. E 100, 032133 (2019), 10.1103/PhysRevE.100.032133.
  * [49] I. Neri, N. Kern and A. Permeggiani, Totally Asymmetric Simple Exclusion Process on Networks, Phys. Rev. Lett. 107, 068702 (2011), 10.1103/PhysRevLett.107.068702.
  * [50] K. Mallick, Shocks in the asymmetry exclusion model with an impurity, J. Phys. A: Math. Gen. 29, 5375 (1996), 10.1088/0305-4470/29/17/013.
  * [51] F. H. Jafarpour, Partially Asymmetric Simple Exclusion Model in the Presence of an Impurity on a Ring, J. Phys. A: Math. Gen. 33, 1797 (2000), 10.1088/0305-4470/33/9/306.
  * [52] F. H. Jafarpour, Exact solution of an exclusion model in the presence of a moving impurity on a ring, J. Phys. A: Math. Gen. 33, 8673 (2000), 10.1088/0305-4470/33/48/307.
  * [53] B. Derrida and M. R. Evans, Bethe ansatz solution for a defect particle in the asymmetric exclusion process, J. Phys. A: Math. Gen. 32, 4833 (1999), 10.1088/0305-4470/32/26/303.
  * [54] C. Boutillier, P. Francois, K. Mallick and S. Mallick, A matrix ansatz for the diffusion of an impurity in the asymmetric exclusion process, J. Phys. A: Math. Gen. 35, 9703 (2002), 10.1088/0305-4470/35/46/301.
  * [55] M. R. Evans, Bose-Einstein condensation in disordered exclusion models and relation to traffic flow, Europhys. Lett., 36, 13 (1996), 10.1209/epl/i1996-00180-y.
  * [56] M. J. Lazo and A. J. Ferreira, Asymmetric exclusion model with impurities, Phys. Rev. E 81, 050104(R) (2010), 10.1103/PhysRevE.81.050104.
  * [57] A. B. Kolomeisky, J. Phys. A: Math. Gen. 31, Asymmetric simple exclusion model with local inhomogeneity, 1153 (1998), 10.1088/0305-4470/31/4/006.
  * [58] S. A. Janowsky and J. L. Lebowitz, Finite-size effects and shock fluctuations in the asymmetric simple-exclusion process, Phys. Rev. A 45, 618 (1992), 10.1103/PhysRevA.45.618.
  * [59] G. Tripathy and M. Barma, Driven lattice gases with quenched disorder: Exact results and different macroscopic regimes, Phys. Rev. E 58, 1911 (1998), 10.1103/PhysRevE.58.1911.
  * [60] T. Chou and G. Lakatos, Clustered Bottlenecks in mRNA Translation and Protein Synthesis, Phys. Rev. Lett. 93, 198101 (2004), 10.1103/PhysRevLett.93.198101.
  * [61] M. Barma, Driven diffusive systems with disorder, Physica A 372, 22 (2006), 10.1016/j.physa.2006.05.002.
  * [62] J. Szavits-Nossan, Disordered exclusion process revisited: some exact results in the low-current regime, J. Phys. A: Math. Theor. 46, 315001 (2013), 10.1088/1751-8113/46/31/315001.
  * [63] B. Waclaw, J. Cholewa-Waclaw and P. Greulich, Totally asymmetric exclusion process with site-wise dynamic disorder, J. Phys. A: Math. Theor. 52, 065002 (2019), 10.1088/1751-8121/aafb8a.
  * [64] U. Basu and P. K. Mohanty, TASEP on a ring with internal degrees of freedom, Phys. Rev. E. 82, 041117 (2010), 10.1103/PhysRevE.82.041117.
  * [65] S. Zeraati, F.H. Jafarpour and H. Hinrichsen, Phase transition in an exactly solvable reaction-diffusion process, Phys. Rev. E. 87, 062120 (2013), 10.1103/PhysRevE.87.062120.
  * [66] M. Ghadermazi and F. H. Jafarpour, A new family of exactly solvable disordered reaction-diffusion systems, J. Stat. Mech. P09023 (2013), 10.1088/1742-5468/2013/09/P09023.
  * [67] N. Shiraishi and K. Matsumoto, Undecidability in quantum thermalization, Nature Communications 12, 5084 (2021), 10.1038/s41467-021-25053-0.
  * [68] K. A. Johnson and R. S. Goody, The Original Michaelis Constant: Translation of the 1913 Michaelis-Menten Paper, Biochemistry 50, 8264 (2011), 10.1021/bi201284u.
  * [69] M. J. Schnitzer, K. Visscher and S. M. Block, Force production by single kinesin motors , Nat. Cell Biol. 2, 718 (2000), 10.1038/35036345.
  * [70] R. Grima, Noise-Induced Breakdown of the Michaelis-Menten Equation in Steady-State Conditions, Phys. Rev. Lett. 102, 218103 (2009), 10.1103/PhysRevLett.102.218103.
  * [71] K. Mallick, S. Mallick and N. Rajewsky, Exact solution of an exclusion process with three classes of particles and vacancies, J. Phys. A: Math. Gen. 32, 8399 (1999), 10.1088/0305-4470/32/48/303.
  * [72] K. Klobas, M. Medenjak and T. Prosen, Exactly solvable deterministic lattice model of crossover between ballistic and diffusive transport, J. Stat. Mech. 123202 (2018), 10.1088/1742-5468/aae853.
  * [73] F. H. L. Essler and L. Piroli, Integrability of one-dimensional Lindbladians from operator-space fragmentation, Phys. Rev. E 102, 062210 (2020), 10.1103/PhysRevE.102.062210.
  * [74] P. Sala, T. Rakovszky, R. Verresen, M. Knap and F. Pollmann, Ergodicity Breaking Arising from Hilbert Space Fragmentation in Dipole-Conserving Hamiltonians, Phys. Rev. X 10, 011047 (2020), 10.1103/PhysRevX.10.011047.
  * [75] R. K. P. Zia, E. L. Præstgaard and O. G. Mouritsen, Getting more from pushing less: Negative specific heat and conductivity in nonequilibrium steady states, Am. J. Phys. 70, 384 (2002), 10.1119/1.1427088.
  * [76] R. Eichhorn, P. Reimann, B. Cleuren and C. Van den Broeck, Moving backward noisily, Chaos 15, 026113 (2005), 10.1063/1.1869932.
  * [77] O. Bénichou, P. Illien, G. Oshanin, A. Sarracino and R. Voituriez, , Phys. Rev. Lett. 113, 268002 (2014), 10.1103/PhysRevLett.113.268002.
  * [78] D. Dhar, Diffusion and drift on percolation networks in an external field, J. Phys. A 17, L257 (1984), 10.1088/0305-4470/17/5/007.
  * [79] P. Baerts, U. Basu, C. Maes and S. Safaverdi, Frenetic origin of negative differential response, Phys. Rev. E 88, 052109 (2013), 10.1103/PhysRevE.88.052109.
  * [80] M. Sellitto, Asymmetric Exclusion Processes with Constrained Dynamics, Phys. Rev. Lett. 101, 048301 (2008), 10.1103/PhysRevLett.101.048301.
  * [81] R. Eichhorn, J. Regtmeier, D. Anselmettib and P. Reimann, Negative mobility and sorting of colloidal particles, Soft Matter 6, 1858 (2010), 10.1039/b918716m.
  * [82] M. Baiesi, C. Maes and B. Wynants, Fluctuations and Response of Nonequilibrium States, Phys. Rev. Lett. 103, 010602 (2009), 10.1103/PhysRevLett.103.010602.
  * [83] U. Basu and C. Maes, Mobility transition in a dynamic environment, J. Phys. A: Math. Theor. 47, 255003 (2014), 10.1088/1751-8113/47/25/255003.
  * [84] J. Cividini, D. Mukamel and H. A. Posch, Driven tracer with absolute negative mobility, J. Phys. A: Math. Theor. 51, 085001 (2018), 10.1088/1751-8121/aaa630.
  * [85] A. K. Chatterjee, U. Basu and P. K. Mohanty, Negative differential mobility in interacting particle systems, Phys. Rev. E 97, 052137 (2018), 10.1103/PhysRevE.97.052137.
  * [86] P. F. Arndt, T. Heinzel and V. Rittenberg, Spontaneous breaking of translational invariance in one-dimensional stationary states on a ring, J. Phys. A: Math. Gen. 31, L45 (1998), 10.1088/0305-4470/31/2/001.