# Counter-flow induced clustering: Exact results

Amit Kumar Chatterjee<EMAIL_ADDRESS>Yukawa Institute for
Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku,
Kyoto 606-8502, Japan  Hisao Hayakawa<EMAIL_ADDRESS>Yukawa
Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho,
Sakyo-ku, Kyoto 606-8502, Japan  Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan

###### Abstract

We analyze the cluster formation in a non-ergodic stochastic system as a
result of counter-flow, with the aid of an exactly solvable model. To
illustrate the clustering, a two species asymmetric simple exclusion process
with impurities on a periodic lattice is considered, where the impurity can
activate flips between the two non-conserved species. Exact analytical
results, supported by Monte Carlo simulations, show two distinct phases, free
flowing phase and clustering phase. The clustering phase is characterized by
constant density and vanishing current of the non-conserved species, whereas
the free flowing phase is identified with non-monotonic density and non-
monotonic finite current of the same. The $n$-point spatial correlation
between $n$ consecutive vacancies grows with increasing $n$ in the clustering
phase, indicating the formation of two macroscopic clusters in this phase, one
of the vacancies and the other consisting of all the particles. We define a
rearrangement parameter that permutes the ordering of particles in the initial
configuration, keeping all the input parameters fixed. This rearrangement
parameter reveals the significant effect of non-ergodicity on the onset of
clustering. For a special choice of the microscopic dynamics, we connect the
present model to a system of run-and-tumble particles used to model active
matter, where the two species having opposite net bias manifest the two
possible run directions of the run-and-tumble particles, and the impurities
act as tumbling reagents that enable the tumbling process.

##  I Introduction

Clustering of non-attractive elements is an intriguing phenomenon occurring in
diverse areas of science and society. It appears in various fields such as
granular materials Liu_1998 ; Biroli_2007 ; Behringer_2019  , vehicular and
pedestrian traffic flows Nagel_1992 ; Helbing_2001 ; Nagatani_2002 ;
Schadschneider_2011 , active matter Mandal_2020 ; Yang_2022 , biology
Sadati_2013 ; Oswald_2017  etc. At the heart of the clustering phenomena lies
a transition between free flowing phase and jammed phase where motion becomes
highly restrictive, as some suitable system parameter is tuned. The tuning
parameter for granular material is packing fraction, it is car or pedestrian
density for traffic flow and self-propulsion force in case of dense active
matter. In fact, a jamming phase diagram has been proposed keeping in mind the
generality of such transitions Liu_1998  and there are experiments conducted
on colloidal particles supporting the concept of this generic phase diagram
such as Ref. Trappe_2001 . The control of clustering is important in daily
life. Indeed, understanding the formation of jamming and finding ways to
transit to unjammed state, has immense importance in traffic science. The
other applications include occurrence of cellular jamming transitions in
cancer Sadati_2013 ; Oswald_2017 . Analysis of clustering, as a phase
transition in interacting many body systems, constitutes an interesting topic
from the physics point of view. Thus, it seems natural to investigate the
clustering phenomena through the lens of statistical mechanics.

One possible mechanism of clustering is the presence of counter-flow in a
system. Counter-flows arise naturally in real life situations, e.g. pedestrian
traffic flow Schadschneider_2011  in narrow lanes, busy pedestrian zones,
crossings between two footpaths etc. In these situations, we often observe
pedestrians moving in opposite directions, automatically create counter-flow
situations. Counter-flows in pedestrian dynamics can result in interesting
phenomena like traffic jam at high densities Muramatsu_1999 ; Tajima_2002 ,
self-organized lane formations Helbing_1995 ; Hoogendoorn_2005 , oscillatory
changes in the dominant direction of motion at bottlenecks Kretz_2006_1 ; as
well as unexpected events like panic and crowd disasters Klingsch_2010 .
Interesting experimental results regarding counter-flows have revealed that
the total flow in a counter-flow situation can be greater than the sum of the
two comparable unidirectional flows Hoogendoorn_2005 ; Kretz_2006 . Numerical
studies show that the existence of counter-flowing agents ignoring traffic
rules can surprisingly lower the possibility of traffic jam Baek_2009 . Apart
from the pedestrian dynamics, counter-flow occurs in nature and has been
incorporated in industries. In fact, counter-current exchange of heat or
chemicals between two oppositely flowing fluids (i.e. counter-flow situation),
has been proved to be much more efficient than co-current exchange of the same
between two parallel flowing fluids. Such principles are often used to devise
heat exchangers Sridhar_2017  used in industries and are found in nature e.g.
salt glands in sea and desert birds Schmidt_1959 ; Proctor_1993 , mammalian
kidneys Gottschalk_1958  etc.

In the context of cluster formation resulting from counter flow, spontaneous
symmetry breaking in a model with two oppositely moving particles with
exchange interaction, has been discussed with periodic boundaries Lee_1997 ;
Arndt_1998  and open boundaries Evans_1995 . Notably, later it has been proved
that the spatial condensation of particles under periodic boundary condition
Arndt_1998 , is not associated with a phase transition in the grand canonical
ensemble Rajewsky_2000 ; Sasamoto_2001 . There exist several models in non-
equilibrium statistical mechanics that exhibit phase transitions in one
dimension Schmittmann_1994 ; Privman_1997 ; Marro_1999 ; Schutz_2001 ;
Henkel_2008 ; Henkel_2010 ; Derrida_1993 ; Evans_1998 ; Lahiri_1997 ;
Evans_2002 ; Mallick_1996 ; Derrida_1999 . However, considerably few one
dimensional non-equilibrium models with simplified local dynamics, are
amenable to exact analytical calculations which provide much insights to the
microscopic origins of the phase transitions Derrida_1993 ; Mallick_1996 ;
Derrida_1999 ; Evans_2000 ; Mukamel_2000 . In context of cluster formations in
exactly solvable disordered systems, Bose-Einstein condensates have been
studied in multi-species asymmetric simple exclusion processes Evans_1996 ;
Krug_1996 , utilizing matrix product ansatz Blythe_2007  and zero range
process Spitzer_1970 ; Evans_2005 . Traffic jam in more realistic traffic
models like bus route models, has been considered Loan_1998_1 ; Loan_1998_2
and its connection to Nagel-Schreckenberg model of traffic flow, has been
explored in details Chowdhury_2000 . There are various other interesting
studies of one dimensional traffic flows under periodic boundary conditions
Nagel_1992 ; Helbing_2001 ; Nagatani_2002 ; Schadschneider_2011 ; Bando_1995 ;
Komatsu_1995 ; Hayakawa_1998 ; Igarashi_1999 . Notably, once a formed cluster
is stable and moves in one-direction, it is similar to the time-crystal
Wilczek_2012 ; Shapere_2012 ; Watanabe_2015 ; Yao_2017 ; Else_2020 . To avoid
such clustering, model with bypassing defects through long range hopping, has
been analyzed Nossan_2011 . Notably, exact analytical results showing
clustering phenomenon has been considerably few.

Recently, clustering of self-propelled objects like bacteria have gained much
attention Berg_2004 ; Brumley_2019 . Such active matter can exhibit motility
induced phase separation by accumulating in regions with slow movement
Tailleur_2008 ; Cates_2015 ; Fodor_2017 , or they can aggregate near chemical
nutrients Brumley_2019 . The motion of self-propelled objects including
bacteria, is often described by run-and-tumble particles (RTPs) Thompson_2011
; Soto_2014 ; Slowman_2016 ; Malakar_2018 ; Mallmin_2019 ; Das_2020 ;
Mukherjee_2022 ; Jose_2022 . Analytical solutions of run-and-tumble models are
few, including exact solution for one and two RTPs Slowman_2016 ; Malakar_2018
; Mallmin_2019 ; Das_2020 ; Mukherjee_2022 , mean field analysis of many
interacting RTPs Dandekar_2020  and approximate solution for restricted RTPs
Mukherjee_2022 . The tumbling dynamics corresponds to occasional change in
direction of motion, which specifically in one dimension, would correspond to
flip between right and left moving objects. Thus, it would be interesting to
ask how to accommodate such flip dynamics of RTPs in simple lattice models
that are exactly solvable, contain many interacting RTPs and exhibit
clustering.

In this work, we provide exact results showing cluster formation on a one
dimensional periodic lattice, when counter-flow and flip dynamics are present
in the system. As a model, we consider the two species asymmetric simple
exclusion process with impurity activated flips ($2$-ASEP-IAF) Chatterjee_2022
. We show the existence of two different phases, namely the free flowing phase
and the clustering phase, and characterize them using observables like average
density, drift current and spatial correlations, obtained analytically and
supported by Monte Carlo simulations. Interestingly, in the counter-flow
situation, with a special choice of the microscopic hop rates, the $2$-ASEP-
IAF can be interpreted as a system of many RTPs in presence of tumbling
reagents. Since the RTPs have been extremely useful to model active matter on
lattices, we would discuss separately the clustering in this special case in
some details. Also, the $2$-ASEP-IAF has plausible mappings to two-lane ASEP
with bridges Chatterjee_2022  and enzymatic chemical reactions Chatterjee_2022
, explained briefly in the next section.

In Section II we define the model and briefly describe its connection to the
two-lane model, enzymatic chemical reactions and RTPs. Section III describes
the steps to obtain the exact matrix product steady state of the model with
infinite dimensional matrices. The model being non-ergodic, we state the
choice of initial configuration considered here in Section IV and calculate
the partition function showcasing the important steps. In Section V we discuss
in details the analytical results for density, current and spatial
correlations which establish the existence of two different phases. The role
of initial configurations on the onset of the clustering, is investigated in
Section VI. The special constraint on the microscopic rates, for which the
$2$-ASEP-IAF can be mapped to a model of run-and-tumble particles, is analyzed
in details in Section VII. We summarize our main results and future directions
in Section VIII. Appendix A provides expressions and brief derivations of some
observables. The explicit form of the density-fugacity relation and a special
case that leads to a closed form solution of the fugacity, are presented in
Appendix B. In Appendix C we clarify on the variations of the initial
configuration considered in the main text, that can give rise to similar
clustering. We show the convergence of our results with system size in
Appendix D. The comparison of our model to the Arndt-Heinzel-Rittenberg model
of counter-flow, is briefly discussed in Appendix E.

##  II Model

We consider two different species ($1$ and $2$) along with impurities ($+$)
and vacancies ($0$) on a one dimensional periodic lattice with $L$ sites,
$i=1,2,\dots,L$. Species $1$ and species $2$ can hop to right with rates
$p_{1}$ and $p_{2}$ respectively, or to left with rates $q_{1}$ and $q_{2}$
respectively, if the target site is vacant. The motion of the impurity is
restricted only towards right, with rate $\epsilon$. The unidirectional motion
of the impurity is a characteristic that differentiates it from the two
species. Further, this assumption turns out to be crucial to obtain the exact
steady state probability distribution. Alongside the hopping dynamics, species
$1$ and $2$ can transform into each other with rates $w_{12}$ and $w_{21}$,
such flipping being activated only in presence of impurity. Following the
nomenclature of Ref. Chatterjee_2022 , we refer this model as $2$-ASEP-IAF,
where ASEP represents asymmetric simple exclusion process and IAF stands for
impurity activated flips. The microscopic dynamics of $2$-ASEP-IAF is
represented as

$\displaystyle 10\,\,\xrightleftharpoons[q_{1}]{p_{1}}\,\,01,\hskip
22.76228pt$ $\displaystyle 20\,\,\xrightleftharpoons[q_{2}]{p_{2}}\,\,02,$
$\displaystyle+0\,\,\stackrel{{\scriptstyle\epsilon}}{{\longrightarrow}}\,\,0+,\hskip
22.76228pt$ $\displaystyle 1+\,\,\xrightleftharpoons[w_{21}]{w_{12}}\,\,2+.$
(1)

The $2$-ASEP-IAF can be restated as a four-state model where each lattice site
can be in either of the four possible states $1$ or $2$ or $+$ or $0$. The
input parameter space for the model contains several parameters, more
precisely given by
$(p_{1},p_{2},q_{1},q_{2},\epsilon,w_{12},w_{21},\rho_{0},\rho_{+})$. The
parameters $\rho_{0}=N_{0}/L$ and $\rho_{+}=N_{+}/L$ are the conserved
densities of the vacancies and the impurities, respectively, with $N_{0}$ and
$N_{+}$ being the number of vacancies and impurities in the system. In this
paper, we aim to investigate the effect of counter-flow on the system, and to
do so, we pick out $q_{1}$ to be the tuning parameter, keeping all other rates
fixed. We fix $p_{2}>q_{2}$ (apart from the discussion in Section VII) and
$\epsilon>0$. Subsequently, if $q_{1}<p_{1}$, all particles have net bias
along the same direction, keeping the system in natural flow situation. On the
other hand, when $q_{1}>p_{1}$, species $1$ has net bias in the direction
opposite to the net bias of species $2$ and impurities, thereby species $1$
opposes the motion of other components and creates a counter-flow situation.
Since the tuning parameter $q_{1}$ can control the flow situation in the
system, we denote it as the counter-flow parameter.

We would like to mention three interesting plausible mappings of the present
model [Eq. (1)] to other systems of interest. These three systems, as briefly
explained below, are (i) run-and-tumble particles, (ii) enzymatic chemical
reaction and (iii) two-lane ASEP.

