# Steady-state selection in multi-species driven diffusive systems

Luigi Cantini  Laboratoire de Physique Théorique et Modélisation (CNRS UMR
8089), CY Cergy Paris Université, 95302 Cergy-Pontoise, France  Ali Zahra
Laboratoire de Physique Théorique et Modélisation (CNRS UMR 8089), CY Cergy
Paris Université, 95302 Cergy-Pontoise, France  Departamento de Matemática,
Instituto Superior Técnico, Av. Rovisco Pais, 1049-001 Lisboa, Portugal

###### Abstract

We introduce a general method to determine the large scale non-equilibrium
steady-state properties of one-dimensional multi-species driven diffusive
systems with open boundaries, generalizing thus the max-min current principle
known for systems with a single type of particles. This method is based on the
solution of the Riemann problem of the associated system of conservation laws.
We demonstrate that the effective density of a reservoir depends not only on
the corresponding boundary hopping rates but also on the dynamics of the
entire system, emphasizing the interplay between bulk and reservoirs. We
highlight the role of Riemann variables in establishing the phase diagram of
such systems. We apply our method to three models of multi-species interacting
particle systems and compare the theoretical predictions with numerical
simulations.

Driven diffusive systems appear in various areas across physics, chemistry,
and theoretical biology [1, 2, 3] and are widely regarded as a fundamental
playground in order to understand the behavior of complex systems away from
thermal equilibrium [4]. A classic illustration of such systems involves
particles moving within a lattice and subject to hard–core exclusion. The
introduction of a bias in their movement, simulating the influence of an
external driving force, builds up macroscopic currents in the stationary
state. A particularly relevant setting consists in putting a one-dimensional
system in contact with boundary particles reservoir, the interplay between
boundary dynamics and bulk driving leading to genuinely out of equilibrium
phenomena such as boundary induced phase transitions [5]. In this case, when
the system presents a single species of particles, a simple general principle
known as the max-min current principle [5, 6, 7, 8] allows to determine the
phase diagram for the steady state current and particle density as a function
of the boundary reservoir densities. Despite the success of this principle in
treating one-dimensional open boundary problems, its generalization to systems
containing several different species of particles has been a long-standing
challenge [9, 10, 11, 12].

The goal of the present paper is to put forward a scheme that permits to
determine the steady state average particle densities and currents of one-
dimensional multi-species driven system with open boundaries. Such a scheme is
based essentially on the sole knowledge of the bulk hydrodynamic behavior of
the model. As a starting point, similarly to the max-min principle, one
supposes the boundary densities to be known. In a systems with $n$ different
particle species, these are denoted by
$\boldsymbol{\rho}^{L}=\\{\rho_{1}^{L},\rho_{2}^{L},\dots,\rho_{n}^{L}\\}$ for
the left boundary and
$\boldsymbol{\rho}^{R}=\\{\rho_{1}^{R},\rho_{2}^{R},\dots,\rho_{n}^{R}\\}$ for
the right boundary. Then the bulk density is determined by the solution of the
associated Riemann problem at the origin (RP0)

$(\boldsymbol{\rho}^{L},\boldsymbol{\rho}^{R})\xrightarrow[]{\text{RP}_{0}}\boldsymbol{\rho}^{B}.$
(1)

As a first argument in support of this claim, we shall show that this
principle is equivalent to to Krug’s max-min current principle when applied to
the case of single-species model. We shall moreover present a further
heuristic justification of it based on a vanishing viscosity regularization of
the associated conservation laws which applies to general multi-species case.

By itself the principle (1) is not enough to determine the bulk densities
since one has at the same time to make sense of the boundary densities. If one
supposes that the boundary currents are functions of the boundary densities
alone, then current conservation through the entire systems provides the
missing conditions to completely determine both bulk and boundary densities.
We apply this scheme to three models, where we have access to the particle
currents as functions of the particle densities (which is necessary in order
to solve numerically the associated Riemann problem): 2-TASEP with arbitrary
bulk hopping rates, hierarchical 2-ASEP and a 3-TASEP. In all these three
model we find good agreement with numerical simulations.

##  I The scheme

The large scale behavior of driven diffusive system consisting of $n$ species
of particles is generally governed by a system of conservation laws