(i) To study active matter on lattices, run-and-tumble particles (RTPs)
constitute useful models Soto_2014 , notably with considerably few analytical
results Slowman_2016 ; Malakar_2018 ; Mallmin_2019 ; Das_2020 ; Mukherjee_2022
; Jose_2022 . The RTPs keep moving along a particular direction until they
tumble i.e. change the direction of motion. The two different species having
net bias towards opposite directions in $2$-ASEP-IAF, can be considered as two
possible orientations of the RTP, only if the counter-flow situation is
considered. Basically the RTP is a two-state particle where its two possible
states or run directions are imitated by species $1$ and $2$ of the $2$-ASEP-
IAF. A particularly relevant case is when $p_{2}-q_{2}=q_{1}-p_{1}$ i.e. the
RTP moves with same speed in both directions with intermediate tumbling, in
accordance with the recently proposed restricted tumbling model Mukherjee_2022
and the continuum active random walk model Jose_2022 . Notably, in our
mapping, the impurities act as the mediators for the tumbling process. This
special case which maps the $2$-ASEP-IAF to RTPs, would be discussed in more
details in Section VII.

(ii) The $2$-ASEP-IAF can be connected to an enzymatic chemical reaction in a
narrow channel of diffusing chemical reagents. In this case, the impurity
represent the enzymes which initiate the reaction between substrate and
products, that amounts to the flip dynamics in the present model
Chatterjee_2022 .

(iii) The $2$-ASEP-IAF can be mapped to a two-lane ASEP, a simple model for
two-lane traffic flow, where the species $1$ and $2$ play the roles of
particles hopping with different rates in the two different lanes and the
impurities in $2$-ASEP-IAF mimic the bridges connecting the two lanes in the
two-lane ASEP Chatterjee_2022 . Thus the flip dynamics between two species
activated by impurities represent the lane change dynamics of the particles in
the two-lane ASEP model.

Due to such connections of the $2$-ASEP-IAF to several other important models,
we would like to study analytically the possibility of cluster formation in
the model in presence of counter-flow, which can possibly provide information
about clustering in the connected models.

##  III Steady state

The probability $P(\left\\{s_{i}\right\\})$ of any configuration
$\left\\{s_{i}\right\\}$ ($s_{i}=1$ or $2$ or $0$ or $+$, denoting the
constituent at site $i$), in the non-equilibrium steady state corresponding to
the microscopic dynamics Eq. (1), is obtained in the following matrix product
form

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

In Eq. (3), the matrix $X_{i}$ represents the component $s_{i}$ at site $i$
and $\delta_{(}.,.)$ is Kronecker delta symbol. Specifically, the matrices
$D_{1},D_{2},E,A$ correspond to species $1$, species $2$, vacancy and
impurity, respectively. The configurations of the system evolve according to
the Master equation

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

where the matrix $M$ is the rate matrix containing transition rates between
configurations and $|P(t)\rangle$ is the column vector whose elements are time
dependent probabilities $P\left(\left\\{s_{i}\right\\},t\right)$ for all
possible configurations $\left\\{s_{i}\right\\}$. In the steady state, the
probabilities $P\left(\left\\{s_{i}\right\\},t\right)$ converge to the time
independent values $P\left(\left\\{s_{i}\right\\}\right)$ mentioned in Eq.
(3). For $2$-ASEP-IAF with two-site local dynamics [Eq. (1)], the rate matrix
can be decomposed 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 $16\times 16$ matrix acting on the pair of
sites $(i,i+1)$ and $I$ is $4\times 4$ identity matrix placed at every site
except the pair $(i,i+1)$. In steady state $M|P\rangle=0$. The steady state
can be achieved using the following two-site flux-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}$ $\displaystyle=$
$\displaystyle\left(E,A,D_{1},D_{2}\right)^{T},$ (7)
$\displaystyle\tilde{\mathbf{X}}$ $\displaystyle=$
$\displaystyle\left(\tilde{E},\tilde{A},\tilde{D}_{1},\tilde{D}_{2}\right)^{T},$
(8)

where $(.)^{T}$ denotes the transpose of the row vector $(.)$ and
$\tilde{E},\tilde{A},\tilde{D}_{1}.\tilde{D}_{2}$ 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,A,D_{1},D_{2}$.
We find that suitable choices for the auxiliary matrices in this case, are

$\displaystyle\tilde{E}=1,\,\,\tilde{A}=0,\,\,\tilde{D}_{1}=0,\,\,\tilde{D}_{2}=0.$
(9)

Using the above choices of the auxiliary matrices in Eq. (6), we arrive at the
following matrix algebra

$\displaystyle p_{1}D_{1}E-q_{1}ED_{1}$ $\displaystyle=$ $\displaystyle
D_{1},$ (10) $\displaystyle p_{2}D_{2}E-q_{2}ED_{2}$ $\displaystyle=$
$\displaystyle D_{2},$ (11) $\displaystyle\epsilon AE$ $\displaystyle=$
$\displaystyle A,$ (12) $\displaystyle w_{12}D_{1}A$ $\displaystyle=$
$\displaystyle w_{21}D_{2}A.$ (13)

We find that Eq. (13) is satisfied by the infinite dimensional matrix
representations given below,

$\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)$ (26) $\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$ (35) $\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},\hskip
14.22636ptd_{2}^{1,1}=w_{12}.$

The notation $(m)_{r}$ corresponds to Pochhammer symbol for rising factorials,
$(m)_{r}:=m(m+1)(m+2)\dots(m+r-1)$. In Eq. (35), the matrices $D_{1}$ and
$D_{2}$, representing the two non-conserved species, are upper triangular
matrices, whose elements involve their corresponding hop rates and the flip
rate from the other species (e.g. $D_{1}$ involves $p_{1},q_{1},w_{21}$). The
impurity is represented by the matrix $A$ which has nonzero elements in the
first row only, the elements being functions of the impurity hop rate
$\epsilon$. The matrix $E$, characterizing vacancy, is a lower shift matrix.

We should mention that the matrix representations for any number of non-
conserved species $\mu$ ($\mu=2$ for the present discussion) i.e. $\mu$-ASEP-
IAF, has been obtained recently in Ref. Chatterjee_2022 . However, we would
like to emphasize that the exact analytical calculations of observables and
their properties strongly depend on the choice of initial configurations,
owing to the non-ergodicity resulting from the microscopic dynamics Eq. (1).
Particularly, the initial configuration considered in Ref. Chatterjee_2022  is
only a special case of the one that would be discussed in the next section
(Section IV). In fact, interestingly, as we would show later, the onset of
clustering and its demarcation from the free flowing phase, crucially depend
on the choice of initial configuration (Section VI).

##  IV Initial configuration and partition function

In spite of the presence of flip-dynamics, the microscopic dynamics in Eq. (1)
preserves certain orderings of the different species and impurities from
initial configuration. It implies that the system is non-ergodic and, can
access only a subspace of the whole configuration space, starting from a
particular initial configuration. To demonstrate the cluster formation in
$2$-ASEP-IAF, we consider the following initial configuration (represented by
corresponding matrices),

$D_{2}A\dots D_{2}A\,\,D_{1}A\dots D_{1}A\,\,D_{1}\dots D_{1}\,\,E\dots E.$
(36)

The dots in $\mathcal{Y}\dots\mathcal{Y}$ ($\mathcal{Y}=D_{2}A,D_{1}A,D_{1}$
and $E$) represent an uninterrupted sequence of the unit $\mathcal{Y}$. Note
that, there are two types of species $1$ particles in Eq. (36). One type is
those which can flip, belongs to the sequence $\mathcal{Y}=D_{1}A$, while the
others are non-flipping species $1$ particles belonging to the sequence
$\mathcal{Y}=D_{1}$. The density of such non-flipping species $1$ particles is
denoted by

$\bar{\rho}:=\frac{\bar{N}}{L},$ (37)

where the number of particles of species $1$ that cannot flip to species $2$
is given by $\bar{N}$. We do not require to introduce such symbol for species
$2$ since the number of non-flipping species $2$ particles is zero in Eq.
(36). A careful look at the initial configuration [Eq. (36)] and the dynamics
[Eq. (1)] reveals that each configuration in the system can be identified as a
sequence of intervals between impurities. Each interval contains one species
$1$ or species $2$ particle and vacancies, except one interval that contains
additional $\bar{N}$ non-flipping species $1$ particles apart from the one
particle that can flip. The species $1$ and species $2$ particles can hop to
right or left within this interval and can flip into each other at the two
boundaries (i.e. impurities) of the interval. Each interval can increase or
decrease in size by the incoming or outgoing vacancies at the right and left
interval boundaries, respectively. This effectively gives rise to a cyclic
motion of the intervals towards right.

For convenience, we calculate the partition function corresponding to the
initial configuration Eq. (36) in the grand canonical ensemble, by associating
the fugacity $z_{0}$ with the vacancies ($0$). Consequently, the partition
function would be

$\displaystyle
Q=\sum_{m_{1}=0}^{\infty}\dots\sum_{m_{N_{+}}=0}^{\infty}\,\,\sum_{\bar{m}_{1}=0}^{\infty}\dots\sum_{\bar{m}_{N_{+}}=0}^{\infty}\,\,\sum_{n_{1}=0}^{\infty}\dots\sum_{n_{\bar{N}}=0}^{\infty}\mathrm{Tr}$
(38)
$\displaystyle\left[\left(\prod_{i=1}^{N_{+}}(D_{1}+D_{2})(z_{0}E)^{m_{i}}A(z_{0}E)^{\bar{m}_{i}}\right)\left(\prod_{j=1}^{\bar{N}}D_{1}(z_{0}E)^{n_{j}}\right)\right].$

The emergence of the products inside the trace in Eq. (LABEL:eq:pf1) can be
understood from the interpretation of each configuration as a sequence of
intervals between impurities discussed above. Since we have $N_{+}$ impurities
in the periodic system, that leads to $N_{+}$ intervals each bounded by
impurities at both ends. This corresponds to the product of $N_{+}$ intervals
in the first term inside the trace in Eq. (LABEL:eq:pf1). Each product in the
first term inside the trace correspond to interval containing one particle
that can flip, whereas the last product term denotes the sequence of $\bar{N}$
non-flipping species $1$ particles inside the one exceptional interval. We
would use the explicit matrix representations from Eq. (35) to calculate the
partition function in Eq. (LABEL:eq:pf1). It is suitable to write down the
matrices in the concise form as

$\displaystyle E$ $\displaystyle=$
$\displaystyle\sum_{\gamma=1}^{\infty}|\gamma+1\rangle\langle\gamma|\Rightarrow
E^{n}=\sum_{\gamma=1}^{\infty}|\gamma+n\rangle\langle\gamma|,$ (40)
$\displaystyle D_{1,2}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{\infty}\sum_{\beta=\alpha}^{\infty}(d_{1,2})_{\alpha,\beta}|\alpha\rangle\langle\beta|,$
(41) $\displaystyle(d_{1,2})_{\alpha,\beta}$ $\displaystyle:=$
$\displaystyle\frac{(\beta-1)!}{(\alpha-1)!(\beta-\alpha)!}\,\frac{q_{1,2}^{\alpha-1}}{p_{1,2}^{\beta-1}},$
(42) $\displaystyle A$ $\displaystyle=$
$\displaystyle\sum_{\delta=1}^{\infty}\frac{1}{\epsilon^{\delta-1}}|1\rangle\langle\delta|,$
(43)

where $\langle k|=(0,0,\dots 1,\dots 0)$ is a standard basis vector with only
non-zero element $1$ at the $k$-th place and $|k\rangle=(0,0,\dots 1,\dots
0)^{T}$ where the superscript $T$ denotes transpose of the vector under
consideration. With the above expressions, we can simplify the matrix strings
in Eq. (LABEL:eq:pf1). For example, we obtain

$\displaystyle\prod_{j=1}^{\bar{N}}D_{1}(z_{0}E)^{n_{j}}=\sum_{n_{1}}\dots\sum_{n_{\bar{N}}}\,\,\sum_{\alpha_{1}}\sum_{\alpha_{2}}\dots\sum_{\alpha_{\bar{N}}}\,\sum_{\beta_{\bar{N}}}\,z_{0}^{n_{1}+\dots+n_{\bar{N}}}\,$
(44)
$\displaystyle(d_{1})_{\alpha_{1},\alpha_{2}+n_{1}}\dots(d_{1})_{\alpha_{\bar{N}-1},\alpha_{\bar{N}}+n_{\bar{N}-1}}\,(d_{1})_{\alpha_{\bar{N}},\beta_{\bar{N}}+n_{\bar{N}}}|\alpha_{1}\rangle\langle\beta_{\bar{N}}|.$
(45)

We incorporate the explicit form of $(d_{1})_{\alpha,\beta}$ from Eq. (43) to
evaluate the sums in Eq. (45). However, it is instructive to perform the above
calculation recursively, i.e. first for single $D_{1}$, then two $D_{1}$,
followed by three $D_{1}$ and finally generalize the result for $\bar{N}$
$D_{1}$-s by noting the trend. On the other hand, after invoking Eq. (43), the
other string of matrices in Eq. (LABEL:eq:pf1) reduces to

$\displaystyle\prod_{i=1}^{N_{+}}(D_{1}+D_{2})(z_{0}E)^{m_{i}}A(z_{0}E)^{\bar{m}_{i}}=$
(46)
$\displaystyle\left(\sum_{m_{1}}z_{0}^{m_{1}}\sum_{\alpha_{1}}\left((d_{1})_{\alpha_{1},1+m_{1}}+(d_{2})_{\alpha_{1},1+m_{1}}\right)|\alpha_{1}\rangle\right)\times$
(47)
$\displaystyle\left(\sum_{m}\sum_{n}\frac{z_{0}^{m+n}}{\epsilon^{n-1}}\sum_{\alpha}\left(\frac{(d_{1})_{\alpha,1+m}}{\epsilon^{\alpha}}+\frac{(d_{2})_{\alpha,1+m}}{\epsilon^{\alpha}}\right)\right)^{N_{+}-1}\times$
(48)
$\displaystyle\sum_{\bar{m}_{N_{+}}}z_{0}^{\bar{m}_{N_{+}}}\sum_{\delta_{N_{+}}}\frac{\langle\delta_{N_{+}}|}{\epsilon^{\delta_{N_{+}}-1+\bar{m}_{N_{+}}}}.$
(49)

Using Eqs. (45) and (49) in Eq. (LABEL:eq:pf1), along with the explicit forms
of $(d_{1,2})_{\alpha,\beta}$ from Eq. (43), we finally arrive at the
following expression of the partition function

$\displaystyle
Q=\left(\frac{1}{1-\frac{z_{0}}{\epsilon}}\right)^{N_{+}}\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}-\frac{z_{0}}{\epsilon}\frac{q_{1}}{p_{1}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}-\frac{z_{0}}{\epsilon}\frac{q_{2}}{p_{2}}}\right]^{N_{+}-1}$
(50)
$\displaystyle\times\prod_{k=1}^{\bar{N}}\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{k-1}}\times\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{\bar{N}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}-\frac{q_{2}}{p_{2}}\frac{z_{0}}{p_{1}}S_{\bar{N}-1}}\right],$
(51)

where

$\displaystyle S_{k}$ $\displaystyle:=$
$\displaystyle\sum_{j=0}^{k}\left(\frac{q_{1}}{p_{1}}\right)^{j}+\left(\frac{q_{1}}{p_{1}}\right)^{k}\frac{q_{1}}{\epsilon}$
(52) $\displaystyle=$
$\displaystyle\frac{p_{1}\left[\left(\frac{q_{1}}{p_{1}}\right)^{k+1}((k+1)(p_{1}-q_{1})-\epsilon)+\epsilon\right]}{(p_{1}-q_{1})\epsilon}.$
(53)

The fugacity $z_{0}$ is computed from the density-fugacity relation

$\rho_{0}=\frac{z_{0}}{L}\frac{d}{dz_{0}}\mathrm{ln}Q.$ (54)

Thus, we have obtained the analytical form of the partition function Eq. (53)
corresponding to the initial configuration Eq. (36).

Figure 1: The figures (a) and (b) show heat-maps for density $\rho_{1}$ and
current $J_{10}$ respectively, for species $1$, in the $p_{1}-q_{1}$ parameter
plane. Both of them, particularly figure (b) clearly exhibits that
$q_{1}=p_{1}$ separates two different phases. Figure (a) also implies the non-
monotonicity of the density. The parameters used are
$L=10^{3},p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

##  V Two phases: analytical results

In this section, we would discuss exact results for observables, showing the
emergence of two different phases in the $2$-ASEP-IAF, with the variation of
the parameter $q_{1}$. Starting from the initial configuration Eq. (36), the
analytical calculations are performed following the footsteps sketched in the
previous section, where we use the explicit representations of the matrices
from Eq. (35). Before entering into the detailed discussions of the
observables, below we briefly summarize our main findings.
(i) The analytical results for one-point (average species densities), two-
point (drift current) and $n$-point functions (correlation between consecutive
vacancies) exhibit the existence of two-different phases: the free flowing
phase and the clustering phase. In Figs. 1(a) and (b), we present the heat
maps for average species density and average drift current (for species $1$)
respectively, in the parameter plane $p_{1}-q_{1}$. Fig. 1(b) prominently
distinguishes two phases separated by the line $q_{1}=p_{1}$. Fig. 1(a) also
shows two distinct phases demarcated by a region around the line
$q_{1}=p_{1}$.
(ii) The average species density (of species $1$), when plotted against
$q_{1}$, remarkably exhibits two different behaviors for $q_{1}<p_{1}$ and
$q_{1}>p_{1}$. For $q_{1}<p_{1}$, defining the free flowing phase, the density
is non-monotonic, whereas it remains constant for $q_{1}>p_{1}$ that
constitutes the clustering phase. Thus the non-conservation plays an important
role in identifying the two different phases through a simple one point
function.
(iii) The average drift current (of species $1$) is non-zero in the free-
flowing phase, followed by a sharp descent toward zero near $q_{1}=p_{1}$, and
remains vanishingly small in the clustering phase.
(iv) The $n$-point correlation function between $n$ consecutive vacancies
increases with increasing $n$ in the counter-flow phase, directly pointing
towards the formation of vacancy cluster, and consequently that of particle
cluster.

###  V.1 Average species density

We first consider the one-point functions, i.e. the average species densities
$\rho_{1}=\langle 1\rangle$ and $\rho_{2}=\langle 2\rangle$ of the non-
conserved species $1$ and $2$, respectively, where $\langle.\rangle$ denotes
ensemble average in the steady state. Since, the total density of the two
species remain constant ($\rho_{1}+\rho_{2}=\bar{\rho}+\rho_{+}$), it would be
sufficient to focus only on the density behavior of one of the species, say
$\rho_{1}$. The formal expression for $\langle 1\rangle$ that has been used to
calculate $\rho_{1}$, is presented in Eq. (LABEL:eq:ob1) of Appendix A. With
the aid of Eqs. (45) and (49), the Eq. (LABEL:eq:ob1) leads to the following
exact formula

$\displaystyle\rho_{1}$ $\displaystyle=$
$\displaystyle\bar{\rho}+\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{0}}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{0}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}\left(1+\frac{q_{2}}{p_{1}}S_{-1}\right)}\right]}$
(56)
$\displaystyle+\frac{1}{L}\frac{\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{\bar{N}}}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{\bar{N}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}\left(1+\frac{q_{2}}{p_{1}}S_{\bar{N}-1}\right)}\right]},$

where $S_{k}$ ($k=0$ or $-1$ or $N$ or $N-1$) is given by Eq. (53). To express
the density in Eq. (56) as a function of the input parameters only, it is
required to solve for $z_{0}$, the density-fugacity relation Eq. (54) which
more explicitly takes the form Eq. (99). However, the closed form solution of
Eq. (99) is difficult to obtain. The solution of $z_{0}$ obtained from
Mathematica, is used to replace it in Eq. (56). We present the corresponding
analytical results (solid line) for $\rho_{1}$, agreeing with Monte Carlo
simulation results (dots) as a function of $q_{1}$ in Fig. 2. It is
fascinating that even the simplest one-point function, namely the average
species density, clearly indicates the existence of two different phases. In
one phase, occurring in the parameter region $q_{1}<p_{1}$, the density
exhibits a non-monotonic behavior. Whereas it remains constant in the other
phase, characterizing the parameter region $q_{1}>p_{1}$. As one increases
$q_{1}$ starting from zero, the hopping of species $1$ particles to left,
becomes increasingly likely. This means lesser chances for species $1$
particles to have impurities as their right nearest neighbors. Consequently,
the flipping of species $1$ to $2$ decreases with increasing $q_{1}$, and
therefore $\rho_{1}$ increases. After reaching the maximum density, $\rho_{1}$
starts decreasing, that results in a steep fall near $q_{1}\lesssim p_{1}$.
For $q_{1}>p_{1}$, the density remains constant indicating that the drift
process no longer can affect $\rho_{1}$, meaning the species $1$ particles
cannot access vacancies due to possible clustering. The point $q_{1}=p_{1}$
demarcating two different phases, is regarded as the transition point. Note
that the Fig. 2 is a cross section of Fig. 1(a) at $p_{1}=0.4$, and further
clarifies the behavior of $\rho_{1}$ observed in the heat-map in Fig. 1(a).

Figure 2: The figure shows the behavior of the average species density as a
function of the counter-flow parameter, obtained from analytical result (brown
solid line) and Monte Carlo simulation (red dots). Two different phases are
apparent, separated by the transition point $q_{1}=p_{1}$. For $q_{1}<p_{1}$,
the density exhibits non-monotonicity whereas for $q_{1}>p_{1}$ the density
remains constant. The inset presents the variation of $z_{0}$ with $q_{1}$.
The fugacity exhibits a sharp decrease towards zero in $q_{1}\lesssim p_{1}$
and remains close to zero for $q_{1}>p_{1}$.

The parameters used are
$L=10^{3},p_{1}=0.4,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$..

To see how the two different behaviors of $\rho_{1}$, non-monotonic (in
$q_{1}<p_{1}$) and constant (in $q_{1}>p_{1}$), are influenced by that of the
fugacity $z_{0}$, we show the variation of $z_{0}$ with $q_{1}$ in the inset
of Fig. 2. The inset shows that $z_{0}$ has a sharp decrease near
$q_{1}=p_{1}$ and it is almost zero in the regime $q_{1}>p_{1}$. Consequently,
we can take the limit $z_{0}\rightarrow 0$ in Eq. (56) which would give us the
working formula for $\rho_{1}$ in the constant density phase. This leads us to
the following formula

$\displaystyle\rho_{1}=\bar{\rho}+\frac{w_{21}}{w_{12}+w_{21}}\rho_{+},\hskip
11.38092pt\mathrm{for}\,\,q_{1}>p_{1}.$ (57)

The Eq. (57) clearly manifests that $\rho_{1}$ essentially becomes independent
of the parameter $q_{1}$ in the regime $q_{1}>p_{1}$. The exact above
expression for the average species density can be obtained in a situation
where we have two species $1$ and $2$ along with impurities in a system, but
no vacancies, meaning the only dynamics is
$1+\xrightleftharpoons[w_{21}]{w_{12}}2+$ and there is no drift since there
are no vacancies. Such analysis strongly points towards the existence of a
single cluster consisting of all the particles accompanied by another cluster
consisting of only vacancies in our system in the parameter region
$q_{1}>p_{1}$. Consequently, this phase is referred as the clustering phase.
On the other hand, for $q_{1}<p_{1}$, the drift process is significant with
flows of particles and vacancies, thereafter named as the free flowing phase.
In the free flowing phase, near the transition point i.e. $q_{1}\lesssim
p_{1}$, we find from Eq. (56) that the fall of $\rho_{1}$ towards the constant
value [Eq. (57)] takes the form below

$\displaystyle\rho_{1}$ $\displaystyle\approx$
$\displaystyle\bar{\rho}+\frac{w_{21}}{w_{12}+w_{21}}\rho_{+}$ (59)
$\displaystyle+\,z_{0}\,\frac{w_{21}w_{12}}{(w_{12}+w_{21})^{2}}\left(A_{1}+A_{2}\right)\hskip
5.69046pt\mathrm{for}\,\,q_{1}\lesssim p_{1},$

where
$A_{1}=\left(\frac{1}{p_{1}}+\frac{1}{\varepsilon}-\frac{1}{p_{2}}-\frac{q_{2}}{p_{2}\varepsilon}\right)\left(\rho_{+}-\frac{1}{L}\right)$
and
$A_{2}=\left(\frac{1}{p_{1}}+\frac{1}{\varepsilon}\right)\left(\bar{\rho}(1-\frac{q_{2}}{p_{2}})+\frac{1}{L}\right)-\frac{1}{p_{2}L}$,
we have neglected higher orders of $z_{0}$. It is evident from Eq. (59) that
the fall of $\rho_{1}$ in the regime $q_{1}\lesssim p_{1}$ is linear in
$z_{0}$ and would depend on how $z_{0}$ falls to zero in $q_{1}\lesssim p_{1}$
as a function of $(p_{1}-q_{1})$.

Figure 3: The figure demonstrates two different phases manifested by the drift
current, obtained from analytical result (orange solid line) and Monte Carlo
simulation (blue dots). In the parameter regime $q_{1}<p_{1}$, the current is
finite whereas it falls abruptly to vanishingly small values at $q_{1}=p_{1}$
and remains so for $q_{1}>p_{1}$. The parameters used are
$L=10^{3},p_{1}=0.4,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

###  V.2 Average drift current

The drift currents for the non-conserved species and the impurity are given by

$\displaystyle J_{10}=p_{1}\langle 10\rangle-q_{1}\langle
01\rangle=z_{0}\rho_{1},$ (60) $\displaystyle J_{20}=p_{2}\langle 10\rangle-
q_{2}\langle 01\rangle=z_{0}\rho_{2},$ (61) $\displaystyle
J_{+0}=\epsilon\langle+0\rangle=z_{0}\rho_{+}.$ (62)

The average steady state densities of the non-conserved species are provided
in Eq. (56) whereas the impurity density $\rho_{+}$ is a conserved quantity.
The fugacity $z_{0}$ obeys Eq. (54), more precisely, is the solution of Eq.
(99). The analytical result for $J_{10}$ has been presented as a heat-map in
Fig. (1)(b) which prominently exhibits the existence of two difference phases,
demarcated by the line $q_{1}=p_{1}$. In Fig. 3, we present a cross section of
Fig. (1)(b) at $p_{1}=0.4$. We observe that the current, starting from non-
zero value in the free flowing phase ($q_{1}<p_{1}$), shows a weak increase.
This is due to the increasing density in this parameter regime (Fig. 2).
However, the slope of the increase of $J_{10}$ (Fig. 3) is much smaller than
that of $\rho_{1}$ (Fig. 2), because $J_{10}=z_{0}\rho_{1}$ and $z_{0}$
decreases with $q_{1}$ in the same parameter regime (inset of Fig. 2). After
reaching the maximum, $J_{10}$ exhibits a sharp fall near the transition point
$q_{1}=p_{1}$, and it remains almost zero (i.e. vanishingly small values) in
the clustering phase ($q_{1}>p_{1}$). This observation is consistent with the
emergence of the cluster formation for $q_{1}>p_{1}$ corresponding to the
approximate situation where the vacancies merely play any role i.e. the drift
process is almost absent. As mentioned in the previous section, this gives
rise to two macroscopic clusters, one consisting of all particles and the
other made of only vacancies. With the drift becoming negligible in the
clustering phase, the density remains constant [Fig. 2, Eq. (57)] and the
current stays vanishingly small [Fig. 3]. The species $2$ and impurity would
have similar characteristics of current, therefore we restrict our
investigation to $J_{10}$ only.