$\partial_{t}\boldsymbol{\rho}+\partial_{x}\boldsymbol{J}=0$ (2)

where the $n$ locally conserved quantities are the coarse-grained particle
densities $\boldsymbol{\rho}(x,t)=(\rho_{1}(x,t),...,\rho_{n}(x,t))$, with
associated currents
$\boldsymbol{J}(\boldsymbol{\rho})=(J_{1}(\boldsymbol{\rho}),..,J_{n}(\boldsymbol{\rho}))$.
When the system is defined on a finite interval $x\in[0,L]$ and coupled to two
reservoirs with densities $\boldsymbol{\rho}^{L}$ and $\boldsymbol{\rho}^{R}$
the system reaches in the limit $t\rightarrow\infty$ a steady state with
uniform bulk densities
$\boldsymbol{\rho}^{B}(\boldsymbol{\rho}^{L},\boldsymbol{\rho}^{R})$. We claim
that for $L\rightarrow\infty$, these bulk densities are determined by solving
a Riemann problem. Such a problem is formulated on an infinite line
$x\in\mathbb{R}$ with an initial condition consisting of two regions of
uniform densities, on the left and on the right of the origin $x=0$

$\boldsymbol{\rho}(x,0)=\boldsymbol{\rho}^{L}\mathds{1}_{x<0}(x)+\boldsymbol{\rho}^{R}\mathds{1}_{x>0}(x)\quad
x\in\mathbb{R}.$

The solution of the Riemann problem is invariant under the rescaling
$(x,t)\rightarrow(\lambda x,\lambda t)$ and therefore takes the form
$\boldsymbol{\rho}(x,t)=\boldsymbol{\rho}(\frac{x}{t})$. In particular, for
$t>0$, $\boldsymbol{\rho}(0,t)$ is independent of time, so we define:
$\boldsymbol{\rho}|_{0}(\boldsymbol{\rho}^{L},\boldsymbol{\rho}^{R}):=\boldsymbol{\rho}(0,t)$
and we call it the solution to the Riemann problem at the origin. Our claim is
that the bulk densities for the open boundary problem with given boundary
densities coincide with the solution at zero of the corresponding Riemann
problem, namely:

$\boxed{\boldsymbol{\rho}^{B}(\boldsymbol{\rho}^{L},\boldsymbol{\rho}^{R})=\boldsymbol{\rho}|_{0}(\boldsymbol{\rho}^{L},\boldsymbol{\rho}^{R})}$
(3)

The exact meaning of the boundary conditions is a mathematically subtle issue
[13, 14, 15]. We define them in an operative way as the densities of the first
and last site of the lattice, meaning that the two boundary sites can be
conceptually considered as part of their nearby reservoirs. Let us be more
specific about the boundary dynamics we shall consider. At each boundary a
particle can either enter or exit the system, or it can change its own
species. If we identify to empty sites as particles of a species $0$, the
dynamics is fully encoded in the rates
$\boldsymbol{\nu}^{L}=\\{\nu^{L}_{i,j},0\leq i\neq j\leq n\\}$ at the left and
$\boldsymbol{\nu}^{R}=\\{\nu^{R}_{i,j},0\leq i\neq j\leq n\\}$ at the right
boundary

$j\xrightarrow[]{\nu^{L}_{i,j}}i\qquad i\xrightarrow[]{\nu^{R}_{i,j}}j$

The boundary densities $\boldsymbol{\rho}^{L}$ and $\boldsymbol{\rho}^{R}$, as
well as the bulk ones are then functions of the boundary rates.

Since the boundary hopping rates are independent of the rest of the system, we
can write the current on a given boundary as a function of the density of that
boundary only

$\begin{split}J_{i}^{L}(\boldsymbol{\rho}^{L})=\sum_{j=1}^{n}\rho_{j}\nu^{L}_{ij}-\rho_{i}\sum_{j=1}^{n}\nu^{L}_{ji}\\\
J_{i}^{R}(\boldsymbol{\rho}^{R})=\rho_{i}\sum_{j=1}^{n}\nu^{R}_{ij}-\sum_{j=1}^{n}\rho_{j}\nu^{R}_{ji}\end{split}$
(4)

In the steady state, we have

$\boldsymbol{J}^{L}(\boldsymbol{\rho}^{L})=\boldsymbol{J}(\boldsymbol{\rho}^{B})=\boldsymbol{J}^{R}(\boldsymbol{\rho}^{R})$
(5)