From the analysis of the average species density and average current, we
understand that the $2$-ASEP-IAF can exist in two distinct phases, the free
flowing phase characterized by non-monotonic density and non-zero current, and
the clustering phase identified by constant density and vanishing current. It
is important to note that the transition occurs at $q_{1}=p_{1}$ meaning that
sufficiently small counter-flow is enough to impose clustering in the system.

###  V.3 Cluster formation

In the previous sections V.1 and V.2, we have inferred the formation of
clusters from the features of the density and the current. In this section, we
discuss the physical origin of the clustering in $2$-ASEP-IAF and provide
direct observational evidence of the cluster formation. To analyze the reason
behind clustering, we remember that every configuration of the system can be
interpreted as a sequence of intervals, as discussed in section IV. Each
interval has two impurities at the interval boundaries and particles (species
$1$ and species $2$) that drift within intervals and flip at the boundaries.
The total number of particles within each interval is conserved. The intervals
can drift towards right and increase or decrease in size as vacancies can
enter into and exit from these intervals. The non-flipping species $1$
particles inside the rightmost interval presented in Eq. (36), form a dense
region against the impurity to their left. This is because the impurity can
only move to right as opposed to the non-flipping species $1$ sequence that
prefers to move left. For each interval, the vacancies can enter the interval
only through right boundary whereas they can exit from the interval only
through left boundary. Thus, whichever way the vacancies are distributed
initially within the intervals, at long time they finally exit through the
left boundaries of all the intervals and aggregate between the left boundary
of leftmost interval and the accumulated species $1$ particles in the
rightmost interval. Such clustering of vacancies also imply the formation of a
cluster of all particles (species $1$ and species $2$ and impurities). Once
the two clusters (particles and vacancies) are formed, the only way for re-
entrance of the vacancies inside the particle cluster is by hopping through
the accumulated non-flipping species $1$ particles. This requires right
hopping of these species $1$ particles which is less probable as they are left
biased. As the sequence of non-flipping species $1$ particles becomes larger,
it takes longer time for a vacancy to travel thorough the whole particle
cluster. It points towards a slow shift of the particle cluster to right and
we expect the velocity of the shift to become smaller with increasing
$\bar{\rho}$ [Eq. (37)]. We denote the average velocity of the particle
cluster by $v_{\mathrm{cl}}$. From the expressions for currents (as product of
density and velocity) in Eq. (62), the fugacity $z_{0}$ [Eq. (54)] is
identified to $v_{\mathrm{cl}}$.

Figure 4: The figure shows the formation of two macroscopic clusters of
particles and vacancies for $q_{1}>p_{1}$, obtained from Monte Carlo
simulation. The clusters remain stationary in the shifted reference frame
$j^{\prime}=j-v_{\mathrm{cl}}\tau$, implying the drift of cluster with
velocity $v_{\mathrm{cl}}$ in the original frame of reference. The cluster
velocity $v_{\mathrm{cl}}$ decreases with increasing $\bar{\rho}$ [Eq. (37)],
shown in the inset. The parameters used are
$L=10^{3},p_{1}=0.4,q_{1}=0.5,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.29$
and $\rho_{0}=1-2\rho_{+}-\bar{\rho}$.

Next we show the clustering and the drift of the particle cluster explicitly.
After the cluster is formed, we start our observation from any configuration
and note the position of the particle cluster (determined by the impurity at
the leftmost end of the particle cluster) at $\tau=0$. In Monte Carlo
simulation, we measure the position of the cluster at different $\tau$. After
time $\tau$, the cluster should shift by an amount of $v_{\mathrm{cl}}\tau$.
It would be convenient to consider a frame of reference where the cluster
remains stationary. Such stationary cluster frame of reference can be
conceived by using the shifted lattice index
$j^{\prime}(\tau):=j-v_{\mathrm{cl}}\tau$. To characterize the position of
particles and vacancies in any configuration, we consider the following
observable in the shifted reference frame:

$\displaystyle C_{j^{\prime}(\tau)}$ $\displaystyle=$ $\displaystyle 1\hskip
14.22636pt\mathrm{if}\,\,s_{j}^{\prime}=1\,\,\mathrm{or}\,\,s_{j}^{\prime}=2\,\,\mathrm{or}\,\,s_{j}^{\prime}=+$
(63) $\displaystyle=$ $\displaystyle 0\hskip
14.22636pt\mathrm{if}\,\,s_{j}^{\prime}=0.$ (64)

From the Monte Carlo simulation results in Fig. 4, we observe that two
macroscopic clusters are formed in the system at any $\tau$, one consisting of
all the particles and the other made up of all the vacancies. Further, the
cluster is observed to form at the same position in the shifted frame of
reference, for different $\tau$. This implies that the cluster indeed moves
with velocity $v_{\mathrm{cl}}$ in the original frame of reference. The slow
shift of the cluster is similar to time crystals Wilczek_2012 ; Shapere_2012 ;
Watanabe_2015 ; Yao_2017 ; Else_2020 . Note that the slow shift in our model
is caused by the unidirectional motion of impurities, and thus the emergence
of this slow shift is not contradictory with the no-go theorem of time-
crystals Watanabe_2015 . The inset of Fig. 4 shows the variation of
$v_{\mathrm{cl}}$ as a function of the density $\bar{\rho}$ [Eq. (37)] of the
non-flipping species $1$ particles. We observe that, with increasing
$\bar{\rho}$, the cumulative effect of non-flipping species $1$ particles
(prone towards left hopping) increases and thereby the cluster velocity
(towards right) decreases.

###  V.4 Spatial correlation

An illuminating way to analytically show the formation of clusters, is to
calculate the $n$-point correlation between consecutive vacancies (it is more
helpful than calculating correlations between particles, because we have
mixture of different species and impurities inside the particle cluster). The
formal expression for the expression for $n$-point correlation function
between $n$ consecutive vacancies, is the following

$C^{[n]}_{0}=\langle 00\dots 0\rangle-\langle 0\rangle^{n}.$ (65)

The superscript $[n]$ in Eq. (65) represents the length of the uninterrupted
sequence of consecutive vacancies (denoted by $0$ in the subscript). In fact,
the above correlation can be obtained recursively starting from the two-point
correlation $C^{[2]}_{0}=\langle 00\rangle-\rho_{0}^{2}$. We obtain the
following expression for the two point correlation (see Appendix A for more
details),

$\displaystyle
C^{[2]}_{0}=\rho_{0}-\rho_{0}^{2}-\left(\rho_{+}-\frac{1}{L}\right)\frac{z_{0}}{\epsilon}-\frac{1}{L}\frac{z_{0}}{p_{1}}\sum_{k=0}^{\bar{N}-1}S_{k}$
(66)
$\displaystyle-\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}X_{1}}{1-X_{1}}+\frac{w_{12}X_{2}}{1-X_{2}}}{\frac{w_{21}}{1-X_{1}}+\frac{w_{12}}{1-X_{2}}}-\frac{1}{L}\frac{\frac{w_{21}Y_{1}}{1-Y_{1}}+\frac{w_{12}Y_{2}}{1-Y_{2}}}{\frac{w_{21}}{1-Y_{1}}+\frac{w_{12}}{1-Y_{2}}},$
(67)

Figure 5: The figure shows abrupt change in the $n$-point correlation between
consecutive vacancies near the transition point $q_{1}=p_{1}$. For
$q_{1}>p_{1}$, the growth of the correlation with increasing $n$ implies the
tendency of vacancies to form larger cluster. Solid lines and dots correspond
to analytical calculation and Monte Carlo simulation results respectively. The
parameters used are
$L=10^{3},p_{1}=0.4,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

where $S_{k}$ is given by Eq. (53) and
$X_{1}=\frac{z_{0}}{p_{1}}+\frac{z_{0}}{\epsilon}\frac{q_{1}}{p_{1}},X_{2}=\frac{z_{0}}{p_{2}}+\frac{z_{0}}{\epsilon}\frac{q_{2}}{p_{2}},Y_{1}=\frac{z_{0}}{p_{1}}S_{\bar{N}},Y_{2}=\frac{z_{0}}{p_{2}}+\frac{z_{0}}{p_{1}}\frac{q_{2}}{p_{2}}S_{\bar{N}-1}$.
The sum $\frac{1}{L}\sum_{k=0}^{\bar{N}-1}S_{k}$ in Eq. (67), and consequently
$C^{[2]}_{0}$, can be obtained in closed form, given by Eq. (87) in Appendix
A. We also provide the successive expressions of $C^{[n]}_{0}$ for different
$n$, along with the corresponding formula for any general $n$ in Eq. (97)
(Appendix A). The variation of $C^{[n]}_{0}$ with the counter-flow parameter
$q_{1}$ for different values of $n$, is presented in Fig. 5. For each $n$ the
correlation, similar to density and current, displays the existence of two
different phases for $q_{1}<p_{1}$ and $q_{1}>p_{1}$. In Fig. 5, we observe a
steep increase in the correlations near $q_{1}\lesssim p_{1}$. Interestingly,
the correlations between consecutive vacancies increase considerably with
increasing $n$ for $q_{1}>p_{1}$. This implies more number of vacancies prefer
to stick together in the parameter regime $q_{1}>p_{1}$. It appears to be a
direct evidence of macroscopic cluster formation of vacancies in the
clustering phase. We would like to use the fact that $z_{0}$ remains almost
zero in the clustering phase (inset of Fig. 2), in Eq. (67). This simplifies
the expression for $C^{[2]}_{0}$ considerably in the clustering phase and
leads to

$C^{[2]}_{0}\approx\rho_{0}-\rho_{0}^{2},\hskip
11.38092pt\mathrm{for}\,\,q_{1}>p_{1}.$ (68)

The above result is consistent with $C^{[2]}_{0}$ being constant for
$q_{1}>p_{1}$ in Fig. 5, with the corresponding value agreeing to the same
given in Eq. (68). In fact, similar arguments apply for general $n$, and Eq.
(97) consequently results in the following particularly simple expression for
$C^{[n]}_{0}$ in the clustering phase

$C^{[n]}_{0}\approx\rho_{0}-\rho_{0}^{n},\hskip
11.38092pt\mathrm{for}\,\,q_{1}>p_{1}.$ (69)

We understand from Eq. (69) that with increasing $n$, $C^{[n]}_{0}$ approaches
to $\rho_{0}$ which is simply the density of the vacancies, and this statement
is also evident from Fig. 5.

##  VI Effect of non-ergodicity on clustering

The discussion of the two different phases up to now, corresponds to the
initial configuration Eq. (36). The features of the free flowing phase and the
clustering phase would remain the same for initial configurations that can be
prepared by permuting the positions of vacancies in Eq. (36). But, owing to
the non-ergodic nature of $2$-ASEP-IAF, it is important to ask about the
effect of the rearrangement of the two species and impurities in the initial
configuration Eq. (36). To answer this, in this section we would like to
investigate the effect of non-ergodicity on the clustering phenomena. We
consider the following variation of the initial configuration Eq. (36)

$D_{2}A..D_{2}A\,\,D_{1}D_{1}A..D_{1}D_{1}A\,\,D_{1}A..D_{1}A\,\,D_{1}..D_{1}\,\,E..E.$
(70)

Figure 6: The figure demonstrates the effect of rearrangement parameter $\eta$
(characterizing non-ergodicity) on the clustering phase. As $\eta$ decreases,
the distinction between the free flowing phase and clustering phase becomes
more and more abrupt. Particularly, the density remains monotonic for large
$\eta\approx 1$, whereas it is non-monotonic for small and moderate values of
$\eta$. Solid lines and dots correspond to analytical calculation and Monte
Carlo simulation results respectively. The parameters used are
$L=10^{3},p_{1}=0.4,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

We emphasize that all the input parameters
$(p_{1,2},\epsilon,w_{12,21},\rho_{0},\rho_{+},\bar{\rho})$ are same both for
Eqs. (36) and (71). The difference between these two initial configurations
lies in the rearrangement of the non-flipping species $1$ particles. In Eq.
(71), we have two types of non-flipping species $1$ particles, isolated (left
$D_{1}$ of unit $D_{1}D_{1}A$, cannot come in contact with another non-
flipping $D_{1}$), and non-isolated (belongs to sequence $D_{1}..D_{1}$). For
Eq. (36), we have only the non-isolated type. However, the total density of
non-flipping $D_{1}$ is $\bar{\rho}$, same for both Eq. (36) and Eq. (71). To
quantify their difference, we denote the fraction of isolated non-flipping
species $1$ particles in Eq. (71) as $\eta$. If we denote the densities of
isolated and non-isolated types of non-flipping species $1$ particles by
$\rho_{\mathrm{iso}}$ and $\rho_{\mathrm{niso}}$ respectively, then $\eta$ is
given by

$\eta=\frac{\rho_{\mathrm{iso}}}{\bar{\rho}}=1-\frac{\rho_{\mathrm{niso}}}{\bar{\rho}},$
(71)

where $\bar{\rho}$ is defined through Eq. (37). The initial configuration in
Eq. (36) corresponds to $\eta=0$. Since the variation of $\eta$ simply
rearranges the non-flipping $D_{1}$-s in the initial configuration, we denote
it as rearrangement parameter. Thus $\eta$ appears as a hallmark of non-
ergodicity. We present the behavior of average species density $\rho_{1}$ as a
function of $\eta$ in Fig. 6, obtained in Eq. (76) (Appendix A), following
similar methods discussed in Section IV. With decreasing $\eta,$ the
distinction between two phases become more evident in Fig. 6 and the fall of
the density to constant value gets sharper. Intriguingly, for small or
moderate values of rearrangement parameter, the density is non-monotonic,
contrary to its monotonic nature for $\eta\approx 1$. Figure 6 thus
illustrates that the onset of clustering strongly depends on $\eta$, and
therefore on the choice of the initial configuration. In reality, often there
are restrictions on the tunable range of the tuning parameter. For the
variable range of the parameter under consideration (e.g. say $q_{1}\in(0,1)$
for some system), Fig. 6 tells us which initial configurations are more prone
to clustering in the steady state and which initial configurations are
suitable to avoid such clustering, for a given value of $q_{1}$. Taking into
consideration the possible mapping of the present model to a narrow two-lane
system with counter-flow, the analysis of diagrams like Fig. 6 might help in
predicting the chances of jamming in the steady state starting from different
initial configurations.

##  VII Analysis of the case $q_{1}-p_{1}=p_{2}-q_{2}$: mapping to run-and-
tumble particles

Figure 7: The figures (a) and (b) show heat-maps for density $\rho_{1}$ and
current $J_{10}$ respectively, in the $p_{1}-q_{1}$ plane, for the special
case $q_{2}-p_{2}=p_{1}-q_{1}$ that maps the $2$-ASEP-IAF to a system of RTPs.
Both figures exhibit the existence of two different phases for $q_{1}>p_{1}$
and $q_{1}<p_{1}$. The current being close to zero for $q_{1}>p_{1}$ in figure
(b) confirms that this parameter region corresponds to the clustering phase,
while $q_{1}<p_{1}$ with non-zero current gives the free flowing phase. The
parameters used are
$L=10^{3},p_{2}=1.0,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

Figure 8: The figures (a) and (b) show two different phases for the RTPs,
through the variation of average species density $\rho_{1}$ and average drift
current $J_{10}$ respectively, as functions of the parameter $q_{1}$. For
$q_{1}>p_{1}$ and $q_{2}<p_{2}$, $\rho_{1}$ is constant and $J_{10}$ is
vanishingly small, characterizing the phase as clustering phase. For
$q_{1}<p_{1}$ and $q_{2}>p_{2}$, we observe the free flowing phase where both
density and current vary considerably with $q_{1}$. The solid lines and the
dots correspond to analytical and Monte Carlo simulation results,
respectively. The parameters used are
$L=10^{3},p_{1}=0.4,p_{2}=1.0,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

Figure 9: The figures (a) and (b) illustrate the effect of the tumbling rate
$w_{12}$ on the two phases. The figure (a) exhibits similar qualitative
behavior (monotonic decrease) of the density $\rho_{1}$ in both phases. The
figure (b) shows that the clustering phase remains unaffected with the change
of $w_{12}$, whereas the current decreases monotonically with $w_{12}$ in the
free flowing phase. The solid lines and the dots correspond to analytical and
Monte Carlo simulation results, respectively. The parameters used are
$L=10^{3},\epsilon=0.1,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

Figure 10: The figures (a) and (b) illustrate the effect of the drift
$\epsilon$ of tumbling reagents (i.e. impurities) on the two phases. Both the
density and current, in figures (a) and (b) respectively, show that the
variation of $\epsilon$ has no effect in the clustering phase, whereas both
the observables increase monotonically with $\epsilon$ in the free flowing
phase. The solid lines and the dots correspond to analytical and Monte Carlo
simulation results, respectively. The parameters used are
$L=10^{3},w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

In this section, we would like to consider the subspace of the parameter space
that satisfies the constraint $q_{1}-p_{1}=p_{2}-q_{2}$. This special case is
worthy of some detailed analysis because of its underlying connection to run-
and- tumble particles.

Self-propelled objects such as bacteria respond to chemical stimulus (e.g.
nutrients or harmful substances) present in the environment through the
process of chemotaxis. Chemotaxis, in a simple language, can be understood as
run-and-tumble motion. Run corresponds to motion along fixed direction and
tumble corresponds to intermediate random reorientation of direction of
motion. Individual swimming flagellated bacteria e.g. E. coli in suspensions
exhibit run-and-tumble dynamics Turner_2000 ; Berg_2004 ; Sidortsov_2017 ;
Beer_2019 . More interestingly, there are run-and-tumble wild bacteria e.g.
Vibrio ordalli in oceans, for which even the chemical stimulus can be dynamic
Brumley_2019 . The creation, destruction and diffusion of the short-lived
chemical stimulus in ocean create the dynamic chemical gradient Brumley_2019 .
Bacteria like E. coli or V. ordalli can accumulate near chemical stimulus to
create clusters Turner_2000 ; Berg_2004 ; Brumley_2019 . Such run-and-tumble
motion and clustering of bacteria can be best described in two dimensional
continuous space, which poses difficulties for exact analysis. Rather, one
dimensional discrete lattice models, although simpler adaptation of the
complex bacterial systems, can capture the run-and-tumble dynamics and cluster
formation. In this connection, it is intriguing to note that the $2$-ASEP-IAF
[Eq. (1)] can be mapped to a system of many RTPs, under the special
circumstance $q_{1}-p_{1}=p_{2}-q_{2}$.

In this case, the net bias of the two species are equal and opposite to each
other, and therefore the two species can be considered as the manifestations
of two possible orientations of RTPs in one dimension. The impurities act as
origins of the tumbling process. The flip between species can be interpreted
as tumbling of directions. Particularly, when the impurity causes flip to a
species that has net bias in the direction same as that of the impurities, the
impurity at that instant can be thought of a attracting reagent or nutrients.
On the other hand, when the resulting species from flip has net bias opposite
to the impurity bias, the impurity at that instant is considered as a
repelling reagent for tumbling process. Note that apart from the RTPs
(constituted by the $D_{1}$-s and $D_{2}$-s in the sequences of
$\mathcal{Y}=D_{1}A,D_{2}A$ in Eq. (36)) and sources of tumbling process (i.e.
impurities $A$ in the sequences of $\mathcal{Y}=D_{1}A,D_{2}A$ in Eq. (36)),
there are particles which do not tumble in the sense that their direction of
net bias remain unchanged (the non-flipping $D_{1}$-s in the sequence of
$\mathcal{Y}=D_{1}$ in Eq. (36)). With this mapping of the $2$-ASEP-IAF to
RTPs, we would like to discuss the behaviors of the average density [Eq. (56)]
and current [Eq. (62)] of the species for the special case
$q_{1}-p_{1}=p_{2}-q_{2}$, as a function of the system parameters.

We present the heat-maps of the average density $\rho_{1}$ and current
$J_{10}$ in Figs. 7(a) and (b), in the $p_{1}-q_{1}$ plane. It should be
mentioned that the mapping discussed in this section forces us to consider
$q_{2}-p_{2}=p_{1}-q_{1}$, i.e. Fig. 7 strictly corresponds to the counter-
flow situation (both for $q_{1}<p_{1}$ and $q_{1}>p_{1}$). This is contrary to
the setup in Fig. 1 where we have natural flow for $q_{1}<p_{1}$ and counter-
flow for $q_{1}>p_{1}$. In spite of this, we observe from Fig. 7 that, even
within the counter-flow scenario, two different phases emerge showing distinct
features for $q_{1}<p_{1}$ and $q_{1}>p_{1}$ that are qualitatively similar to
Fig. 1. For further analysis, in Figs. 8(a) and (b), we present cross sections
of Figs. 7(a) and (b), respectively, at fixed $p_{1}$. Both Figs. 8(a) and (b)
exhibit two different phases, the free flowing phase and the clustering phase.
For $q_{1}<p_{1}$ and $q_{2}>p_{2}$, both the density and current vary non-
monotonically with $q_{1}$, which we understand to be the free flowing phase.
Whereas, for $q_{1}>p_{1}$ and $q_{2}<p_{2}$, both $\rho_{1}$ and $J_{10}$
remain constant, specifically the current is vanishingly small identifying the
phase as the clustering phase. In comparison to Fig. 2 where $\rho_{1}$
attains maximum value in the free flowing phase, the maximum density for Fig.
8(a) is achieved in the clustering phase. The cluster formation owes to the
cumulative effect emerging from the presence of the non-tumbling sequence of
$\mathcal{Y}=D_{1}$ in the initial configuration Eq. (36). We should mention
that there is disagreement to some extent between our exact results and the
Monte Carlo simulation results near $q_{1}=p_{1}$ in the clustering phase
observed in Fig. 8(a), whereas both results match well for all other values of
$q_{1}$. This discrepancy, as per our current understanding, is due to the
restriction of insufficient ensemble average near the transition point.

It is quite natural to ask the effect of tumbling rates ($w_{12}$ and
$w_{21}$) on the two phases. To investigate this, we present the variation of
$\rho_{1}$ and $J_{10}$ as functions of the tumbling rate $w_{12}$ in Figs.
9(a) and (b), respectively, using our analytical findings (solid lines) from
Eqs. (56) and (62), supported by Monte Carlo simulation results (dots). In
each figure, we consider two situations, one corresponding to $q_{1}<p_{1}$
and $q_{2}>p_{2}$, and the other to $q_{1}>p_{1}$ and $q_{2}<p_{2}$.
Particularly, Fig. 9(b) clearly shows the current remains vanishingly small
for any value of $w_{12}$, in the case $q_{1}>p_{1}$ and $q_{2}<p_{2}$,
suggesting that the formed cluster is stable to change of the tumbling rates.
On the other hand, for $q_{1}<p_{1}$ and $q_{2}>p_{2}$, the current of species
$1$ decreases monotonically with increasing tumbling rate $w_{12}$ (species
$1$ to species $2$) as expected. However, the fact that the change of tumbling
rate does not have any effect on the formed clusters, is not visible from the
density characteristics in Fig. 9(a) which shows similar qualitative nature
for both the free flowing and clustering phase. We would also like to explore
the effect of the drift process of the impurities, that act as attracting
reagent (nutrient) or repelling reagent causing the tumbles, on the two
phases. The observables $\rho_{1}$ and $J_{10}$ in Figs. 10(a) and (b),
respectively, exhibit that the clustering phase is undisturbed by the change
of the impurity drift rate $\epsilon$, where as both the average density and
current increases with increasing $\epsilon$ in the free flowing phase. Thus
we have shown that, starting from initial configuration Eq. (36), the system
of RTPs (along with attracting or repelling reagents causing tumbles) can
transit between free flowing phase and clustering phase, with the variation of
the control parameter $q_{1}$ that tunes the net biases of the two tumbling
directions. Interestingly, the formed clusters are stable against any change
of the tumbling rates $w_{12}$ (and $w_{21}$) and reagent (i.e. impurity
causing tumble) drift rate $\epsilon$.

##  VIII Summary

We have demonstrated the formation of clusters induced by counter-flow in a
non-ergodic system. To illustrate this, we have considered two species
asymmetric simple exclusion process along with impurities. Apart from the
drift of the species and impurities, additionally the impurities activate
flips between the two species. The exact analytical results obtained for the
observables (average species density, drift current, spatial correlation),
show two distinct phases, free flowing phase and clustering phase, as
functions of the counter-flow parameter $q_{1}$, which basically controls the
flow direction of species $1$ by tuning its net bias. The free flowing phase
is specified by non-monotonic density and non-monotonic finite current,
whereas the clustering phase is characterized by constant density and
vanishing current. The clustering situation can be thought as an equivalent
system almost devoid of vacancies, where the drift dynamics is almost absent.
This is compatible with the vanishing current and the constant density in the
clustering phase. Specifically, the growth of $n$-point spatial correlation
between $n$ consecutive vacancies with increasing $n$, directly implies the
accumulation of vacancies to form one macroscopic cluster, along with another
macroscopic cluster formed by all the particles. The heat maps of density and
current in $p_{1}-q_{1}$ plane show that the two phases are demarcated by the
line $q_{1}=p_{1}$ ($p_{1}$ and $q_{1}$ being right and left hop rates of
species $1$, respectively). This is further clarified by the sharp descents of
density and current and the sharp ascent in spatial correlation between
vacancies near $q_{1}=p_{1}$. The slow drift of the cluster in the clustering
phase can be observed by tracking individual configurations over long time.
The effect of non-ergodicity on the system is characterized through a
rearrangement parameter, which only permutes the positions of some particles
while keeping all the input parameters fixed. We observe that the choice of
initial configuration affects the onset of clustering significantly. In fact,
for certain initial configurations, we do not see signatures of clustering for
finite values of $q_{1}$ and the corresponding densities increase
monotonically in contrast to the non-monotonic densities for initial
configurations showing cluster formations. Interestingly, for a special case
of the microscopic dynamics when the net bias of the two species are equal and
opposite to each other, we map the $2$-ASEP-IAF to a system of RTPs. The
species with net bias to right (left) can be considered as right (left)
running RTPs, where the impurities act as tumbling reagents that cause the
tumbling of RTPs, equivalent to the flip process in $2$-ASEP-IAF. Notably,
although this mapping is valid only in the counter-flow situation, we still
observe two different phases, the free flowing phase and the clustering phase.
We further find that the clustering phase remains stable against the variation
of the tumbling rate and drift of the tumbling reagents. We believe that our
analysis supported by exact analytical results enlightens the understanding of
clustering phenomena.