In conclusion, eqs.(3,5) provide a system of equation enabling to determine
the bulk and boundary densities of the system.

###  I.1 Reformulation of the max-minximal Current Principle

A first argument in favor of the principle eq.(3) is the fact that in the case
of a single species of particle it coincides with Krug’s max-min current
principle. According to this principle, the steady-state current is obtained
as: [5, 7, 8, 16]

$j=\begin{cases}\max_{\rho\in[\rho^{R},\rho^{L}]}J(\rho)&\text{if
$\rho^{L}>\rho^{R}$}\\\ \min_{\rho\in[\rho^{L},\rho^{R}]}J(\rho)&\text{if
$\rho^{L}<\rho^{R}$}\\\ \end{cases}$ (6)

Let’s compare this result with what one would obtain by applying eq.(3). Let’s
start with the case where $\rho^{R}>\rho^{L}$, which corresponds to a minimum
current phase. When considering the associated Riemann problem we can assume
the current $J$ to be a convex function of the density in the interval
$[\rho^{L},\rho^{R}]$, otherwise one has to replace it with its the convex
hull in the interval $[\rho^{L},\rho^{R}]$ [17]. The solution to the Riemann
problem can be expressed as a function of $u=\frac{x}{t}$:

$\rho(u)=\rho^{L}{\bf 1}_{u<v(\rho^{L})}+\rho^{R}{\bf
1}_{u>v(\rho^{R})}+v^{-1}(u){\bf 1}_{v(\rho^{L})<u<v(\rho^{R})}$ (7)

where $v(\rho):=\frac{dJ}{d\rho}$. To compare the solution at zero with the
density predicted by the minimum current phase, we can identify three cases:
1) If $v(\rho^{L})>0$, then the solution at zero has a value of $\rho^{L}$,
and simultaneously, the minimum $\min_{\rho\in[\rho^{L},\rho^{R}]}(J(\rho))$
is reached at $\rho^{L}$. In this case, the bulk has the same density as the
left boundary, which we refer to as the left-induced phase.
2) If $v(\rho^{R})<0$, then the solution at zero has a value of $\rho^{R}$,
and simultaneously, the minimum $\min_{\rho\in[\rho^{L},\rho^{R}]}(J(\rho))$
is attained at $\rho^{R}$. This is referred to as a right-induced phase.
3) If neither of the two previous statements is true, there exists, due to the
monotonicity of the derivative, a unique value
$\rho^{B}\in[\rho^{L},\rho^{R}]$ for which $v(\rho^{B})=0$. This value
corresponds to both the Riemann solution at zero and the minimum
$\min_{\rho\in[\rho^{L},\rho^{R}]}(J(\rho))$. We refer to this situation as
the bulk-induced phase.

When $\rho^{R}<\rho^{L}$, a similar reasoning can be applied, but we replace
$J(\rho)$ with its concave hull over the interval $[\rho^{R},\rho^{L}]$. So we
conclude that the max-min current principle and the eq.(3) give the same
answer.

As an example, in the case of a single-species TASEP, we have
$v(\rho^{B})=1-2\rho^{B}$. When $v>0$, we have $\rho^{B}<\frac{1}{2}$, which
corresponds to the low-density phase, and the bulk is left-induced. The high-
density regime corresponds to a right-induced bulk density. The maximal
current phase, where $\rho=\frac{1}{2}$, is not induced from either the left
or the right.

###  I.2 Multiple Conserved Quantities

In this section we shall provide a plausibility argument for eq.(3). It will
be by no means a proof of that equation, but more support will come from the
comparison with simulations, discussed in the next section. Our argument is
based on a vanishing viscosity approach. This involves adding a diffusive
component to the current such that the total current, which remains constant
in the steady state, is given by:

$\boldsymbol{J}^{total}=\boldsymbol{J}(\boldsymbol{\rho})-\epsilon
D(\boldsymbol{\rho})\frac{\partial\boldsymbol{\rho}}{\partial x}$ (8)