The model studied here, having resemblance to two-lane traffics and RTPs in
active matter, points towards analytical understanding of traffic jams and
clustering active matter systems. Further careful and thorough investigations
are required to establish such connections. It would be interesting to study
variations of the local microscopic dynamics considered here, that can produce
dynamical ways to get rid of the clustering phase. It would be important to
explore the effects of boundary conditions, e.g. allowing entries and exits of
selective particles or all particles, on the different phases of the system.

Acknowledgements.- This work is partially supported by the Grants-in-Aid for
Scientific Research (JSPS KAKENHI Grant No. JP21H01006). A.K.C. gratefully
acknowledges postdoctoral fellowship from the YITP. The numerical calculations
have been done on the cluster Yukawa-21 at the YITP.

##  Appendix A Derivation of observables

In this appendix, we sketch the steps for calculating the average species
densities and $n$-point correlation function between consecutive vacancies.
These steps essentially follow the methods discussed in obtaining the
partition function in Section IV. The average density $\rho_{1}$ of species
$1$, can be written as

$\displaystyle\rho_{1}=\langle 1\rangle$ $\displaystyle=$
$\displaystyle\bar{\rho}+\frac{\left(\rho_{+}-\frac{1}{L}\right)}{Q}\sum_{m_{1}=0}^{\infty}\dots\sum_{n_{\bar{N}}=0}^{\infty}\,\,\mathrm{Tr}\left[D_{1}(z_{0}E)^{m_{1}}A(z_{0}E)^{\bar{m}_{1}}\,\prod_{i=2}^{N_{+}}(D_{1}+D_{2})(z_{0}E)^{m_{i}}A(z_{0}E)^{\bar{m}_{i}}\,\,\prod_{j=1}^{\bar{N}}D_{1}(z_{0}E)^{n_{j}}\right]$
$\displaystyle+\frac{1}{L}\frac{1}{Q}\sum_{m_{1}=0}^{\infty}\dots\sum_{n_{\bar{N}}=0}^{\infty}\,\,\mathrm{Tr}\left[\prod_{i=1}^{N_{+}-1}(D_{1}+D_{2})(z_{0}E)^{m_{i}}A(z_{0}E)^{\bar{m}_{i}}\,\,D_{1}(z_{0}E)^{m_{N_{+}}}A(z_{0}E)^{\bar{m}_{N_{+}}}\,\,\prod_{j=1}^{\bar{N}}D_{1}(z_{0}E)^{n_{j}}\right].$

The term $\bar{\rho}$ appears directly due to the initial density $\bar{\rho}$
of the non-flipping species $1$ particles. In the second part with pre-factor
$\left(\rho_{+}-\frac{1}{L}\right)$, we place at least one $D_{1}$ in a
flipping term whereas any other flipping term can have $D_{1}$ or $D_{2}$, the
density of such terms is $\left(\rho_{+}-\frac{1}{L}\right)$. The last part
contributes due to the $D_{1}$ that is placed in the last flipping term after
which the non-flipping $D_{1}$-s start. To proceed, we would use the explicit
matrix representations in Eq. (35) and the expressions derived in Eqs. (45)
and (49) to evaluate the matrix strings in Eq. (LABEL:eq:ob1). This leads us
to the final formula for the average species density in Eq. (56). We can
follow similar procedure to calculate the average species densities for the
general initial configuration Eq. (71) with nonzero $\eta$ in Section VI. In
this case, we get,

$\displaystyle\rho_{1}(\eta)$ $\displaystyle=$
$\displaystyle\bar{\rho}+\left(\rho_{+}-\frac{1}{L}-\eta\bar{\rho}\right)\frac{\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{0}}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{0}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}\left(1+\frac{q_{2}}{p_{1}}S_{-1}\right)}\right]}+\eta\bar{\rho}\frac{\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{1}}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{1}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}\left(1+\frac{q_{2}}{p_{1}}S_{0}\right)}\right]}$
(75)
$\displaystyle+\frac{1}{L}\frac{\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{(1-\eta)\bar{N}}}}{\left[\frac{w_{21}}{1-\frac{z_{0}}{p_{1}}S_{(1-\eta)\bar{N}}}+\frac{w_{12}}{1-\frac{z_{0}}{p_{2}}\left(1+\frac{q_{2}}{p_{1}}S_{(1-\eta)\bar{N}-1}\right)}\right]},$
$\displaystyle\rho_{2}(\eta)$ $\displaystyle=$
$\displaystyle\bar{\rho}+\rho_{+}-\rho_{1}(\eta),$ (76)

It is straightforward to check that $\eta=0$ in Eq. (76) gives the density in
Eq. (56).

The formal expression for two-point nearest neighbor correlation between
vacancies is

$C^{[2]}_{0}=\langle 00\rangle-\rho_{0}^{2}.$ (77)

It is difficult to calculate $\langle 00\rangle$ directly using the matrix
representations, rather it is easier to use the following conservation

$\langle 00\rangle+\langle 01\rangle+\langle 02\rangle+\langle
0+\rangle=\rho_{0}.$ (78)

Using Eq. (78) into Eq. (77), we get

$C^{[2]}_{0}=\rho_{0}-\rho_{0}^{2}-\langle 01\rangle-\langle 02\rangle-\langle
0+\rangle.$ (79)

Evaluating $\langle 01\rangle$, $\langle 02\rangle$ and $\langle 0+\rangle$
using the matrix representations, we finally obtain the following expression
for $C_{00}$ as

$\displaystyle
C^{[2]}_{0}=\rho_{0}-\rho_{0}^{2}-\left(\rho_{+}-\frac{1}{L}\right)\frac{z_{0}}{\epsilon}-\frac{1}{L}\frac{z_{0}}{p_{1}}\sum_{k=0}^{\bar{N}-1}S_{k}$
(80)
$\displaystyle-\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}X_{1}}{1-X_{1}}+\frac{w_{12}X_{2}}{1-X_{2}}}{\frac{w_{21}}{1-X_{1}}+\frac{w_{12}}{1-X_{2}}}-\frac{1}{L}\frac{\frac{w_{21}Y_{1}}{1-Y_{1}}+\frac{w_{12}Y_{2}}{1-Y_{2}}}{\frac{w_{21}}{1-Y_{1}}+\frac{w_{12}}{1-Y_{2}}},$
(81)

where

$\displaystyle X_{1}$ $\displaystyle=$
$\displaystyle\frac{z_{0}}{p_{1}}+\frac{z_{0}}{\epsilon}\frac{q_{1}}{p_{1}},$
(82) $\displaystyle X_{2}$ $\displaystyle=$
$\displaystyle\frac{z_{0}}{p_{2}}+\frac{z_{0}}{\epsilon}\frac{q_{2}}{p_{2}},$
(83) $\displaystyle Y_{1}$ $\displaystyle=$
$\displaystyle\frac{z_{0}}{p_{1}}S_{\bar{N}},$ (84) $\displaystyle Y_{2}$
$\displaystyle=$
$\displaystyle\frac{z_{0}}{p_{2}}+\frac{z_{0}}{p_{1}}\frac{q_{2}}{p_{2}}S_{\bar{N}-1}.$
(85)

Clearly, the two-point correlation would have a closed form if the sum
$\sum_{k=0}^{\bar{N}-1}S_{k}$ has a closed form, and indeed this can be
evaluated as

$\displaystyle\frac{1}{L}\sum_{k=0}^{\bar{N}-1}S_{k}=\bar{\rho}\frac{p_{1}}{(p_{1}-q_{1})}$
(86)
$\displaystyle+\frac{1}{L}\left[1-\left(\frac{q_{1}}{p_{1}}\right)^{\bar{N}}\right]\left(\frac{p_{1}-q_{1}}{\epsilon}-1\right)\frac{p_{1}q_{1}}{(p_{1}-q_{1})^{2}}.$
(87)

So, we have evaluated the two-point correlation function $C^{[2]}_{0}$
exactly. Of course, $z_{0}$ has to be calculated from the density-fugacity
relation. Importantly in Eq. (87), note that we could scale $\bar{N}$ by
system-size $L$ properly in the first term so that it becomes a function of
$\bar{\rho}$, but this is not possible in case of the second term where
$\bar{N}$ appears in the power.

In fact, using the result of $C^{[2]}_{0}$, we can calculate $C^{[3]}_{0}$ and
then $C^{[4]}_{0}$ using $C^{[3]}_{0}$, and so on. For simplified notations,
we denote $C_{0}^{[n]}$ as the correlation between consecutive $n$ vacancies
[Eq. (65)], then we obtain in iterative way,

$\displaystyle C_{0}^{[2]}$ $\displaystyle=$
$\displaystyle\rho_{0}-\rho_{0}^{2}-\left(\rho_{+}-\frac{1}{L}\right)\frac{z_{0}}{\epsilon}-\frac{1}{L}\frac{z_{0}}{p_{1}}\sum_{k=0}^{\bar{N}-1}S_{k}-\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}X_{1}}{1-X_{1}}+\frac{w_{12}X_{2}}{1-X_{2}}}{\frac{w_{21}}{1-X_{1}}+\frac{w_{12}}{1-X_{2}}}-\frac{1}{L}\frac{\frac{w_{21}Y_{1}}{1-Y_{1}}+\frac{w_{12}Y_{2}}{1-Y_{2}}}{\frac{w_{21}}{1-Y_{1}}+\frac{w_{12}}{1-Y_{2}}},$
(88) $\displaystyle C_{0}^{[3]}$ $\displaystyle=$
$\displaystyle{C_{0}^{[2]}+\rho_{0}^{2}-\rho_{0}^{3}}-\left(\rho_{+}-\frac{1}{L}\right)\left(\frac{z_{0}}{\epsilon}\right)^{2}-\frac{1}{L}\left(\frac{z_{0}}{p_{1}}\right)^{2}\sum_{k=0}^{\bar{N}-1}\left(S_{k}\right)^{2}-\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}{X_{1}^{2}}}{1-X_{1}}+\frac{w_{12}{X_{2}^{2}}}{1-X_{2}}}{\frac{w_{21}}{1-X_{1}}+\frac{w_{12}}{1-X_{2}}}-\frac{1}{L}\frac{\frac{w_{21}{Y_{1}^{2}}}{1-Y_{1}}+\frac{w_{12}{Y_{2}^{2}}}{1-Y_{2}}}{\frac{w_{21}}{1-Y_{1}}+\frac{w_{12}}{1-Y_{2}}},$
(90) $\displaystyle C_{0}^{[4]}$ $\displaystyle=$
$\displaystyle{C_{0}^{[3]}+\rho_{0}^{3}-\rho_{0}^{4}}-\left(\rho_{+}-\frac{1}{L}\right)\left(\frac{z_{0}}{\epsilon}\right)^{3}-\frac{1}{L}\left(\frac{z_{0}}{p_{1}}\right)^{3}\sum_{k=0}^{\bar{N}-1}\left(S_{k}\right)^{3}-\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}{X_{1}^{3}}}{1-X_{1}}+\frac{w_{12}{X_{2}^{3}}}{1-X_{2}}}{\frac{w_{21}}{1-X_{1}}+\frac{w_{12}}{1-X_{2}}}-\frac{1}{L}\frac{\frac{w_{21}{Y_{1}^{3}}}{1-Y_{1}}+\frac{w_{12}{Y_{2}^{3}}}{1-Y_{2}}}{\frac{w_{21}}{1-Y_{1}}+\frac{w_{12}}{1-Y_{2}}},$
(92) $\displaystyle\dots$ $\displaystyle=$ $\displaystyle\dots$ (94)
$\displaystyle C_{0}^{[n]}$ $\displaystyle=$
$\displaystyle{C_{0}^{[n-1]}+\rho_{0}^{n-1}-\rho_{0}^{n}}-\left(\rho_{+}-\frac{1}{L}\right)\left(\frac{z_{0}}{\epsilon}\right)^{n-1}-\frac{1}{L}\left(\frac{z_{0}}{p_{1}}\right)^{n-1}\sum_{k=0}^{\bar{N}-1}\left(S_{k}\right)^{n-1}$
(97)
$\displaystyle-\left(\rho_{+}-\frac{1}{L}\right)\frac{\frac{w_{21}{X_{1}^{n-1}}}{1-X_{1}}+\frac{w_{12}{X_{2}^{n-1}}}{1-X_{2}}}{\frac{w_{21}}{1-X_{1}}+\frac{w_{12}}{1-X_{2}}}-\frac{1}{L}\frac{\frac{w_{21}{Y_{1}^{n-1}}}{1-Y_{1}}+\frac{w_{12}{Y_{2}^{n-1}}}{1-Y_{2}}}{\frac{w_{21}}{1-Y_{1}}+\frac{w_{12}}{1-Y_{2}}}.$

Thus we have obtained the analytical formulae for average species densities,
drift currents and $n$-point correlation between consecutive vacancies. The
results in Eq. (97) have been used to present the behavior of the correlations
with the variation of the counter flow parameter $q_{1}$ in Fig. 5.