Here, $\epsilon>0$ and $D(\boldsymbol{\rho})$ is a positive-definite matrix.
Since the conservation laws become locally scalar in the directions of the
eigenvectors of the Jacobian $\frac{\partial J_{i}}{\partial\rho_{j}}$ we
assume that this property extends to the viscous case, implying that
$D(\boldsymbol{\rho})$ commutes with the Jacobian. This assumption ensures a
mathematically stable regularization scheme for the boundary problem.

For the rest of the argument we shall assume that the conservation laws eq.(2)
admit $n$ independent Riemann variables
$\boldsymbol{\xi}=(\xi_{1},\dots,\xi_{n})$. These are functions of the
densities $\boldsymbol{\xi}(\boldsymbol{\rho})$, that ”diagonalize” the
conservation equations eq.(2), in the sense

$\partial_{t}\xi_{i}(x,t)+v_{i}(\boldsymbol{\xi})\partial_{x}\xi_{i}(x,t)=0,$

where it can be shown that the speeds $v_{k}$ are the eigenvalues of the
Jacobian matrix $\frac{\partial J_{i}}{\partial\rho_{j}}(\boldsymbol{\rho})$.
We remark that the existence of the Riemann variables is ensured for $n=1,2$
(for $n=1$ the Riemann variable is the density itself). Now, rewriting eq.(8)
in terms of the Riemann variables we get the ordinary differential equation:

$\frac{\partial\boldsymbol{\xi}}{\partial
x}=\epsilon^{-1}M^{-1}D^{-1}(J(\boldsymbol{\xi})-J^{total}):=F(\boldsymbol{\xi})$
(9)

where $M_{ij}=\frac{\partial\rho_{i}}{\partial\xi_{j}}$. In the limit
$\epsilon\rightarrow 0$ we have as expected $J(\boldsymbol{\xi})=J^{total}$ on
all the system, with the possible exception of microscopic regions close to
the boundaries. This means that the bulk value $\boldsymbol{\xi}^{B}$
represents a stationary point of the ODE (9), $F(\boldsymbol{\xi}^{B})=0$. In
order to determine the relation between the bulk and boundary values of each
Riemann variable, we linearize the ODE around the the stationary point. It is
not difficult to show that the Jacobian matrix $\frac{\partial
F}{\partial\boldsymbol{\xi}}$ is diagonal at the stationary point
$\frac{\partial
F_{i}}{\partial\xi_{j}}(\boldsymbol{\xi}^{B})=\epsilon^{-1}d_{i}^{-1}v_{i}\delta_{ij}$,
where $d_{i}>0$ are the eigenvalues of the diffusion matrix $D$. An
illustrative example of the field associated to the ODE for a two-component
system is in figure 2

  * •

When $v_{i}<0$, then $\xi_{i}(x)$ experiences exponential decay towards the
stationary bulk value. The decay rate is given by
$\mu_{i}=\epsilon^{-1}d_{i}^{-1}v_{i}$ . In this scenario, the bulk stationary
value is attained on the left side after a boundary layer of typical size
$1/\mu_{i}$, which is proportional to $\epsilon$. On the right boundary, the
system simply extends the bulk behavior, indicating a right-induced phase.

  * •

When $v_{i}>0$, using a similar argument we can infer that $\xi_{i}$ is
induced from the left, and the boundary layer is located on the right.

  * •

When $v_{i}=0$; the size of boundary layer diverges for finite $\epsilon$. The
flow of the ODE in the direction of the associated eigenvector indeed ceases
to be exponential and becomes rather polynomial. The bulk is therefore not
induced by any boundary, however, it belongs to the manifold
$v_{i}(\boldsymbol{\xi})=0$. We say that we are in a bulk-induced phase for
$\xi_{i}$.

This is the same result one would obtain by considering the solution of
Riemann problem at the origin. Let’s point out that the idea of looking at the
signs of eigenvalues governing the phase transition in multi-species driven
diffusive systems has already been discussed in [18], however without
reference to the Riemann variables.

##  II Application to multi-components interacting particles systems

In this section we consider three different driven diffusive systems. The
first two contain each two species of particles. More specifically the first
one is the $2$–TASEP introduced in [19, 20], while the second one is a
hierarchical $2$–species ASEP. The third model is a particular case of
$3$–species TASEP. For all this models we compare numerical simulations with
the predictions of the system of equations eq.(3) and eq.(5).

This system of equations cannot be solved analytically therefore we make use
of an iterative procedure: we begin by selecting random initial densities for
the boundaries. Then, we determine the bulk density using equation 3, which
provides information about the current. Subsequently, we calculate the
boundary densities by inverting equation 4. We continue this iteration process
between the boundaries and the bulk until convergence is achieved. However, it
is worth noticing that this algorithm may encounter cyclic trajectories. To
prevent this issue, we introduce a damping parameter $\gamma$, which should be
chosen sufficiently small. The updated equation becomes:
$\boldsymbol{x}^{n+1}=\gamma\boldsymbol{f}(\boldsymbol{x}^{n})+(1-\gamma)\boldsymbol{x}^{n}$
Here, $\boldsymbol{x}^{n}$ represents the set of variables after the $n$-th
iteration, and $\boldsymbol{f}$ represents the set of functions governing the
iterations.

###  II.1 2-TASEP with arbitrary hopping rates

This first model is a two-species generalization of TASEP, it consists of two
types of particles, denoted by $\bullet$ and $\circ$, (empty sites are denoted
by $\ast$). The hopping rates in the bulk are :

$\bullet\ast\xrightarrow[]{\beta}\ast\bullet\qquad\ast\circ\xrightarrow[]{\alpha}\circ\ast\qquad\bullet\circ\xrightarrow[]{1}\circ\bullet$

while the only non vanishing boundary rates we consider are
$\nu^{L/R}_{\bullet\ast},\nu^{L/R}_{\ast\circ},\nu^{L/R}_{\bullet\circ}$. The
currents for this model have been calculated in [21] and used in [22] in order
to study its hydrodynamic behavior and in particular to solve the
corresponding Riemann problem. Let’s recall the expression of the currents:

$\displaystyle
J_{\circ}(\rho_{\circ},\rho_{\bullet})=z_{\alpha}(z_{\beta}-1)+\rho_{\circ}(z_{\alpha}-z_{\beta})$
(10) $\displaystyle
J_{\bullet}(\rho_{\circ},\rho_{\bullet})=z_{\beta}(1-z_{\alpha})+\rho_{\bullet}(z_{\alpha}-z_{\beta})$
(11)

where $z_{\alpha}\in[0,\min(1,\alpha)]$ and $z_{\beta}\in[0,\min(1,\beta)]$
are solution of the saddle point equations

$\displaystyle\frac{\rho_{\circ}}{z_{\alpha}}+\frac{\rho_{\bullet}}{z_{\alpha}-1}+\frac{1-\rho_{\circ}-\rho_{\bullet}}{z_{\alpha}-\alpha}=0$
(12)
$\displaystyle\frac{\rho_{\bullet}}{z_{\beta}}+\frac{\rho_{\circ}}{z_{\beta}-1}+\frac{1-\rho_{\circ}-\rho_{\bullet}}{z_{\beta}-\beta}=0.$
(13)

The variables $z_{\alpha},z_{\beta}$ happen to be the Riemann variables for
this model [22]. In figure 1 (left) we reported two examples of simulations of
the $2$-TASEP on a lattice of size $L=100$ and with different values of the
model parameters. We see that the numerical result agrees very well with the
theoretical prediction obtained through the iterative solution of eqs.(3,5).
The convergence of the iterative procedure is reported on the right of the
same figure.

Figure 1:  On the left, two examples of Monte-Carlo simulation of the density
profile for 2-TASEP (continuous lines) along with the corresponding Riemann
variables (dashed lines) for a lattice of size $L=100$. The horizontal
segments represent the predicted values. On the right, the evolution of
densities for the iterative algorithm with damping $\gamma=0.01$ (up to 1000
iterations). Parameters values for top diagrams: $\alpha=0.5$, $\beta=1.5$,
$(\nu_{\bullet\circ}^{R},\nu_{\ast\circ}^{R},\nu_{\bullet\ast}^{R})=(0.29,0.08,0.07)$,
$(\nu_{\bullet\circ}^{L},\nu_{\ast\circ}^{L},\nu_{\bullet\ast}^{L})=(0.24,0.04,0.12)$.
For the bottom diagrams: $\alpha=0.4$, $\beta=0.7$,
$(\nu_{\bullet\circ}^{R},\nu_{\ast\circ}^{R},\nu_{\bullet\ast}^{R})=(0.5,0.1,0.8)$,
$(\nu_{\bullet\circ}^{L},\nu_{\ast\circ}^{L},\nu_{\bullet\ast}^{L})=(0.1,0.2,0.5)$.