##  Appendix B Density-fugacity relation: solution for special cases

Here we would like to state the explicit form of the density-fugacity
relation, calculated from the partition function in Eq. (51). This relation is
used to solve the fugacity $z_{0}$ as a function of the input parameters.
Consequently, we can replace the corresponding value of $z_{0}$ in the
expressions of the observables e.g. in Eqs. (56), (62) and (67) so that they
become functions of the input parameters only. The formal expression for the
density-fugacity relation is given in Eq. (54). Using Eq. (51) in Eq. (54), we
have the following explicit form of the density-fugacity relation to solve,

$\displaystyle\frac{\rho_{+}}{1-\frac{z_{0}}{\epsilon}}+\frac{z_{0}\left(\rho_{+}-\frac{1}{L}\right)}{w_{21}(1-z_{0}X^{\prime}_{2})+w_{12}(1-z_{0}X^{\prime}_{1})}\left[w_{21}X^{\prime}_{1}\frac{1-z_{0}X^{\prime}_{2}}{1-z_{0}X^{\prime}_{1}}+w_{12}X^{\prime}_{2}\frac{1-z_{0}X^{\prime}_{1}}{1-z_{0}X^{\prime}_{2}}\right]$
(98)
$\displaystyle+\frac{\frac{z_{0}}{L}}{w_{21}(1-z_{0}Y^{\prime}_{2})+w_{12}(1-z_{0}Y^{\prime}_{1})}\left[w_{21}Y^{\prime}_{1}\frac{1-z_{0}Y^{\prime}_{2}}{1-z_{0}Y^{\prime}_{1}}+w_{12}Y^{\prime}_{2}\frac{1-z_{0}Y^{\prime}_{1}}{1-z_{0}Y^{\prime}_{2}}\right]+\frac{1}{L}\sum_{k=1}^{\bar{N}}\frac{1}{1-\frac{z_{0}}{p_{1}}S_{k-1}}=\rho_{0}+\rho_{+}+\bar{\rho},$
(99)

with $X^{\prime}_{1,2}=X_{1,2}/z_{0}$ and $Y^{\prime}_{1,2}=Y_{1,2}/z_{0}$,
where $X_{1,2}$ and $Y_{1,2}$ follow Eq. (85), $S_{k}$ is given by Eq. (53).
The reason behind such rescaling by $z_{0}$, is simply to express
$X^{\prime}_{1,2}$ and $Y^{\prime}_{1,2}$ as functions of the input parameters
only. In general, we have solved Eq. (99) in Mathematica to get $z_{0}$ for a
given set of input parameters. Also note that the complexity of the equation
increases with increasing number of non-flipping species $1$ particles
$\bar{N}$, in terms of the highest degree of $z$ present in the polynomial of
$z_{0}$ in Eq. (99). For some special cases with specific choices of the hop-
rates, one can obtain closed form solutions for the fugacity $z_{0}$.

A particularly simple case corresponds to $\bar{N}=0$. This also implies
$Y^{\prime}_{1}=X^{\prime}_{1}$ and $Y^{\prime}_{2}=X^{\prime}_{2}$, which can
be seen directly from Eq. (85) with the help of Eq. (53). We further consider
the special situation of $X^{\prime}_{1}=X^{\prime}_{2}$. The density-fugacity
relation Eq. (99) simplifies to

$\frac{z_{0}\rho_{+}X^{\prime}_{1}}{1-z_{0}X^{\prime}_{1}}+\frac{\rho_{+}}{1-\frac{z_{0}}{\epsilon}}=\rho_{0}+\rho_{+}.$
(100)

The above equation has the following solution

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

To better understand the constraint on the hop rates for which we have got the
exact solution Eq. (101), we explore the situation
$X^{\prime}_{1}=X^{\prime}_{2}$, which basically boils down to

$\epsilon=\frac{q_{2}p_{1}-q_{1}p_{2}}{p_{2}-p_{1}}.$ (102)

The above subspace of hop rates can create both natural flow and counter-flow
situations and also includes the very special case $p_{1}=p_{2}$ and
$q_{1}=q_{2}$. We should mention that, even without the assumption
$X^{\prime}_{1}=X^{\prime}_{2}$, we have a quartic equation in $z_{0}$ that
can be solved exactly in Mathematica. However, the solution of $z_{0}$ in that
case is too lengthy to include here.

Another noteworthy point is, the fugacity $z_{0}$ actually equals to the
cluster velocity $v_{\mathrm{cl}}$ discussed in the main text (see inset of
Fig. 4). This is evident from the current-density relation Eq. (62), which can
be considered as $J=v_{\mathrm{cl}}\rho$ in the clustering phase. In fact, one
can check the inset of Fig. 4 in the main text, obtained from Monte Carlo
simulations, can be reproduced by calculating $z_{0}$ for the corresponding
set of input parameters.

##  Appendix C A comment regarding the initial configuration

We have considered step-like initial configuration [Eq. (36)] in the main
text. To elaborate, initially all the particles (both species and impurities)
occupy consecutive lattice sites with no vacancy between them. Starting from
such step-like initial configuration, in the free flowing phase, the vacancies
get randomly distributed between the particles. On the other hand, in the
clustering phase, any steady state configuration remains step-like, with the
particle cluster shifting slowly to right.

Naturally, the question arises if the cluster can be formed from an initial
configuration which is not step-like, rather there are vacancies distributed
between particles. The answer is yes. If we start from a non-step-like initial
configuration given below,

$C(0)=D_{2}E^{m_{1}}AE^{\bar{m}_{1}}\dots
D_{2}E^{m_{N_{+}/2}}AE^{\bar{m}_{N_{+}/2}}\,\,D_{1}E^{n_{1}}AE^{\bar{n}_{1}}\dots
D_{1}E^{n_{N_{+}/2}}AE^{\bar{n}_{N_{+}/2}}\,\,D_{1}E^{r_{1}}\dots
D_{1}E^{r_{\bar{N}}},$ (103)

where the total number of vacancies is

$\sum_{i=1}^{N_{+}/2}\left(m_{i}+\bar{m}_{i}+n_{i}+\bar{n}_{i}\right)+\sum_{j=1}^{\bar{N}}r_{k}=N_{0}.$
(104)

Figure 11: The figure demonstrates the step-like cluster formation in the
steady state, starting from a non-step-like initial configuration where
vacancies are randomly distributed among particles. Here we consider $\tau=0$
as the observation time of cluster in steady state, whereas the initial
configuration corresponds to time $t=0$. The parameters used are
$L=10^{3},p_{1}=0.4,q_{1}=0.5,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.29,\rho_{0}=0.4$.

Note that, the above initial configuration contains exact same ordering of
particles just like the one [Eq. (36)] studied in the main text and the total
number of vacancies are also the same for both configurations, only difference
being the initial configuration [Eq. (36)] in the main text is step-like,
whereas Eq. (103) is non-step-like. Since the ordering of vacancies actually
do not matter, both of these initial configurations lead to the same
configuration sub-space in the steady state. Thereby all the characteristics
of the system in the steady state remain same for both of these initial
configurations. Thus we expect to see the clustering phenomena starting from
initial configuration Eq. (103) just like we did for the one in the main text
(see Fig. 4). Indeed, in Fig. 11, we observe that the macroscopic cluster is
formed in the steady state, while the initial configuration is non-step-like.
So, the only advantageous and satisfactory thing about Eq. (103), is the fact
that the cluster is formed in a dynamic way from non-step-like initial
configuration.

##  Appendix D Variation of observables with system size $L$

In this appendix, we present the behaviors of average species density and
drift current as function of $q_{1}$, for different system sizes $L$.

Figure 12: The figure (a) illustrates the convergence of the fugacity behavior
with increasing system size $L$, plotted against the parameter $q_{1}$. The
figure (b) how the approach of $z_{0}$ towards this convergence with varying
$L$, at the transition point $q_{1}=p_{1}$. The exact numerical data obtained
from the solution of $z_{0}$ is fitted with the form $z_{0}=cL^{\alpha}$,
resulting in $\alpha\approx-0.76$. The parameters used are
$p_{1}=0.4,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

Figure 13: The figures (a) and (b) show convergence of $\rho_{1}$ and $J_{10}$
respectively, with increasing system size $L$. The insets (in log scale) of
(a) and (b) show the variations of $\rho_{1}$ and $J_{10}$ respectively with
$L$ at the transition point $q_{1}=p_{1}$. They follow the relations
$\rho_{1}=\bar{\rho}c^{\prime\prime}L^{\delta}$ and
$J_{10}=c^{\prime}L^{\gamma}$ with exponent values $\delta\approx-0.078$ and
$\gamma\approx-0.88$ obtained from fitting. The parameters used are
$p_{1}=0.4,p_{2}=1.0,q_{2}=0.5,\epsilon=0.1,w_{12}=1.0,w_{21}=0.1,\rho_{+}=0.25,\rho_{0}=0.4$.

.

Note that, the density $\rho_{1}$ should be an intensive quantity independent
of $L$ and the steady state current should also be a constant independent of
$L$. In fact, like the fugacities e.g. temperature, chemical potential which
are independent of system size, the fugacity $z_{0}$ should also be a constant
independent of system size. Thus, $z_{0}$, $\rho_{1}$ and $J_{10}$ here must
show convergence with increasing $L$. Indeed, the fugacity obtained from the
exact numerical data by solving Eq. (99), show converging behavior with
increasing system size in Fig. 12(a). At the transition point $q_{1}=p_{1}$,
we investigate the approach of $z_{0}$ towards convergence as a function of
$L$, by fitting the data with the functional form $z_{0}=cL^{\alpha}$ in Fig.
12(b). Also, both Figs. 13(a) and (b) exhibit convergence of the corresponding
behaviors of $\rho_{1}$ and $J_{10}$ with increasing system size. In Fig.
13(b), the current shows rapid convergence in both the free flowing phase and
clustering phase for $L=10^{3}$ and more. Notably, in Fig. 13(a), the density
shows faster convergence with increasing $L$ in the free flowing phase, in
comparison to that of the clustering phase. We should mention that the
convergence value of $\rho_{1}$ in the clustering phase in Fig. 13(a), agrees
with the simplified formula obtained in Eq. (57). In the insets of Figs. 13(a)
and (b), we examine the approach of $\rho_{1}$ and $J_{10}$ respectively
towards convergence, with varying system size $L$. We fit the corresponding
data with the relations of the form $\rho_{1}=c^{\prime\prime}L^{\delta}$ and
$J_{10}=c^{\prime}L^{\gamma}$. Keeping in mind the relation [Eq. (62)] between
current, density and fugacity of the form $J_{10}=z_{0}\rho_{1}$, the exponent
$\delta$ can be determined from the relation $\delta=(\gamma-\alpha)$. Our
findings in Figs. 12 and 13(b) predict $\delta$ to be $(\gamma-\alpha)=-0.12$.
This is close to but not exactly equal to the $\delta$ value ($\approx-0.078$)
obtained in the inset of Fig. 13(a), because of the finite system sizes used
in the discussion.

##  Appendix E Arndt-Heinzel-Rittenberg model of counter-flow: comparisons

In this section, we compare the microscopic dynamics of our model with the
dynamics of the Arndt-Heinzel-Rittenberg (AHR) model that is known to exhibit
three different phases Arndt_1998 ; sArndt_1998_1 ; sArndt_1999 . The AHR
model considers positive (say, species $1$) and negative (say, species $2$)
particles along with vacancies ($0$) on a one dimensional periodic lattice and
the follow the dynamical rules given below

$\displaystyle
10\,\,\stackrel{{\scriptstyle\lambda}}{{\longrightarrow}}\,\,01,\hskip
22.76228pt02\,\,\stackrel{{\scriptstyle\lambda}}{{\longrightarrow}}\,\,20,\hskip
22.76228pt12\,\,\xrightleftharpoons[1]{q}\,\,21.$ (105)

A straightforward comparison of the AHR model in Eq. (105) with our model in
Eq. (1) reveal the following factors: (i) our model allows a non-conserving
flip dynamics that activates transformations between species $1$ and $2$,
which is absent in the AHR model for which each microscopic dynamics maintains
particle number conservation of every species. Also, as a consequence, our
model requires minimum four species in total (species $1$, species $2$,
impurity and vacancy) to operate both drift and flip dynamics, while the AHR
model deals with three species in total (species $1$, species $2$, and
vacancy). (ii) AHR model has an exchange dynamics that allows the two species
to exchange their positions. Such exchange dynamics is absent in our case.
(iii) Our model is non-ergodic in contrast to the ergodic nature of the AHR
model.

The exchange rate $q$ is considered as the tuning parameter for the AHR model
and the density $\rho$ of the two species are taken to be equal Arndt_1998 ;
sArndt_1999 . Three different phases, namely pure phase, mixed phase and
disordered phase are observed, as $q$ is varied Arndt_1998 ; sArndt_1999 .
When $q<1$, the species $1$ is more probable to reside at left of species $2$
which again likes to be at left of vacancies. This leads to a complete species
segregation with three types of blocks each purely consisting of one species
(either $1$ or $2$ or $0$), thereby referred as pure phase. For $1<q<q_{c}$
(where $q_{c}$ depends on $\lambda$ and $\rho$), a condensate is formed that
has both species $1$ and species $2$ mixed up, accompanied by a fluid
consisting of vacancies and some particles of the two species. This phase is
known as the mixed phase. There is no species segregation or condensate
formation for $q>q_{c}$, which is the disordered phase. To observe these three
phases, two point functions like drift current and correlations between
different species, have been used Arndt_1998 ; sArndt_1999 .