####  II.1.1 Phase diagram

Following the discussion in Section I.2 we partition the phase space of the
bulk densities of this model in phases, characterized by the sign of the
functions $v_{k}(\mathbf{z}^{B})$. This a priori results in 9 phases for a
two-component system, however, hyperbolicity of the corresponding conservation
laws implies that some phases are forbidden as illustrated in the following
table

| $v_{\alpha}<0$ | $v_{\alpha}=0$ | $v_{\alpha}>0$
---|---|---|---
$v_{\beta}<0$ | $RR$ | $BR$ | $LR$
$v_{\beta}=0$ | $\times$ | $BB$ | $LB$
$v_{\beta}>0$ | $\times$ | $\times$ | $LL$

In the preceding table the first letter represents the state of $z_{\alpha}$:
L: left induced, R: right induced, B: bulk induced. The second letter is for
the state of $z_{\beta}$. The symbol $\times$ is for a forbidden phase. See
figure 2 for the result of this partitioning for the values
$\alpha=0.8,\beta=0.9$ of the bulk exchange rates.

Figure 2: Phase diagram of a 2-species TASEP ($\alpha=0.8,\beta=0.9$). The
signs on the left correspond to the velocities $v_{\alpha}$ and $v_{\beta}$ in
order. On the right, we have an example of the ODE flow exhibiting a sink
singularity in the left-induced phase and a saddle point in the mixed-induced
phase.

Numerical evidence for this diagram is reported in figure 3, where the results
of simulations are shown together with theoretical predictions with varying
parameter $\nu^{L}_{\bullet\ast}$ and all the other parameters fixed. We
notice that $z_{\beta}^{B}$ coincides with $z_{\beta}^{L}$ within the region
where $v_{\beta}<0$, and they split in the region where $v_{\beta}=0$. At the
same time $z_{\alpha}^{L}$ coincides with $z_{\alpha}^{B}$ for both regions
since $v_{\alpha}<0$.

Figure 3: Bulk and boundary densities (left) and the corresponding Riemann
variables (right) of 2-TASEP as a function of the $\nu_{\bullet\ast}^{L}$. The
crosses represent the numerical simulations, while the lines are the
theoretical predictions. For the green shaded region $v_{\beta}>0$, while for
the yellow shaded section $v_{\beta}=0$ (in both regions $v_{\alpha}<0$).
Figure 4: Bulk and boundary densities of the 2-ASEP as a function of the
$\nu_{\bullet\ast}^{L}$. The crosses represent the numerical simulations,
while the lines are the theoretical predictions. For the green shaded region
$v_{\bullet}>0$, while for the yellow shaded section $v_{\bullet}<0$ (in both
regions $v_{\circ}<0$). Parameters values:
$q=0.5,\nu^{L}_{\ast\circ}=0.9,\nu^{L}_{\bullet\circ}=\nu^{L}_{\ast\bullet}=1,\nu^{L}_{\circ\ast}=\nu^{L}_{\circ\bullet}=0.6,\nu^{R}_{\circ\ast}=\nu^{R}_{\bullet\circ}=0.1,\nu^{R}_{\circ\bullet}=\nu^{R}_{\ast\circ}=0.3,\nu^{R}_{\bullet\ast}=0.4,\nu^{R}_{\ast\bullet}=0.8$.

###  II.2 2-species ASEP and 3-species TASEP

We have considered other two models for which we have access to the exact
expressions of the hydrodynamic current as functions of the densities,.

The first model, a $2$-species ASEP, contains two species of particles and the
following bulk exchange rates:

$\nu_{ij}=\begin{cases}1&\text{if }\quad i>j\\\ q&\text{if }\quad
i<j\end{cases}$ (14)

where we have chosen the following order on the species: $\bullet>\ast>\circ$.

Although the stationary measure for a uniform state is not a product measure,
yet, it’s straightforward to write the currents-density relations since each
of the $\bullet$ and $\circ$ particles dynamics can be decoupled in the bulk:

$\begin{split}J_{\bullet}=(1-q)\rho_{\bullet}(1-\rho_{\bullet})\\\
J_{\circ}=(q-1)\rho_{\circ}(1-\rho_{\circ}).\end{split}$ (15)

From these equations it is immediate that the densities are also Riemann
variables for this model. However, the dynamics of the two species cannot in
general be decoupled on the boundaries, making the max-min principle not
applicable in this case.

The last model we have considered, a $3$-species TASEP, contains particles
with labels $(1,2,3,4)$, where the type $4$ can be seen as empty sites, and
bulk hopping rates:

$ij\xrightarrow[]{\nu_{ij}}ji\qquad\nu_{ij}=\begin{cases}0&\text{if }\quad
i>j\\\ \nu_{12}&\text{if }\quad(i,j)=(1,2)\\\ \nu_{34}&\text{if
}\quad(i,j)=(3,4)\\\ 1&\text{otherwise}\end{cases}$ (16)

The particle currents of this model can be derived from those of the
$2$-TASEP, $J_{\circ/\bullet}(\rho_{\circ},\rho_{\bullet},\alpha,\beta)$, by
making some particle identifications. Firstly, the particles $4$ and $3$ can
be seen as $\circ$, $1$ as $\bullet$ and $2$ as $\ast$, for
$\alpha=1,\beta=\nu_{12}$. Secondly, $1$ and $2$ can be seen as $\bullet$, $3$
as $\ast$ and $4$ as $\circ$ with $\alpha=\nu_{34},\beta=1$. Using densities
of particles of species $1,2$ and $4$ as independent variables one finds

$\begin{split}J_{1}&=J_{\bullet}(1-\rho_{1}-\rho_{2},\rho_{1},1,\nu_{12})\\\
J_{2}&=J_{\bullet}(\rho_{4},\rho_{1}+\rho_{2},\nu_{34},1)-J_{1}\\\
J_{4}&=J_{\circ}(\rho_{4},\rho_{1}+\rho_{2},\nu_{34},1).\end{split}$ (17)

In figure 4 and 5 we report the results for the bulk and boundary densities of
these models, obtained though simulations of a system of size $L=100$, along
with the theoretical predictions. One boundary parameter is varied
($\nu^{L}_{\bullet\ast}$ in the $2$-ASEP and $\nu^{L}_{12}$ in the $3$-TASEP)
while all the other parameters are fixed. Similarly to the case of the
$2$–TASEP seen in the previous section, we find good agreement.

Figure 5: Bulk and boundary densities of 3-TASEP as a function of the
parameter $\nu_{12}^{L}$. The crosses represent the numerical simulations,
while the continuous lines are the theoretical predictions. Parameters values:
$\nu_{12}=0.6,\nu_{34}=1.3,\nu^{L}_{13}=0.3,\nu^{L}_{14}=0.9,\nu^{L}_{23}=0.5,\nu^{L}_{24}=0.8,\nu^{L}_{34}=1,\nu^{R}_{12}=0.2,\nu^{R}_{13}=0.4,\nu^{R}_{14}=0.2,\nu^{R}_{23}=0.6,\nu^{R}_{24}=0.7,\nu^{R}_{34}=0.4$.

###  II.3 Conclusion

In conclusion, this paper introduces a method which allows to determine the
steady state average particle densities and currents of one-dimensional multi-
species driven system with open boundaries. The method, rooted in the bulk
hydrodynamic behavior of the model, extends the max-min principle applicable
to single-species models [5, 6, 7, 8]. By comparing our method’s predictions
with numerical simulations across three models, we observed good agreement.
Our analysis of bulk hydrodynamic conservation laws enables us to predict the
phase diagram, which becomes more intelligible when considering the behavior
of the Riemann variables of the model (when they exist).

The major open question pertains to the method’s domain of validity,
particularly in establishing precise definitions of boundary densities for
more general boundary conditions. The heuristic argument in favor of our
method rests on the existence of a complete set of Riemann variables in bulk
dynamics. Therefore, exploring models with more than two species, lacking this
completeness, and subjecting our method to such models, presents an intriguing
avenue for future research.

## ACKNOWLEDGMENTS

We thank Gunter Schütz for useful discussions. The work of A. Zhara has been
partially funded by the ERC Starting Grant 101042293 (HEPIQ) and completed
while he was a member of LPTM.