Interestingly, in our model, due to the presence of non-conserving flip
dynamics, we have even simpler one point function like average species
densities among observables of interest, and indeed the average species
densities clearly show the existence of two different phase, the free flowing
phase and the clustering phase. The free flowing phase in our model is similar
to the disordered phase of AHR model. On the other hand, for the specific
choice of initial configuration considered here and due to the flip dynamics,
species segregation in pure form is not possible in the clustering phase.
Rather, we have two macroscopic clusters, one consisting of only vacancies and
the other consisting of all kinds of particles (species $1$ and species $2$
and impurities). Although the mixing up of different species and impurities
inside the particle cluster has resemblance to the mixed phase condensate of
AHR model, we do not have a background fluid in our case. Rather, for any
$q_{1}>p_{1}$, we have only two clusters in the clustering phase, a vacancy
cluster and a particle cluster, with the particle cluster drifting along right
with considerably small velocity that depends on the density of the non-
flipping species $1$ particles in the system. Another noteworthy point in our
analysis is the rearrangement parameter whose variation captures the effect of
non-ergodicity on the clustering phenomenon, there is no such counterpart in
the ergodic AHR model.

Note that exact analysis in Refs. Rajewsky_2000  and Sasamoto_2001  later
revealed that there is actually no phase transition between mixed and
disordered phase in the AHR model, in the thermodynamic limit within the grand
canonical ensemble framework. This conclusion is associated with the existence
of extremely long but still finite correlation lengths in the infinite system.

There is an alternative approach to compare the dynamics of the two models,
although the key points of the comparative analysis between them remain the
same. Our tuning parameter has been the counter-flow parameter $q_{1}$, which
is a part of the drift dynamics. To treat the tuning parameter $q$ of the AHR
model on an equivalent footing, one can relabel species $1$, species $2$ and
vacancy in AHR model as vacancy, species $1$ and species $+$ respectively.
Thus one arrives at an alternative version of the AHR model given below

$\displaystyle
0+\,\,\stackrel{{\scriptstyle\lambda}}{{\longrightarrow}}\,\,+0,\hskip
22.76228pt+1\,\,\stackrel{{\scriptstyle\lambda}}{{\longrightarrow}}\,\,1+,\hskip
22.76228pt10\,\,\xrightleftharpoons[q]{1}\,\,01.$ (106)

Here we see, in comparison to our model Eq. (1), the species $2$ is absent and
we do not term $+$ as impurity because there is no flip dynamics at all. So
$+$ drifts to left only, species $1$ drifts to right or left with rates $q$
and $1$ respectively, with an additional exchange of positions between $+$ and
$1$. Clearly $q>1$ here corresponds to natural flow situation and $q<1$ refers
to the counter-flow situation.

The exact steady state of AHR model in Eq. (105) has been obtained in matrix
product from, which has a two dimensional representation in the limit
$q\rightarrow\infty$ and infinite dimensional representations in general
sArndt_1999 , which have different structures in comparison to the infinite
dimensional matrices in our case Eq. (35).

## References

  * (1) A. Liu, and S. R. Nagel, Nature 396, 21 (1998).
  * (2) G. Biroli, Nature Physics 3, 222 (2007).
  * (3) R. P. Behringer, and B. Chakraborty,
Rep. Prog. Phys. 82, 012601 (2019).

  * (4) K. Nagel, and M. Schreckenberg,
J. Phys. I France 2, 2221 (1992).

  * (5) D. Helbing, Rev. Mod. Phys. 73, 1067 (2001).
  * (6) T. Nagatani, Rep. Prog. Phys. 65, 1331 (2002).
  * (7) A. Schadschneider, D. Chowdury and K. Nishinari, Stochastic transport in complex system: From molecules to vehicles (Elsevier, Amstrdam, 2011).
  * (8) R. Mandal, P. J. Bhuyan, P. Chaudhuri, C. Dasgupta, and M. Rao, Nat Commun. 11, 2581 (2020).
  * (9) J. Yang, R. NI, and J. P. Ciamarra,
Phys. Rev. E 106, L012601 (2022).

  * (10) M. Sadati, N. T. Quazvini, R. Krishnan, C. Y. Park, and J. J. Fredberg, Differentiation 86, 121 (2013).
  * (11) L. Oswald, S. Grosser, D. M. Smith, and J. A. Käs,
J. Phys. D: Appl. Phys. 50, 483001 (2017).

  * (12) V. Trappe, V. Prasad, L. Cipelletti, P. N. Segre, and D. A. Weitz, Nature 411, 772 (2001).
  * (13) M. Muramatsu, T. Irie, and T. Nagatani,
Physica A 267, 487 (1999).

  * (14) Y. Tajima, K. Takimoto, and T. Nagatani,
Physica A 313, 709 (2002).

  * (15) D. Helbing, and P. Molnár, Phys. Rev. E 51, 4282 (1995).
  * (16) S.P. Hoogendoorn, S. Luding, P.H.L. Bovy, M. Schreckenberg, and D.E. Wolf (Eds.), Traffic and Granular Flow ’03, Springer, Berlin, Heidelberg, New York, (2005).
  * (17) T. Kretz, M. Wölki, and M. Schreckenberg,
J. Stat. Mech. (2006) P02005.

  * (18) W.W.F. Klingsch, C. Rogsch, A. Schadschneider, and M. Schreckenberg (Eds.), Pedestrian and Evacuation Dynamics 2008, Springer, Heidelberg, Dordrecht, London, New York, 2010.
  * (19) T. Kretz, A. Grünebohm, M. Kaufman, F. Mazur, and M. Schreckenberg, J. Stat. Mech. (2006) P10001.
  * (20) S. K. Baek, P. Minnhagen, S. Bernhardsson, K. Choi, and B. J. Kim, Phys. Rev. E 80, 016111 (2009).
  * (21) K. Sridhar, and K. Bicha, IJSETR 06, 0638 (2017).
  * (22) K. Schmidt-Nielsen, Scientific American. 200, 109 (1959).
  * (23) N. S. Proctor, and P J. Lynch, Manual of Ornithology, Yale University Press (1993).
  * (24) C. W. Gottschalk, and M. Mylle,
Science 128, 3324 (1958).

  * (25) H-W Lee, V. Popkov, and D. Kim,
. Phys. A: Math. Gen. 30, 8497 (1997).

  * (26) P. F. Arndt, T. Heinzel, and V. Rittenberg,
J. Phys. A: Math. Gen. 31, L45 (1998).

  * (27) M .R. Evans, D. P. Foster, C. Godreche, and D. Mukamel, Phys. Rev. Lett. 74, 208 (1995).
  * (28) N. Rajewsky, T. Sasamoto, and E. R. Speer,
Physica A 279, 123 (2000).

  * (29) T. Sasamoto, and D. Zagier,
J. Phys. A: Math. Gen. 34, 5033 (2001).

  * (30) 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.
  * (31) V. Privman, Nonequilibrium Statistical Mechanics in One Dimension (Cambridge University Press, 1997).
  * (32) J. Marro and R. Dickman, Nonequilibrium Phase Transition in Lattice Models (Cambridge University Press, Cambridge, 1999).
  * (33) G. M. Schütz, in Phase Transitions and Critical Phenomena Vol. 19, pp. 1-251, edited by C.Domb and J. L. Lebowitz (Academic Press, London, 2001) .
  * (34) M. Henkel, H. Hinrichsen, and S. Lübeck, Non-Equilibrium Phase Transitions: Volume 1: Absorbing Phase Transitions (Springer, Netherlands, 2008).
  * (35) M. Henkel, and M. Pleimling, Non-Equilibrium Phase Transitions: Volume 2: Ageing and Dynamical Scaling Far from Equilibrium (Springer, Netherlands, 2010).
  * (36) B. Derrida, S. A. Janowsky, J. L. Lebowitz and E. R. Speer, J. Stat. Phys. 73, 813 (1993).
  * (37) M. R. Evans, Y. Kafri, H. M. Koduvely, and D. Mukamel, Phys. Rev. Lett. 80, 425 (1998).
  * (38) R. Lahiri, S. Ramaswamy,
Phys. Rev. Lett. 79, 1150 (1997).

  * (39) M. R. Evans, Y. Kafri, E. Levine, and D. Mukamel,
J. Phys. A: Math. Gen. 35, L433 (2002).

  * (40) K. Mallick, J. Phys. A: Math. Gen. 29, 5375 (1996).
  * (41) B. Derrida, and M. R. Evans,
J. Phys. A: Math. Gen. 32, 4833 (1999).

  * (42) M .R. Evans, Braz. J. Phys. 30, 42 (2000).
  * (43) D. Mukamel, in Soft and Fragile Matter: Nonequilibrium Dynamics, Metastability and Flow, edited by M. E. Cates and M. R. Evans (IOP Publishing, Bristol, UK, 2000).
  * (44) M. R. Evans, Europhys. Lett. 36, 13 (1996).
  * (45) J. Krug, and P. A. Ferrari,
J. Phys. A: Math. Gen. 29, L465 (1996).

  * (46) R. A. Blythe, and M. R. Evans,
J. Phys. A: Math. Theor. 40, R333 (2007).

  * (47) F. Spitzer, Adv. in Math. 5, 246 (1970).
  * (48) M. R. Evans, and T. Hanney,
J. Phys. A: Math. Gen. 38, R195 (2005).

  * (49) O. J. O. Loan, M. R. Evans, and M. E. Cates,
Europhys. Lett., 42, 137 (1998).

  * (50) O. J. O. Loan, M. R. Evans, and M. E. Cates,
Phys. Rev. E 58, 1404 (1998).

  * (51) D. Chowdhury, and R. C. Desai,
Eur. Phys. J. B 15, 375 (2000).

  * (52) M. Bando, K. Hasebe, A. Nakayama, A. Shibata, and A. Sugiyama, Phys. Rev. E 51, 1035 (1995).
  * (53) T. S. Komatsu and S. Sasa, Phys. Rev. E 52, 5574 (1995).
  * (54) H. Hayakawa and K. Nakanishi,
Phys. Rev. E 57, 3839 (1998).

  * (55) Y. Igarashi, K. Itoh, K. Nakanishi, K. Ogura and K. Yokokawa, Phys. Rev. Lett. 83, 718 (1999).
  * (56) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * (57) A. Shapere, and F. Wilczek,
Phys. Rev. Lett. 109, 160402 (2012).

  * (58) H. Watanabe and M. Oshikawa,
Phys. Rev. Lett. 114, 251603 (2015).

  * (59) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * (60) D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao,
Annual Review of Condensed Matter Physics. 11, 467 (2020).

  * (61) J. Szavits-Nossan, and K. Uzelac,
J. Stat. Mech. (2011) P05030.

  * (62) H. C. Berg, E. coli in Motion (Springer-Verlag, LLC, New York, 2004).
  * (63) D. R. Brumley, F. Carrara, A. M. Hein, Y. Yawata, S. A. Levin, and R. Stocker, PNAS 116, 10792 (2019).
  * (64) J. Tailleur, and M. E. Cates,
Phys. Rev. Lett. 100, 218103 (2008).

  * (65) M. E. Cates, and J. Tailleur,
Annu. Rev. Condens. Matter. Phys. 6, 219 (2015).

  * (66) E. Fodor, and M. C. Marchetti,
Physica A: Statistical Mechanics and its Applications 504, 106 (2018).

  * (67) A. G. Thompson, J. Tailleur, M. E. Cates, and R. A. Blythe, J. Stat. Mech. (2011) P02029.
  * (68) R. Soto and R. Golestanian,
Phys. Rev. E 89, 012706 (2014).

  * (69) A. B. Slowman, M. R. Evans, and R. A. Blythe, Phys.Rev. Lett. 116, 218101 (2016).
  * (70) K. Malakar, V. Jemseena, A. Kundu, K. Vijay Kumar, S. Sabhapandit, S. N. Majumdar, S. Redner, and A. Dhar, J. Stat. Mech. (2018) 043215.
  * (71) E. Mallmin, R. A. Blythe, and M. R. Evans,
J. Stat.Mech. (2019) 013204.

  * (72) A. Das, A. Dhar, and A. Kundu,
J. Phys. A: Math. Theor. 53, 345003 (2020).

  * (73) I. Mukherjee, A. Raghu, and P. K. Mohanty,
arXiv:2208.05937 (2022).

  * (74) S. Jose, D. Mandal, M. Barma, and K. Ramola, Phys.Rev. E 105, 064103 (2022).
  * (75) R. Dandekar, S. Chakraborti, and R. Rajesh,
Phys. Rev.E 102, 062111 (2020).

  * (76) A. K. Chatterjee, and H. Hayakawa, SciPost Phys. 14, 016 (2023).
  * (77) L. Turner, W. S. Ryu, and H. C. Berg,
J. Bact. 182, 2793 (2000).

  * (78) M. Sidortsov, Y. Morgenstern, and A. Be’er,
Phys. Rev. E 96, 022407 (2017).

  * (79) A. Be’er, and G. Ariel, Mov. Ecol. 7, 9 (2019).
  * (80) P. F. Arndt, T. Heinzel, and V. Rittenberg,
J. Phys. A: Math. Gen. 31, 833 (1998).

  * (81) P. F. Arndt, T. Heinzel, and V. Rittenberg,
J. Stat. Phys. 97, 1 (1999).