## References

  * Chou  _et al._ [2011] T. Chou, K. Mallick, and R. K. P. Zia, Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport, Reports on progress in physics 74, 116601 (2011).
  * Blythe and Evans [2007] R. A. Blythe and M. R. Evans, Nonequilibrium steady states of matrix-product form: a solver’s guide, Journal of Physics A: Mathematical and Theoretical 40, R333 (2007).
  * Fang  _et al._ [2019] X. Fang, K. Kruse, T. Lu, and J. Wang, Nonequilibrium physics in biology, Reviews of Modern Physics 91, 045004 (2019).
  * Schmittmann and Zia [1995] B. Schmittmann and R. K.-P. Zia, Statistical mechanics of driven diffusive systems, Phase transitions and critical phenomena 17, 3 (1995).
  * Krug [1991a] J. Krug, Boundary–induced phase transitions in driven diffusive systems, Physical review letters 67, 1882 (1991a).
  * Krug [1991b] J. Krug, Steady state selection in driven diffusive systems, in  _Spontaneous formation of space-time structures and criticality_ (Springer, 1991) pp. 37–40.
  * Popkov and Schütz [1999] V. Popkov and G. M. Schütz, Steady-state selection in driven diffusive systems with open boundaries, EPL (Europhysics Letters) 48, 257 (1999).
  * Hager  _et al._ [2001] J. Hager, J. Krug, V. Popkov, and G. Schütz, Minimal current phase and universal boundary layers in driven diffusive systems, Physical Review E 63, 056110 (2001).
  * Rákos and Schütz [2004] A. Rákos and G. Schütz, Exact shock measures and steady-state selection in a driven diffusive system with two conserved densities, Journal of statistical physics 117, 55 (2004).
  * Popkov [2004] V. Popkov, Infinite reflections of shock fronts in driven diffusive systems with two species, Journal of Physics A: Mathematical and General 37, 1545 (2004).
  * Bonnin  _et al._ [2021] P. Bonnin, I. Stansfield, M. C. Romano, and N. Kern, Two-species tasep model: from a simple description to intermittency and travelling traffic jams, arXiv preprint arXiv:2102.02486 (2021).
  * Gupta  _et al._ [2023] A. Gupta, B. Pal, and A. K. Gupta, Interplay of reservoirs in a bidirectional system, Physical Review E 107, 034103 (2023).
  * Bardos  _et al._ [1979] C. Bardos, A.-Y. LeRoux, and J.-C. Nédélec, First order quasilinear equations with boundary conditions, Communications in partial differential equations 4, 1017 (1979).
  * Dubois and Le Floch [1988] F. Dubois and P. Le Floch, Boundary conditions for nonlinear hyperbolic systems of conservation laws, Journal of Differential Equations 71, 93 (1988).
  * Mazet and Bourdel [1986] P. Mazet and F. Bourdel, Analyse numérique des équations d’euler pour l’étude des écoulements autour de corps élancés en incidence, CERT Report (1986).
  * [16] S. Katz, J. L. Lebowitz, and H. Spohn, Nonequilibrium steady states of stochastic lattice gas models of fast ionic conductors,  34, 497.
  * Osher [1983] S. Osher, The riemann problem for nonconvex scalar conservation laws and hamilton-jacobi equations, Proceedings of the American Mathematical Society 89, 641 (1983).
  * Popkov and Salerno [2011] V. Popkov and M. Salerno, Hierarchy of boundary-driven phase transitions in multispecies particle systems, Physical Review E 83, 011130 (2011).
  * Derrida [1996] B. Derrida, Statphys-19: 19th IUPAP Int,  _Conf. on Statistical Physics (Xiamen 1996) ed BL Hao (Singapore: World Scientific)_ , (1996).
  * Mallick [1996] K. Mallick, Shocks in the asymmetry exclusion model with an impurity, Journal of Physics A: Mathematical and General 29, 5375 (1996).
  * Cantini [2008] L. Cantini, Algebraic bethe ansatz for the two species asep with different hopping rates, Journal of Physics A: Mathematical and Theoretical 41, 095001 (2008).
  * Cantini and Zahra [2022] L. Cantini and A. Zahra, Hydrodynamic behavior of the two-tasep, Journal of Physics A: Mathematical and Theoretical 55, 305201 (2022).