|  |
---|---|---

École Doctorale EM2PSI (ED 405)

THÈSE DE DOCTORAT
spécialité : physique théorique

Soutenue le 29 Novembre 2022
Ali Zahra

Multi-Species Generalization of the Totally Asymmetric Simple Exclusion
Process

Integrability and Hydrodynamic Aspects

Présentée en vue de l’obtention du grade de DOCTEUR
Dirigée par : Luigi Cantini

Jury de soutenance

Kirone Mallick | Directeur de recherche | CEA (IPhT-Saclay) | Rapporteur
---|---|---|---
Gunter M. Schütz | Professeur | IST (Universidade de Lisboa) | Rapporteur
Flora Koukiou | Professeur | CNRS (CY Université) | Examinateur
Sylvain Prolhac | Maitre de conférence | IRSAMC (Université Paul Sabatier) | Examinateur
Filippo Colomo | Chargé de recherche | INFN (Sezione di Firenze) | Examinateur
Jean Avan | Directeur de recherche | CNRS (CY Université) | Examinateur
Luigi Cantini | Maitre de conférence | CNRS (CY Université) | Directeur de thèse

###### Abstract

Exclusion processes in one dimension first appeared in the 70s and have since
dragged much attention from communities in different domains: stochastic
processes, out of equilibriums statistical physics, and more recently
integrable systems. While the state of the art for a single species totally
asymmetric simple exclusion process (TASEP) can be described, from different
aspects as mature, much less is known when multiple interacting species are
present. Using tools from integrable systems and hydrodynamics in the first
place and stochastic processes in the second place, this work attempts to
study the behavior of a novel version of the model with different species of
particles having hierarchical dynamics that depend on arbitrary parameters.
While Burger’s equation famously represents the hydrodynamic limit of TASEP
with a single species, we present a counterpart coupled system of PDE
representing the hydrodynamic limit for a model with two species. The
solutions of these PDEs display a rich phenomenology of solutions best
characterized through the underlying normal modes. We discuss the associated
Riemann problem and validate our results with numerical simulations. This
system with two species can be used as a toy model for studying driven
diffusive systems with open boundaries. Using heuristics, we present results
suggesting a general principle governing the boundary induced phase diagram of
systems with multiple coupled driven conserved quantities, generalizing thus
the extremal current principle known for the case of a single driven quantity.
The integrability side of our study is mainly concerned with developing a
formalism allowing the computation of the finite-time probability distribution
of particle positions on the 1D lattice, generalizing therefore known results
for TASEP and other multi-species models. We finally study the behavior and
the impact of a single second class impurity initially located at the
interface separating two regions of different densities of first class
particles. Different limit shapes are deduced and observed. Using tools from
probability theory, we generalize the asymptotic speed properties of the
impurity for a regime of the hopping parameters.

Résumé

Les processus d’exclusion à une dimension sont apparus pour la première fois
dans les années 70 et ont depuis attiré beaucoup d’attention de la part des
communautés dans différents domaines : processus stochastiques, physique
statistique hors équilibre, et plus récemment systèmes intégrables. Alors que
l’état de l’art pour un processus d’exclusion simple totalement asymétrique
(TASEP) d’une seule espèce peut être décrit, sous différents aspects comme
mature, on en sait beaucoup moins lorsque plusieurs espèces en interaction
sont présentes. En utilisant des outils issus des systèmes intégrables et de
l’hydrodynamique en premier lieu et des processus stochastiques en second
lieu, ce travail tente d’étudier le comportement d’une nouvelle version du
modèle avec différentes espèces de particules ayant une dynamique hiérarchique
qui dépend de paramètres arbitraires. Alors que l’équation de Burger
représente la limite hydrodynamique de TASEP avec une seule espèce, nous
présentons un système couplé d’EDP représentant la limite hydrodynamique pour
un modèle avec deux espèces. Les solutions de ces EDP présentent une riche
phénoménologie de solutions mieux caractérisée par les modes normaux sous-
jacents. Nous discutons du problème de Riemann associé et validons nos
résultats par des simulations numériques. Ce système à deux espèces peut être
utilisé comme un modèle jouet pour étudier les systèmes diffusifs pilotés avec
des bords ouverts. En utilisant des heuristiques, nous présentons des
résultats suggérant un principe général régissant le diagramme de phase induit
par les frontières des systèmes avec de multiples quantités conservées
couplées, généralisant ainsi le principe du courant extrémal connu pour le cas
d’une seule quantité entraînée. L’aspect intégrabilité de notre étude concerne
principalement le développement d’un formalisme permettant le calcul de la
distribution de probabilité en temps fini des positions des particules sur le
réseau à 1D, généralisant ainsi les résultats connus pour TASEP et d’autres
modèles multi-espèces. Nous étudions enfin le comportement et l’impact d’une
seule impureté de seconde classe initialement située à l’interface séparant
deux régions de densités différentes de particules de première classe.
Différentes formes limites sont déduites et observées. En utilisant des outils
de la théorie des probabilités, nous généralisons les propriétés de vitesse
asymptotique de l’impureté pour un régime des taux.

Acknowledgment

First and foremost, I would like to thank my supervisor Luigi Cantini. Working
with him has been simply great. Our meetings have always been a source of
inspiration to me. I can’t be grateful enough to him for what I learned during
my Ph.D. Besides all of the scientific side, his support and kindness are
exceptional.

I am thankful to all the members of jury for having kindly accepted to
evaluate this work. Most of them had to make a long trip to physically attend
my defense. I would like to thank the two reporters who put a remarkable
effort into reading my manuscript and writing the reports. Their comments and
suggestions have been very useful for improving the quality of this
dissertation. I want to thank all the members of LPTM, who made this lab such
a great environment both on the professional and social levels. The lunch
breaks with Jean Avan and Genvieve Rollet are always rich in culture and
humor. I’m indebted to both of them for the generous support they offered to
me on multiple occasions. A particular thanks go to Andreas Honecker who has
been always very kind and helpful starting from my Master year and throughout
the following years. I had plenty of pleasure sharing teaching duties with Guy
Trambly de Laissardière, Jean Philippe Kownacki, Claire Pinette, Geneviève
Rollet and Andreas Honecker. They were always generous with their insightful
pedagogical hints. I would like to thank the administrator of our lab Sylvie
Villemin who is always there to help with a big smile. I want to thank my
friends and family for their continuous encouragement and support. This text
has been linguistically checked and refined thanks to the effort of my friends
Laurence Verges, Dovile Jankauskaite, Ibrahim Saideh and Marta Pedrosa García-
Moreno. Finally, I want to thank the doctoral school EM2PSI for their
financial support through the doctoral contract.

###### Contents

  1. 0 Introduction
  2. 1 Introduction to conservation laws
    1. 1.1 Scalar conservation laws
      1. 1.1.1 Introduction
      2. 1.1.2 Method of characteristics
      3. 1.1.3 Weak solution, Rankine-Hugoniot condition
      4. 1.1.4 Non-unicity of weak solutions
      5. 1.1.5 Hopf’s treatment of the Burgers equation
      6. 1.1.6 Kružkov entropy condition
      7. 1.1.7 Relation to Hamilton-Jacobi equation
      8. 1.1.8 Riemann problem
    2. 1.2 Hyperbolic Systems of Conservation Laws
      1. 1.2.1 A linear system
      2. 1.2.2 Weak solutions, the Rankine-Hugoniot condition
      3. 1.2.3 The shock curves
      4. 1.2.4 Admissibility conditions
      5. 1.2.5 Rarefaction Curves
      6. 1.2.6 T-curves
      7. 1.2.7 The Riemann Problem
      8. 1.2.8 Riemann Variables
      9. 1.2.9 Temple Class Systems
  3. 2 Hydrodynamic behavior of the two–TASEP
    1. 2.1 Introduction
    2. 2.2 Currents
      1. 2.2.1 The $z$ variables
      2. 2.2.2 Behavior at the boundary of the physical domain
    3. 2.3 Conservations laws
      1. 2.3.1 The cases $\alpha=\beta=1$ and $\alpha+\beta=1$
      2. 2.3.2 The general case: Riemann variables
      3. 2.3.3 Shocks
      4. 2.3.4 Riemann’s problem
    4. 2.4 Monte Carlo simulations
    5. 2.5 Conclusion
  4. 3 Integrable tools for the exclusion process
    1. 3.1 Exclusion process on the line
      1. 3.1.1 Exact solution for TASEP on the line
      2. 3.1.2 Adding second class particles
      3. 3.1.3 Multispecies exclusion process with arbitrary hopping rates
    2. 3.2 Exclusion process on the ring
      1. 3.2.1 Coordinate Bethe Ansatz for a defect in the ring
      2. 3.2.2 Algebraic Bethe Ansatz for The Exclusion Process on the ring
  5. 4 Boundary-induced phase transitions in multi-species driven diffusive systems
    1. 4.1 Introduction
      1. 4.1.1 Outlines and main result:
    2. 4.2 Extremal current principle revisited
      1. 4.2.1 The Riemann problem perspective
      2. 4.2.2 Vanishing viscosity approach:
    3. 4.3 The case of a multi-species driven diffusive system
      1. 4.3.1 Proof of the principle
    4. 4.4 Applications
      1. 4.4.1 Open boundaries 2-TASEP with arbitrary hopping rates
      2. 4.4.2 Limits of the method and open questions
  6. 5 Effect of a single second class particle
    1. 5.1 Introduction
      1. 5.1.1 Notations
      2. 5.1.2 Invariant measure for TASEP
      3. 5.1.3 Convergence of density field:
      4. 5.1.4 Harris graphical representation
      5. 5.1.5 Basic tool: Coupling
      6. 5.1.6 Second class particle with unity rates in a rarefaction fan
      7. 5.1.7 Matrix Product Ansatz for second class particle on the ring
    2. 5.2 A defect in a step initial profile
      1. 5.2.1 1-0 initial condition
      2. 5.2.2 Non escaping particles
      3. 5.2.3 Density field profile and the second class particle
      4. 5.2.4 $\nu-\mu$ step initial configuration
      5. 5.2.5 The case of $\nu>\mu$
      6. 5.2.6 The case of $\nu<\mu$
      7. 5.2.7 A uniform vanishing density of second class particles:
    3. 5.3 Speed process of a defect in a step initial configuration
      1. 5.3.1 Probability distribution of a second class particle of arbitrary rates in a step initial configuration
      2. 5.3.2 Appendix

##  Chapter 0 Introduction

Statistical physics at equilibrium is one of the most impressive success
stories in physics, it allows us to explain the properties of matter
surrounding us. Its mathematical foundations are well established too [1] [2].
Given a Hamiltonian system, one can find the probability distribution of
microscopic states as the one that maximizes the entropy, so for a system
coupled to a reservoir, this probability is given by the Boltzmann-Gibbs
ensemble

$P_{eq}(C)=\frac{1}{Z}e^{-H(C)/k_{B}T}$ (1)

This allows in principle to compute all physical quantities such as free
energy and correlation functions. To perform these computations, one often
needs approximation methods such as the mean-field approach, the
renormalization group, and series expansion. Analytical exact expressions are
possible only for a minor number of models giving them a major role in the
theory [3]. A prominent pioneer example is the Ising model in 2D, that was
solved by Lars Onsager in 1944 [4] and for which critical exponents were
computed exactly for the first time and were different from the mean-field
ones. This had a major impact in understanding the the critical behavior
around phase transition in equilibrium statistical physics and paved the way
for the emergence of the idea of universality where this behavior depends in
many situations only on the dimensionality and the symmetries of interactions
[3]. However, most of the collective phenomena going on in nature are out of
equilibrium. If one defines systems out of equilibrium as merely the
complementary set of systems at equilibrium, then this set is so vast that it
is not reasonable to expect it to have some common theoretical features. So we
are usually led to work within a particular setting. For instance, in closed
quantum systems, we typically consider particular schemes of a time-dependent
Hamiltonian that makes the problem tractable [5] common examples include
periodic driving, where the Hamiltonian has a time periodicity
$H(t)=H(t+\tau)$ [6]. Another one is quantum quenching, for which a system is
prepared at the ground state of a time-independent Hamiltonian, and at some
instant, we suddenly turn on a perturbing Hamiltonian and observe the dynamic
evolution of the system till thermalization [7]. These systems are often hard
to analyze analytically, and even numerically using only a classical computer,
they rather require quantum simulators such as the ones based on ultracold
atoms in optical lattices [8]. Most importantly, quantum systems in nature are
usually not closed but rather coupled to an environment, they thus exhibit
decoherence and their effective behavior collapses in many situations to non-
Hamiltonian stochastic evolution, see chapter 8 of [9] for details. This
brings us to the focus of this dissertation, which is the stochastic systems
out of equilibrium. In particular, we will be considering Markovian systems,
which are as well adapted to classical Hamiltonian systems at the mesoscopic
time scale. In such systems the stochastic evolution depends only on the
current configuration of the system, and not on its history, in other words,
these systems don’t have intrinsic memory. Mathematically, the relevant
information about the system is reduced to the set of transition rates between
the different microscopic configurations, denote $w_{C^{\prime}\rightarrow C}$
the hopping rate from the configuration $C^{\prime}$ to the configuration $C$.
Once these rates are known, one can write the evolution equation of the
probability distribution over the configuration space 111Assume this space is
countable. The more general framework is briefly mentioned in chapter 5

$\frac{dP(C)}{dt}=\sum_{C^{\prime}\neq C}\underbrace{w_{C^{\prime}\rightarrow
C}P(C^{\prime})}_{\text{gain}}-\sum_{C^{\prime}\neq
C}\underbrace{w_{C\rightarrow C^{\prime}}P(C)}_{\text{loss}}$

This is called the master equation. It can be written in a compact form:
$\frac{dP}{dt}=MP$, where $M$ is called the Markov matrix.
$M_{C,C^{{}^{\prime}}}:=w_{C^{\prime}\rightarrow C}$ for the off-diagonal
elements, and $M_{C,C}=-\sum_{C^{\prime}\neq C}w_{C\rightarrow C^{\prime}}$
for the diagonal ones. Starting form some initial probability distribution
over the configurations, and evolving in time with a Markov operator, the
system relaxes in time to a stationary state for which the probability weights
are static. 222For an infinite system, these weights might not be
normalisable, and it’s more accurate to speak about an invariant measure, as
it will be explained in chapter 3 If $\pi$ is the this stationary
distribution, it should verify $M\pi=0$. In other words, this is the
eigenvector corresponding to the zero eigenvalue. The stochastic structure of
the Markov matrix makes it so that all the other eigenvalues have a negative
real part, so they correspond to decaying modes. The eigenvalue with the
largest non zero real part provides an estimation (through its inverse) of the
typical relaxation time of the system, which is a relevant physical observable
quantity. The Markovian framework is adapted for both equilibrium and out
equilibrium systems. In the equilibrium case, although the Boltzmann-Gibbs
doesn’t tell us about the transition rates, it is always possible to assume
detailed balance, meaning that there is no net probability current between any
two configurations at equilibrium, this is expressed as:

$P_{eq}(C)w_{C\rightarrow
C^{\prime}}=P_{eq}(C^{\prime})w_{C^{\prime}\rightarrow C}$

Given a system satisfying detailed balance, if we record its time evolution,
we can’t tell in which direction the film is played. So, this is equivalent to
time reversibility. The existence of a distribution verifying the detailed
balance can be expressed as a restriction on the elements of the Markov
matrix, known as the Kolmogorov criteria 333 Although Kolmogorov implies the
existence of detailed balance, the opposite implication is valid only for
irreducible Markov chain, i.e. chains for which any state is accessible from
any other state (not necessarily directly). ,which states that around any
closed cycle of states, there is no net flow of probability, for example, for
any three configurations, we should have: $w_{1\rightarrow 2}w_{2\rightarrow
3}w_{3\rightarrow 1}=w_{1\rightarrow 3}w_{3\rightarrow 2}w_{2\rightarrow 1}$,
[10]. The simplest way to create a system out of equilibrium is to take a
system in equilibrium and to perturb it slightly so that it is driven out of
its equilibrium distribution, which breaks the detailed balance. Linear
response theory is adapted to deal with this situation [11].It applies
typically to a system with a small gradient of thermodynamic variables
inducing purely diffusive currents. Reciprocity relations over the elements of
the diffusion matrix were revealed by Onsager based on the local microscopic
time reversibility of the interactions and were the origin of his Nobel Prize
in chemistry in 1968. Another basic situation to be out of equilibrium is to
have a configuration $C$ such that $w_{C\rightarrow C^{\prime}}=0$ for all
$C{{}^{\prime}}$ and $w_{C^{\prime}\rightarrow C}\neq 0$ for some
$C^{{}^{\prime}}$, the state $C$ is called an absorbing state, once the system
reaches it, it cannot get out of it, this creates a uni-directional
probability currents towards the absorbing state and ensures being out of
equilibrium. An example for these systems is a model of the spread of an
epidemic, a recovered population would be an absorbing state. Although,
universal behaviors have been observed for absorbing state phase transitions,
most notably the direct percolation universality class where universal
critical exponents were observed [12] [13]. However, this model is not exactly
solvable, so the critical exponents are known only approximately through
numerical means. A remarkable setting where we both have exactly solvable
models with non-equilibrium steady state(NESS) [14], is the driven diffusive
systems, they can be thought of,for instance, as a gas of charged particles
with a driving electromagnetic force breaking the space isotropy and inducing
a permanent current even when the system is homogeneous [15]. They can be of
course defined on the continuous space and for an arbitrary dimension.
However, we will only consider lattice gas models defined on the lattice in
1D, this choice is justified by the availability of an exactly solvable toy
model, that is as well relevant for applications. This model is the Asymmetric
Simple Exclusion Process (ASEP), which is considered a paradigmatic model for
driven diffusive systems in general and transport models in 1D in particles.
It’s often as well described as the Ising model for out-of-equilibrium
statistical physics444The same claim is made for the directed percolation,
however, integrability makes the analogy to the Ising model more relevant for
ASEP. Let’s provide a definition and review briefly its most elementary
properties. These properties can be found in details in few classical reviews
on the topic, for instance: [16] [17].

#### The Asymmetric Simple Exclusion Process

This model is defined as a gas of identical particles on the $\mathbb{Z}$
lattice. Each particle is a random walker in continuous time, it hops forward
at a rate $p$, but only if the site in front of it is empty, and can hop
backward at a rate $q$, only if the site behind is empty. So there is a
hardcore exclusion between the particles that leads to a maximum number of one
particle per site. The initial motivation and context of the appearance of the
model will be mentioned latter in this introduction. A simple particular case
is when $p=q$ then the model is called SSEP (Symmetric Simple Exclusion
Process), this model is not out of equilibrium, it’s easy to understand that
there is no average current and that detailed balance is conserved, however,
this particular case is still relevant either from a mathematical point of
view where it has been historically the first case to be solved exactly by
mapping it to spin chains, or sometimes it can be seen as a critical system
where we have a transition between non-equilibrium and equilibrium. Another
particular case that is interesting for many reasons is when the particles
move only in one direction, take for instance $q=0$, we speak here about the
Totally Asymmetric Simple Exclusion Process (TASEP). We can set $p=1$ by a
change of the scale of time. This particular case allows, in many cases, for
exact computations that are much harder for the general ASEP, so it provided
the simplest out-of-equilibrium exactly solvable model.

Figure 1: The corner growth process as an interpretation of TASEP

What adds to the interest of TASEP is that it has different interpretations.
For instance, it can be mapped to a surface growth model. This was first
pointed out by Rost [18]. Let a 2D surface with an upper boundary represented
by an affine continuous function $h(x,t)$ defined up to an additive constant
and verifying:

$h(j+\frac{1}{2},t)-h(j-\frac{1}{2},t)=\begin{cases}-1\quad&\text{if the site
at $j$ is occupied}\\\ 1\quad&\text{if the site at $j$ is empty}\\\
\end{cases}$ (2)

One can understand that the time evolution of the TASEP corresponds to a
random growth process of the surface, figure 1.

$(a)$$q$$p$$\times$$\times$$(b)$$\alpha$$\delta$$Left$$Res.$$\beta$$\gamma$$Right$$Res.$$(c)$
Figure 2: Different boundary conditions for the exclusion process. (a)
Periodic boundary condition where the sites are identifiable to
$\mathbb{Z}/L\mathbb{Z}$, $L$ being the number of sites. (b) ASEP on the line;
the lattice is $\mathbb{Z}$. (c) Open boundary conditions, where we have a
finite number of site. The system is coupled to a reservoir on the left where
particles can hop inside at a rate $\alpha$ if the first site is empty. If it
is not, then the occupying particle can escape the system at a rate $\delta$.
Similar mechanism occurs on the right but with different rates.

Another interesting interpretation of TASEP is in terms of queuing theory: The
particles can be though of as servers and the voids as clients, so each server
has a queue of clients in front of it. When a particle jumps, a client is
served, this client will queue up in the queue belonging to the following
server, and waits again for its turn. This image was exploited by [19] who
found the invariant measure for TASEP. More details are in the introduction of
chapter 5.

Different types of boundary conditions are possible for ASEP, each has its own
interest. They are illustrated in figure 2. Let’s look at the most classical
properties ASEP on the ring, and TASEP with open boundary conditions.

##### ASEP on the ring

The periodic boundary condition is the simplest. it’s almost a trivial, yet
pedagogical exercise to find the stationary state for ASEP with periodic
boundary conditions. Consider a configuration composed of $l$ blocks of
particles, where a block is a set of adjacent particles surrounded by voids.
The system can leave or join the configuration only by the front or the backs
of a block:

$\circ\bullet\bullet\bullet\circ$$q$$p$$\circ\bullet\bullet\circ\bullet$$\bullet\circ\bullet\bullet\circ$$p$$q$

Now it’s easy to understand that if we chose a uniform probability
distribution for the configurations (each configuration has a probability
$P_{eq}$), then each of the escaping rate and the entering rate will be equal
to $l(p+q)P_{eq}$ which leads to a stationary system. If there are $M$
particles and $N$ sites, each configuration will have the probability:
$P_{eq}=1/\binom{N}{M}$. Now the current can be found exactly by choosing a
site and counting the number of configurations such that this site is occupied
and followed by a void or the other way around:

$J=pE(\bullet\circ)-qE(\circ\bullet)=(p-q)\dfrac{\binom{N-2}{M-1}}{\binom{N}{M}}=(p-q)\frac{M}{N}(\frac{M-N}{M-1})\rightarrow(p-q)\rho(1-\rho)$

where the limit is taken for infinite $N$ and $M$ and keeping a fixed ratio
$\rho:=\frac{M}{N}$ which is the average density. We notice that this
expression in the limit of a large system is the same as one obtained by a
mean field. We will briefly see in chapter 5 that this is due to the fact that
the product measure is invariant for ASEP in an infinite system.

##### Hydrodynamic behaviour of TASEP

Now imagine a system with an average local coarse-grained density changing
over space and time $\rho(x,t)$, regardless of the boundaries,consider the
TASEP case, we can write a conservation equation associated with the previous
expression of the current,

$\partial_{t}\rho+(1-2\rho)\partial_{x}\rho=0$

This equation is called the non-viscous Burgers equation. Its more precise
meaning will be given later in this introduction. However, note that $1-2\rho$
expresses the speed of the front wave around the density $\rho$. Note that it
is a decreasing function of the density. If we have an increasing initial
profile of density over space, then the upper parts will move faster than the
lower parts, creating an even steeper profile, till we finally reach a
discontinuous profile that is called a shock, cite 3. This shock is not
static, if the density on its left is $\rho^{L}$ and on its right is
$\rho^{R}$ then the speed of the shock is given $1-\rho^{L}-\rho^{R}$, as we
will see in chapter 2. On the other hand, if the initial profile is decreasing
as a function of the space, then its slope will get even smaller, and the
solution will stay regular, more details will be provided in chapter 2.

$x$$\rho$0.5$t=2$$t=1$$t=0$ Figure 3: Illustation of the formation of the
shocks in Burgers equation as a result of the group velocity $v=1-2\rho$

##### TASEP with open boundaries

Consider TASEP with open boundaries where particles can hop inside the system
from the left at rate $\alpha$ if the first site is empty, and can leave the
system from the right at rate $\beta$. The left boundary behaves as if has a
density $\rho^{R}=\alpha$, and the right boundary behaves as if it has a
density $\rho^{L}=1-\beta$. Now, it is possible to sketch most of the behavior
of the system using a heuristic hydrodynamic approach based on the front wave
speed $v=1-2\rho$. Note first that if $\alpha<\frac{1}{2}$ then
$v^{L}=1-2\alpha>0$, so there is a kinetic wave at density $\alpha$ trying to
penetrate the system from the left. If $\beta<\frac{1}{2}$ then
$v^{R}=1-2(1-\beta)=2\beta-1<0$, so now the kinetic wave of density $1-\beta$
is trying to enter from the right. Now, we can distinguish four cases, figure
4

$\alpha$$\beta$$\frac{1}{2}$$1$$1$MCLDHD Figure 4: Boundary induced phase
diagram of TASEP. LD: Low-density phase. HD: high-density phase. MC: Maximal
current phase. The blue line represents boundaries where a phase order phase
transition occurs. The orange line represents a second-order phase transition.
The bulk density can be regarded as the order parameter

  * •

$\alpha<\frac{1}{2}$ and $\beta>\frac{1}{2}$, only the wave from the left is
entering the system, and it will reach the bulk, so we have a system dominated
by a density $\alpha<\frac{1}{2}$, this phase is called a low-density phase
(LD)

  * •

$\alpha>\frac{1}{2}$ and $\beta<\frac{1}{2}$, the opposite of the previous
situation, the bulk density will be $1-\beta>\frac{1}{2}$. This is called a
high-density phase (HD)

  * •

$\alpha<\frac{1}{2}$ and $\beta<\frac{1}{2}$, both of kinetic waves are
entering the system, so will create a shock that moves at a speed
$1-\rho^{L}-\rho^{R}=\beta-\alpha$. If this speed is positive then the left
boundary dominates the bulk, extending the low-density phase. Otherwise, the
right boundary wins and the system is in the high-density phase.

  * •

$\alpha>\frac{1}{2}$ and $\beta>\frac{1}{2}$, then both of the waves are
leaving the system, creating a phase where the current is maximal (MC) and the
bulk density is $\frac{1}{2}$

This qualitative hydrodynamic approach has been confirmed by an exact solution
[20] by solving recursion relations of the probability profile on the size of
the system. The boundary induced phase transition for TASEP is one of the
simplest for driven diffusive systems, and it paved the way for developing a
more general principle describing the boundary-induced phase transitions of
any system with a single driven quantity [21] [22], [23], [24]. We will see
how it will be generalized in chapter 4 for systems with multiple coupled
driven quantities.

##### A brief historical perspective:

Let’s now stop at the most prominent stations during the lifetime of the model
that was an inspiration to our work:

  1. 1.

In 1968, MacDonald, Gibbs, and Pipkin first proposed ASEP [25]in the context
of transport in biology modeling the situation of multiple enzymes copying
sequentially from the same DNA template. They actually introduced a more
general version of ASEP where the exclusion rule extends over $L$ neighboring
sites rather than just one, so particles can’t have a distance less than $L$
555Or equivalently, as the original formulation, they are not particles, but
segments of length $L$. They found using a mean-field analysis, the expression
of the current for a uniform $\rho$ density system. In addition, they sketched
the main features of the behavior of ASEP with open boundaries using the
hydrodynamic approach.

  2. 2.

In 1978, Alexander and Holstein [26] mapped the master equation for SSEP to
the Heisenberg spin chain. Since this spin chain was diagonalized exactly by
Bethe in 1931, this sparked the interest of the integrability community in
particle systems. The mapping was later extended to various other one-
dimensional reaction-diffusion processes. Check [27] for an early review which
highlighted the underlying Heck algebra that is common among the evolution
operators of that family of models. In 1992 Gwa and Spohn [28] mapped ASEP to
XXZ spin chain, which allowed them to diagonalize the Markov matrix using
Bethe Ansatz and to estimate the relaxation time. More details will be
provided in chapter 3.

  3. 3.

In 1981, Rost [18] noticed that if time and space are scaled in the same way
(in other words, you compress the space and accelerate the time with the same
large factor), the corresponding density profile converges to a deterministic
limit shape given by Burgers equation. The limit shape is properly defined for
the height function through a hydrodynamic scaling:

$\rho(x,t)=\partial_{x}\lim_{\epsilon\rightarrow 0}\epsilon
h(x\epsilon^{-1},t\epsilon^{-1}),$ (3)

where the macroscopic density is the physical solution of Burgers equation.
666 As we will see in the next chapter, a weak form of Burgers equation admits
unstable solutions that we refer to as non-physical Although this result is
true for any initial condition, Rost proved it with a step initial profile,
all negative sites are occupied, and all positive sites are empty. This
initial profile plays a role comparable to that of quench in quantum out-of-
equilibrium systems. The limit shape with this initial condition is:

$\rho(x,t)=\begin{cases}1\quad\text{if}\quad x<-t\\\
\frac{1}{2}(1-\frac{x}{t})\quad\text{if}\quad-t<x<t\\\ 0\quad\text{if}\quad
x>t\end{cases}$ (4)

And the height will have a limit shape:

$\lim_{t\rightarrow\infty}\frac{h(vt,t)}{t}=\begin{cases}|v|\quad\quad\quad\quad\text{if}\quad\quad|v|<1\\\
\frac{1}{2}(v^{2}+1)\quad\text{if}\quad-1<v<1\\\ \end{cases}$ (5)

Although the Burgers equation existed much earlier, it was the first time the
exclusion process was proposed as a microscopic description of the Burgers
equation.

  4. 4.

In 1991, P.A. Ferrari, while studying the shocks fluctuation of ASEP [29],
introduces a second class particle, this particle jumps as a normal particle
when the following site is empty: $20\rightarrow 02$, but the normal particle
(named first class) see it as void, and can thus swap with it: $12\rightarrow
21$. The second-class particle can’t overtake the first-class particle, hence
the terminology. This particle was introduced as a means to identify
microscopically the shock. It was inspired by the basic coupling technique
introduced by Ligget [30] for TASEP. In the same year Ferrari, Kipnis and
Saada proved that if a second-class particle is added to the origin of a
rarefaction fan, it will choose a random asymptotic speed within the available
ones with a uniform measure. Soon the second-class particle attracts the
attention of a wider audience. In 1993, Derrida, Janowsky, Lebowitz, and Speer
[31] determine the shock profile as seen from the perspective of a second-
class particle. Not much later, the second-class particle acquires an interest
on its own besides its role as a theoretical mean. In 1996, Derrida [32] and
Mallick [33] generalized this concept into a defect or impurity, which is a
second-class particle that can jump with an arbitrary hopping rate and can be
taken over with another arbitrary rate. This was the birth of a new model: the
two species TASEP.

  5. 5.

In 1997 Schütz [34] obtained for TASEP on the infinite line, an exact
expression for the conditional probability
$P(x_{1},...,x_{N};t|y_{1},...,y_{N};0)$ of $N$ particles being at positions
$\\{x_{1},...,x_{N}\\}$ at time $t$ given their initial positions
$\\{y_{1},...,y_{N}\\}$ at time zero. The result is obtained using Bethe
Ansatz and was expressed as an $N\times N$ determinant. This was an early
result connecting TASEP to the domain of integrable probability. Tracy and
Widom generalized it to ASEP [35] [36]. Latter Tasep with second class
particles was treated by Chatterjee and Schütz [37]. This will be reviewed and
extended in chapter 3.

  6. 6.

In 1993, the exact phase diagram for of TASEP with open boundaries was derived
in two independant papers: Schütz and Domany [20] solved recursion relations
on the size of the system for the stationary state, allowing its explicit
expression, and discussed the phase diagram in terms of the dynamics a domain
wall. The second paper is [38] where Derrida, Evans, Hakim and Pasquier
determined this stationary state using a Matrix Product Ansatz (MPA)
formulation. This opened the door for a series of cases where a non-
equilibrium steady state is expressed in a matrix product form. For a
pedagogical review check [16]. A brief explanation of MPA will be provided in
chapter 5.

  7. 7.

In 1999, Johansson [39] revealed a connection to Random Matrix Theory (RMT)
that triggered an impressive quantity of subsequent investigations. This is
important to out-of-equilibrium statistical physics in particular because RMT
has been a gold mine for universal behaviors. Let’s state the main result:

$\lim_{t\rightarrow\infty}Prob\big{(}h(vt,t)>\underbrace{\frac{1+v^{2}}{2}t}_{\text{Limit
shape}}-s\underbrace{\frac{(1-v^{2})^{2/3}}{2^{1/3}}t^{1/3}}_{\text{Fluctuations}}\big{)}=F_{GUE}(s)$
(6)

where $F_{GUE}$ is the cumulative Tracy–Widom distribution, precisely, it is
the distribution of the rescaled largest eigenvalue $\lambda_{max}$ of random
matrix sampled from the Unitary Gaussian Ensemble. If $n\times n$ is the size
of the matrix, this eigenvalue grows as $\sqrt{2n}$ and fluctuates with a
standard deviation of $n^{-\frac{1}{6}}$. Then we have:
$F_{GUE}(s):=\lim_{n\rightarrow\infty}Prob((\lambda_{max}-\sqrt{2n})\sqrt{2}n^{\frac{1}{6}}<s)$.
The term $\frac{1+v^{2}}{2}$ represents the speed of the growth. The most
important information in that equation is the exponent in $t^{1/3}$. It has
been previously conjectured that this growth model belongs to the KPZ
universality class, and thus fluctuates as $t^{1/3}$ however, it was the first
time this was proven and the only model for which it was proven rigorously. To
make sure that 6 is appreciated correctly, one can compare to the central
limit theorem, where $t^{1/3}$ plays the role of $t^{1/2}$ and $F_{GUE}$ plays
the role of the integral of a Gaussian. The poof of the previous result is
based on combinatorics, there is a correspondence with the problem of the
distribution of the length of the longest increasing subsequence in a random
permutation that has the same limiting distribution [40].

##### From single-species to multi-species TASEP

The exclusion process, and TASEP in particular, is far from being only a
mathematical model in 1D that theoretical physicists get excited about. It is
used as a mesoscopic vehicle model in the field of traffic flow, for instance,
the phase diagram for TASEP with open boundaries is celebrated in the traffic
literature [41] [42]. However, this model is too idealized to be applied to
real systems, It’s rather only suited for one-lane, one-direction identical
vehicles on a homogeneous freeway with no accidents. We all know how is it in
daily situations. A similar narrative can be made regarding intracellular
transport in biology where the exclusion process is still relevant, for
instance in molecular motor proteins moving along micro-tubules filaments [43]
[44], but again, real transport in cellular biology involves complex phenomena
not counted by TASEP, such as the existence of multiple types of molecules
transported on the same filament. See figure 5 for an example of transport
molecular motors in neurons, studying this transport phenomenon has been
relevant for the understanding of brain function, development, and disease
[45].

Figure 5: Molecular transport on an axon, the main nerve fiber of a neuron,
featuring different types of molecular motors, adapted from [45]

Hence the need for a model taking into account the presence of different types
of particles that have different rates and that can swap between each other.
In this model the exclusion rule is still valid as well as the local update.
Different species swap with arbitrary entra-species rates:

$(\bullet{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}\rightarrow{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}\bullet)\;\text{with
rate
}\>\tau_{\bullet{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}}$

For this model to be exactly solvable, some restrictions have to be obeyed by
the rates as we will see in chapter 3. In the particular case of two species +
void, it’s enough to assume hierarchy for the model to be exactly solvable,
meaning that if we denote the particles as $0,1,2$, the only possible swaps
are: $10\rightarrow 01$, $20\rightarrow 02$, $12\rightarrow 21$ with arbitrary
rates. This model was first introduced by Derrida [32] and Mallick [33] for a
single second-class particle. Cantini later found the currents for an
arbitrary number of defects [46]. This model, besides its applications,
represents a much richer spectrum of phenomenology compared to TASEP, even
when the simplest questions are asked. The objective of this dissertation is
to be a building block for the bulk of knowledge for the 2-species exclusion
process.

##### Novelties of this work:

  * •

Addressing the hydrodynamic behavior of two species TASEP. Although the
rigorous convergence to the limit shape is a mathematically subtle question,
we will rather make use of the integrability of the model that provides the
currents and solves coupled conservation partial differential equations. The
solutions are substantially more complex and rich than the TASEP one. This is
the subject of our publication [47] which is included as in chapter 2.

  * •

We provide in chapter 3 a framework allowing the calculations of finite time
conditional probability for the position of a finite number of particles of
multiple species. The formalism can be though of as a stochastic vertex model
and leads explicit formulas in particular situations, generalizing the work of
Schütz [34] and Schütz et al. [37]

  * •

We investigate in chapter 4 a method that allows to determine the steady state
of a driven diffusive system with multiple driven coupled quantities,
generalizing thus the extremum current principle proposed by Krug [21], Schütz
and others [22], [23], [24]. This method is operational even for models where
the stationary measure is not a product measure, completing thus other method
proposed in [48] [49] [50]. We apply this formalism to multiple particle
models, 2-TASEP being one of them.

  * •

In Chapter 5, we treat the question of the interaction between a defect
particle and a density field for TASEP on the line with a Riemann initial
condition. Besides the different phenomenology encountered, we expand the
proof of the uniform density for the asymptotic speed for the case of a step
initial profile.

For the unfamiliar reader, the first chapter is dedicated to providing all the
necessary tools from the domain of conservation laws. This is required for the
second chapter as well as the fourth one.

Each of chapters 3,4 and 5 will be the core of a future separate publication.

##  Chapter 1 Introduction to conservation laws

In 1757, Leonhard Euler wrote in his memoir ”Principes généraux du mouvement
des fluides” an equation for the conservation of momentum and another for the
mass. These equations were among the first partial differential equations ever
written [51, 52] and raised the initial problems that led later to the
development of the domain of conservation laws with widespread applications in
physics and chemistry. From a mathematical point of view, they are often
qualified as hyperbolic due to their wavelike solutions. Yet, they are famous
for having shocks singular solutions, requiring mostly an ad-hoc mathematical
framework and placing them often in the last chapter of PDE textbooks. Despite
being an old subject, research is still active in the domain [53]. Although
the space multi-dimensional conservation laws are nowadays an exciting
frontier of research, we restrict our presentation to 1D space, focusing
mainly on the aspects related to the needs of the other chapters.

This essay starts with a discussion of scalar conservation laws in section
1.1, with an emphasis on the techniques that are generalizable to non-scalar
systems with multiple coupled conserved quantities. In particular, the
stability conditions for the weak solution are discussed in details. Burgers
equation is used as a toy example for the scalar laws, the version used here
is $u_{t}+uu_{x}=0$ which is slightly simpler than the TASEP one but
completely equivalent. A flavor of the vanishing viscosity method is given in
section 1.1.5 during Hopf’s treatment of the Burgers equation. This method
will be relevant to chapter 4 when dealing with multiple conservation laws in
a system with open boundaries. In section 1.2, we review the most classical
features of systems of conservation laws, this provides the background
necessary for chapter 2 where we solved a system with two conserved quantities
resulting from a scaled two species TASEP model. We finish this section with a
brief discussion of a particular family of conservation laws known as the
Temple class, which has a curious connection with integrable models that we
briefly investigate on the hydrodynamic level.

Despite that in section 1.1.7, we present the Hopf-Lax formula that allows
formally to treat a wide class of initial conditions, the focus is later given
only to the Riemann initial condition. The relevance of the Riemann problem
can be compared to quenching in a quantum system; it’s a popular procedure
that provides insights into the dynamical behavior of the system and has been
used in chapter 2 for the 2-species TASEP.

This chapter is largely mathematical and mostly based on classical texts [54,
55, 56], [55], [56], [57], [58], [59], [60], [61].

###  1.1 Scalar conservation laws

####  1.1.1 Introduction

In this part we consider a single unknown function: $u(x,t):\mathbb{R\times
R^{+}}\rightarrow\mathbb{R}$ that represents a density satisfying a
conservation laws with initial data at $t=0$:

$\begin{split}&u_{t}+(f(u))_{x}=0\\\ &u(x,0)=u_{0}(x)\end{split}$ (1.1)

with $f$ in the class $\mathcal{C}^{1}$, representing the flux of $u$. We are
concerned here in investigating the solutions of this problem in the most
general setting. To give an initial flavor, although not representative of
general solutions, let’s start with the trivial case of a linear flux.

##### A linear flux:

The most simple case is when $f(u)=cu$. The equation becomes:

$u_{t}+cu_{x}=0$ (1.2)

This means that the directional derivative of $u$ in the direction $(1,c)$ is
zero, so $u$ is constant over the lines $x=ct+x_{0}$:

$u(ct+x_{0},t)=u(x_{0},0)=u_{0}(x_{0})$

With a change of variable we have:

$u(x,t)=u_{0}(x-ct)$ (1.3)

So the initial profile will be just moving at a constant speed $c$. Note that
if we write the general equation in the form: $u_{t}+f^{\prime}(u)u_{x}=0$ and
consider an initial profile that is uniform with an infinitesimal perturbation
around $\tilde{u}$ $u_{0}(x)\approx\tilde{u}$ then it will evolve translating
with the speed $f^{\prime}(\tilde{u})$, we call this speed, the speed of
perturbations.

If the flux is non-linear, the differential equation is said to be quasi-
linear (A fully non-linear equation requires a non-linearity of the highest
derivative: i.e. $u_{x}$ or $u_{t}$). In the next paragraph, we remind a
general method that is used not only for conservation laws but for a wider
class of non-linear first order PDE.

####  1.1.2 Method of characteristics

The method of characteristics consists of partitioning the variables’ space
into a family of curves where the PDE transforms into a system of ODEs
(Ordinary differential equation) on the curves. It is adapted for the general
class of non-linear first order equations, i.e. equations of the form:
$H(Du,u,\textbf{x})=0$ defined on an open domain $\textbf{x}\in
U\subset\mathbb{R}^{n}$ and subject to a boundary condition $u=g$ on a curve
$\Gamma\subset\partial U$. The characteristics are the three functions of a
real parameter:

$\begin{split}&\textbf{x}(s)\\\ &z(s):=u(\textbf{x}(s))\\\
&\textbf{p}(s):=Du(\textbf{x}(s))\end{split}$ (1.4)

Deriving $F$ with respect to $x_{i}$ gives:

$\sum_{j}\frac{\partial H}{\partial p_{j}}u_{x_{j}x_{i}}+\frac{\partial
H}{\partial z}u_{x_{i}}+\frac{\partial H}{\partial x_{i}}=0$ (1.5)

We can identify the first term with
$\dot{p}_{i}(s)=\sum_{j}u_{x_{i}x_{j}}\dot{x}_{j}(s)$ providing that we
identify $\frac{\partial H}{\partial p_{j}}$ with $\dot{x}_{j}$. So finally,
we have this system of ODE:

$\begin{split}&\dot{\textbf{x}}=D_{p}H\\\
&\dot{\textbf{p}}=-D_{x}H-\textbf{p}D_{z}H\\\
&\dot{z}=\textbf{p}D_{p}H\end{split}$ (1.6)

(Note that if we forget about the third equation and the second term of the
second equation, we get Hamilton-Jacobi equations. We will come back to this
later). This equivalence between the PDE and the set of ODEs is formally valid
for regular solutions $u\in\mathcal{C}^{2}$. Under this condition, the Cauchy
problem of the ODE has a unique solution for sufficiently regular $H$
(Lipschitz) providing that the boundary condition is compatible with the
characteristics.

##### Application to the scalar conservation law

For our purpose, it’s quite simple: $\textbf{x}=(x,t)$,
$\textbf{p}=(u_{x},u_{t})$, $H(u_{t},u_{x},u,x,t)=u_{t}+f^{\prime}(u)u_{x}$
The characteristics are:

$\begin{split}&\dot{t}=0\\\ &\dot{x}=f^{\prime}(z)\\\ &\dot{z}=0\end{split}$
(1.7)

They form a closed system. We can obviously use the time as a parameter:

$\begin{split}&\frac{dx}{dt}(t)=f^{\prime}(u(x,t))\\\
&\frac{d}{dt}u(x(t),t)=0\end{split}$ (1.8)

So $u$ is constant all over the characteristics, and they are simply straight
lines:

$\begin{split}&x(t)=f^{\prime}(u_{0}(x_{0}))t+x_{0}\\\
&u(x(t),t)=u_{0}(x_{0})\end{split}$ (1.9)

If the initial condition is smooth then the solution at time $t>0$ is still
smooth as far as the characteristics don’t intersect, so we have a classical
$\mathcal{C}^{1}$ solution for $0\leq t<T$. If $f^{\prime}$ is
$\kappa_{1}$-Lipschitz and $u_{0}$ is $\kappa_{2}$-Lipschitz then the first
intersection of characteristics will appear at $T=\kappa_{1}\kappa_{2}$. At
the points $(x,t)$ of the intersection of characteristics, the value of the
solution is not defined since different values carried from different
characteristics contradict. The limit of $u(x,t)$ at the point of intersection
of characteristics depends on the path, so the solution forms a finite
discontinuity, figure 1.1.

$x$$t$$x_{0}$A characteristic transporting the density $u(x_{0})$ Figure 1.1:
Illustration of singularity formation through characteristic intersections

####  1.1.3 Weak solution, Rankine-Hugoniot condition

Clearly, we need a weaker interpretation of the equation that takes into
account discontinuous solutions. A possible way is to consider the equation in
the distribution sense, so $u$ would be a distribution acting on a Schwartz
space (a space of test functions $\phi(x,t)$ of the class
$\mathcal{C}^{\infty}$ with compact support)

$\int_{0}^{\infty}\int_{-\infty}^{\infty}(u_{t}(x,t)+f(u(x,t))_{x})\phi(x,t)dxdt=0$
(1.10)

The relevance of this writing is that it allows performing the integration by
part:

$\int_{0}^{\infty}\int_{-\infty}^{\infty}u(x,t)\phi_{t}(x,t)+f(u(x,t))\phi_{x}(x,t)dxdt=0$
(1.11)

This equation may be called the integral form of the conservation law. It
doesn’t impose a regularity restriction on the solutions. Functions verifying
this equation are called weak solutions.

Let’s consider now a situation where we have a solution $u(x,t)$ that is
regular all over $\mathbb{R}\times\mathbb{R}^{+}$ except on some continuous
path $\Gamma$ parameterized by $x=s(t)$ (so we assume that there is a single
finite discontinuity at each instant). We are interested in describing the
behavior of this path. Let’s assume as well the following limits exist:

$\begin{split}&\lim_{x\xrightarrow{<}s(t)}u(x,t):=u^{L}(t)\\\
&\lim_{x\xrightarrow{>}s(t)}u(x,t):=u^{R}(t)\end{split}$ (1.12)

The path divides the domain $\mathbb{R}\times\mathbb{R}^{+}$ into two
subdomains, one located on its left $\Omega^{L}$ and another on its right
$\Omega^{R}$. figure 1.2

$x$$t$$\vec{n}$$\Omega^{L}$$\Omega^{R}$$\Gamma$ Figure 1.2: A singular path
propagating in the $(x,t)$ domain

We can decompose eq. (1.11) into :

$\int_{\mathbb{R}\times\mathbb{R}^{+}}u\phi_{t}+f(u)\phi_{x}dxdt=\int_{\Omega^{R}}u\phi_{t}+f(u)\phi_{x}dxdt+\int_{\Omega^{L}}u\phi_{t}+f(u)\phi_{x}dxdt=0$
(1.13)

Let’s integrate by part the first term:

$\int_{\Omega^{L}}u\phi_{t}+f(u)\phi_{x}dxdt=\int_{\Omega^{L}}u_{t}\phi+f(u)_{x}\phi
dxdt-\int_{\partial\Omega^{L}}u\phi n^{t}+f(u)\phi n^{x}dxdt\\\ $ (1.14)

$=-\int_{\Gamma}(u^{L}n^{t}+f(u^{L})n^{x})\phi dxdt$ (1.15)

Where $(n^{x},n^{t})$ is a unit vector normal to the boundaries

We can treat the term on the right in a similar fashion except that we will
have a minus sign from $(n^{x},n^{t})$ since we will use the same normal
vector as previously. so finally, we get:

$\int_{\Gamma}((u^{L}-u^{R})n^{t}+(f(u^{L})-f(u^{R}))n^{x})\phi dxdt=0$ (1.16)

Which means:

$(u^{L}-u^{R})n^{t}+(f(u^{L})-f(u^{R}))n^{x}=0$ (1.17)

knowing $(n^{x},n^{t})$ allows us to have the tangent to $\Gamma$, which gives
the derivative:

$\frac{ds}{dt}=\frac{-n^{t}}{n^{x}}$ (1.18)

This is nothing but the speed of the shock, let’s note it $\sigma$. So we
reach the famous Rankine-Hugoniot formula:

$(u^{L}(t)-u^{R}(t))\sigma(t)=f(u^{L}(t))-f(u^{R}(t))$ (1.19)

There is a simple way to grasp this identity, simply by imagining the shock as
a level of water in a 2D tank that has a source and a sink. Each side
represents the rate of filling of the tank expressed in two different ways.
The problem with weak solutions is that they are not always unique as we will
see in what follows.

####  1.1.4 Non-unicity of weak solutions

While the strong form of the conservation equation doesn’t always have a
solution, the weak form might have more than one solution with the same
initial data. Let’s give the Burgers equation as an example:

$u_{t}+uu_{x}=0$ (1.20)

With the initial condition:

$u_{0}(x)=\mathds{1}_{x>0}(x)$ (1.21)

The flux for this equation is $f(u)=u^{2}/2$. It admits two weak solutions: a
regular one, called a rarefaction fan:

$u(x,t)=\frac{x}{t}\mathds{1}_{0<x<t}(x)+\mathds{1}_{x>t}(x)$ (1.22)

And a shock with a speed $\frac{1}{2}$

$u(x,t)=\mathds{1}_{x>\frac{t}{2}}(x)$ (1.23)

One can argue that the second solution is not stable. Consider for instance
this small (in the sens of $\mathcal{L}^{1}$) perturbation of the initial
profile:

$u_{\epsilon}(x)=\frac{x}{\epsilon}\mathds{1}_{0<x<\epsilon}(x)+\mathds{1}_{x>\epsilon}(x)$
(1.24)

It’s clear by looking at the characteristics that this profile will evolve in
time like the first solution (the fan) and will thus divert from the shock
solution corresponding to $\epsilon=0$. We will see later more formal notions
of stability of solutions.

If we now change the initial condition to this one:

$u_{0}(x)=\mathds{1}_{x<0}(x)$ (1.25)

Then we have this weak solution:

$u(x,t)=\mathds{1}_{x<\frac{t}{2}}(x)$ (1.26)

However, this is a stable one: if we apply a similar perturbation to the
initial data as previously:

$u_{\epsilon}(x)=\frac{\epsilon-x}{\epsilon}\mathds{1}_{0<x<\epsilon}(x)+\mathds{1}_{x<0}(x)$
(1.27)

then it is easy to understand that the effect of the perturbation vanishes
quickly and the evolution will continue as a shock. One can describe this
solution as physical in contrast to the non-physical previous one. Formal
criteria allowing to classify solutions according to their admissibility is
the subject of much literature, often called admissibility conditions, or
entropy conditions, we will revise some of them later.

####  1.1.5 Hopf’s treatment of the Burgers equation

The Burgers equation originates from fluid mechanics. It is the simplest non-
trivial example of a scalar conservation law. A more realistic version of it
includes the viscosity:

$u_{t}+uu_{x}=\nu u_{xx}$ (1.28)

Assume the initial data $u_{0}$ to be bounded. The viscous term expresses a
dependence of the flux on the gradient of $u$ in addition to its dependence on
$u$. In other words, it takes into account the diffusion. What Eberhard Hopf
showed in his paper [62] is that unlike the non-viscous Burgers equation, the
viscous one does not suffer from a lack of unicity of weak solutions, and it’s
only when we set $\nu=0$ where we can encounter this issue of multiple
solutions. Among these multiple solutions, only one is actually resulting from
taking the limit $\nu\rightarrow 0$. This sets a reasonable admissibility
condition for the non-viscous Burgers equation. Solutions verifying this
condition are often called viscous solutions, and the method is called: the
vanishing viscosity method. Because its mathematical and physical importance,
we will be synthesizing Hopf’s paper.

Let first $U$ be an integral of $u$ (sometimes called the hight function as a
reference to the growth process):

$U(x,t)=\int_{0}^{x}u(y,t)dy$ (1.29)

Then:

$U_{x}=u$ (1.30)

We can express Burgers equation in terms of $U$:

$U_{t}=\int_{0}^{x}u_{t}dy=\int_{0}^{x}uu_{x}-\nu
u_{xx}dy=\frac{1}{2}u^{2}-\nu u_{x}$ (1.31) $U_{t}+\frac{1}{2}U_{x}^{2}=\nu
U_{xx}$ (1.32)

We will see how this equation can be understood as the Hamilton-Jacobi
equation. Let’s consider the change of variable:

$\phi=e^{-\frac{U}{2\nu}}$ (1.33)

With elementary operations, we can show that (1.32) get reduced to the heat
equation:

$\phi_{t}=\nu\phi_{xx}$ (1.34)

The solution of this equation is unique and obtained by a convolution of the
initial condition: $\phi(x,0)=e^{\frac{-U(x,0)}{2\nu}}$ with the solution of a
Dicac initial condition:

$\begin{split}\phi(x,t)=(K_{0}(.,t)\ast\phi(.,0))(x)&=\frac{1}{2\sqrt{\pi\nu
t}}\int_{\mathbb{R}}\exp(-\frac{(x-y)^{2}}{{4\nu t}})\phi(y,0)dy\\\
&=\int_{\mathbb{R}}\exp(\frac{-G(x,y,t)}{2\nu})dy\end{split}$ (1.35)

where:

$K_{0}(x,t)=\frac{1}{2\sqrt{\pi\nu t}}e^{-\frac{x^{2}}{4\nu t}}$ (1.36)

and:

$G(x,y,t)=\frac{(x-y)^{2}}{2t}+U_{0}(y)$ (1.37)

We are interested in solving equation (1.28)

$u^{\nu}=-2\nu\frac{\phi_{x}}{\phi}=\frac{\int_{\mathbb{R}}\frac{(x-y)}{{t}}\exp(\frac{-G}{2\nu})dy}{\int_{\mathbb{R}}\exp(\frac{-G}{2\nu})dy}$
(1.38)
$=\frac{x}{t}-\frac{1}{t}\frac{\int_{\mathbb{R}}y\exp(\frac{-G}{2\nu})dy}{\int_{\mathbb{R}}\exp(\frac{-G}{2\nu})dy}$
(1.39)

We are looking for the limit $\lim_{\nu\rightarrow 0}u^{\nu}(x,t)$ If $G$ is
regular ($\mathcal{C}^{2}$) with a unique minimum, we can use the saddle point
method, by expanding $G$ to the second order in the neighborhood of the
minimum. We get:

$\lim_{\nu\rightarrow 0}u^{\nu}(x,t)=\frac{x-y^{0}(x,t)}{t}$

Where $y^{0}(x,t)$ is the point where $G(x,y,t)$ reaches its minimum for the
$y$ variable.

In case we have multiple points for which $G$ reaches its minimum then this
limit is not defined, and we need a couple of extra tools. Let’s denote
$y^{*}$ as the maximum and $y_{*}$ as the minimum of the set of points where
$G$ reaches its minimum. One can show the following:

  1. 1.

$x_{1}<x_{2}\implies y^{*}(x_{1},t)\leq y_{*}(x_{2},t)$

  2. 2.

$y_{*}(x-0,t)=y_{*}(x,t)$ and $y^{*}(x+0,t)=y^{*}(x,t)$

  3. 3.

$y_{*}(\infty,t)=\infty$ and $y^{*}(-\infty,t)=-\infty$

These properties can be proved rigorously, but we will instead contend with an
intuitive understanding. First let’s notice that $G(x,y,t)$ attains its
minimum at the same values as the function :

$G(0,y,t)-\frac{x}{t}y$ (1.40)

Now we can understand the previous properties with the help of figure 1.3

$y$$G(0,y,t)$$y_{*}(x,t)$$y^{*}(x,t)$slope $\frac{x}{t}$ Figure 1.3:
Illustration of the behavior of the set of minimums of $G$.

###### Theorem

$\frac{x-y^{*}(x,t)}{t}=\limsup_{\nu\rightarrow
0}u^{\nu}(x,t)\leq\liminf_{\nu\rightarrow
0}u^{\nu}(x,t)=\frac{x-y_{*}(x,t)}{t}$ (1.41)

###### Remarks

  * •

This formula is valid even when G is not smooth.

  * •

For every $x\in\mathbb{R}$ except on a countable set (forming the shocks) we
have

$\lim_{\nu\rightarrow
0}u^{\nu}(x,t)=\frac{x-y^{*}(x,t)}{t}=\frac{x-y_{*}(x,t)}{t}$ (1.42)

  * •

If $G(0,y,t)$ is convex then $u(.,t)$ is continuous.

###### A first example

Let’s apply this formalism to find the solution of The Burgers equation for in
initial step function:

$u_{0}(x)=\mathds{1}_{x>0}(x)$ (1.43)

So we need to find where (1.40) attains its minimum:

$G(0,y,t)-\frac{x}{t}y=\frac{y^{2}}{2t}+y\mathds{1}_{y>0}(y)-\frac{x}{t}y$
(1.44)

It is smooth everywhere except at $y=0$. If $x<0$ then the minimum is at
$y=x$. If $x>t$ then the minimum is $x-t$.

$y^{*}(x,t)=y_{*}(x,t)=\begin{cases}x\quad if\quad x<0\\\ x-t\quad if\quad
x>t\\\ 0\quad if\quad 0<x<t\\\ \end{cases}$ (1.45)

So the solution $\frac{x-y^{*}}{t}$ becomes

$u(x,t)=\begin{cases}1\quad if\quad x<0\\\ x-t\quad if\quad x>t\\\
\frac{x}{t}\quad if\quad 0<x<t\\\ \end{cases}$ (1.46)

Which is the expected solution.

###### A second example

Let’s consider the other classical initial profile:

$u_{0}(x)=\mathds{1}_{x<0}(x)$ (1.47)

In this case we have:

$G(0,y,t)=\frac{y^{2}}{2t}+y\mathds{1}_{y<0}(y)$ (1.48)

We can see from the figure that the locus of the minimum of
$G(0,y,y)-x\frac{y}{t}$ will not be continuous with respect to $x$. The
discontinuity can be found simply: $x=\frac{t}{2}$

$\begin{split}&y^{*}(\frac{1}{2},t)=1\\\ &y_{*}(\frac{1}{2},t)=0\end{split}$
(1.49) $y^{*}(x,t)=y_{*}(x,t)=\begin{cases}x\quad\quad if\quad
x<\frac{t}{2}\\\ x-t\quad if\quad x>\frac{t}{2}\\\ \end{cases}$ (1.50)

And finally the expected solution

$u(x,t)=\begin{cases}1\quad if\quad x<\frac{t}{2}\\\ 0\quad if\quad
x>\frac{t}{2}\\\ \end{cases}$ (1.51)

We notice that this formalism identifies and select naturally the stable
solution of Burgers equation among the weak ones. We will be visiting in what
follows more formal admissibility conditions that apply more generally.

####  1.1.6 Kružkov entropy condition

The work of Kružkov in the 70’s [63] represents a major development in the
understanding of conservation laws. Let $u$ be a smooth solution (so in the
strong sense) to the problem:

$\begin{split}&u_{t}+(f(u))_{x}=0\\\ &u(x,0)=u_{0}(x)\end{split}$ (1.52)

with $f\in\mathcal{C}^{1}$, Let $\eta$ be a positive convex $\mathcal{C}^{1}$
function. The claim is that $\eta(u)$ will be a conserved quantity under the
evolution of $u$. This is quite easy to show:

$\eta(u)_{t}=\eta^{\prime}(u)u_{t}=\eta^{\prime}(u)(-(f(u))_{x})=-\eta^{\prime}(u)f^{{}^{\prime}}(u)u_{x}$
(1.53)

So if we define $q$ so that:

$q^{{}^{\prime}}(u)=\eta^{\prime}(u)f^{{}^{\prime}}(u)$ (1.54)

We get:

$\eta(u)_{t}+q_{x}=0$ (1.55)

$\eta(.)$ is called an entropy function and $q$ is the associated entroy flux.

This property seems quite counter-intuitive at a first glance, especially
because there are few assumptions on $\eta$, but the key point is that it is
true only for the ”strong” evolution, and can actually be understood
geometrically with the help of the figure 1.4. This property allows us to
distinguish a strong solution from a weak solution, and thus will allow us to
establish admissibility criteria for selecting the viscous solution among the
weak ones, as we will see in what follows.

$u$$x$$f^{{}^{\prime}}(u)$$f^{{}^{\prime}}(u)$ (a) A regular profile

$u$$x$$f^{{}^{\prime}}(u)$$f^{{}^{\prime}}(u)$$\frac{f(u^{L})-f(u^{R})}{u^{L}-u^{R}}$
(b) A singular profile

Figure 1.4: For a regular profile(on the left), the front of an infinitesimal
slice moves at the same speed as its back. This property is conserved by the
application of $\eta$. On the other side, for a weak solution, the ”mass”
moves between different slices, this is responsible for the violation of the
conservation of $\eta(u)$

##### Entropic Admissibility condition

$u$ is said to be an entropy solution if for all entropy functions $\eta$ with
the corresponding flux $q$, this inequality is verified:

$\eta(u)_{t}+q_{x}\leq 0$ (1.56)

It’s fairly easy to show that this is a necessary condition for any viscous
solution: Consider a viscous solution $u=\lim_{\epsilon\rightarrow
0}u^{\epsilon}$, where $u^{\epsilon}$ verifies:

$u^{\epsilon}_{t}+f(u^{\epsilon})_{x}=\epsilon u^{\epsilon}_{xx}$ (1.57)

Then by multiplying both sides by $\eta^{{}^{\prime}}(u)$ we have:

$\eta^{{}^{\prime}}(u^{\epsilon})u^{\epsilon}_{t}+\eta^{{}^{\prime}}(u^{\epsilon})f^{\prime}(u^{\epsilon})u^{\epsilon}_{x}=\epsilon\eta^{{}^{\prime}}(u^{\epsilon})u^{\epsilon}_{xx}$
(1.58)

$\eta(u^{\epsilon})_{t}+q^{\prime}(u^{\epsilon})u^{\epsilon}_{x}=\epsilon(\eta(u^{\epsilon})_{xx}-\eta^{{}^{\prime\prime}}(u^{\epsilon})(u^{\epsilon}_{x})^{2})$
(1.59)

Since $\eta^{{}^{\prime\prime}}(u^{\epsilon})(u^{\epsilon}_{x})^{2}\geq 0$

$\eta(u^{\epsilon})_{t}+q^{\prime}(u^{\epsilon})u^{\epsilon}_{x}\leq\epsilon\eta(u^{\epsilon})_{xx}$
(1.60)

so in the limit $\epsilon\rightarrow 0$ we get:

$\eta(u)_{t}+q_{x}\leq 0$ (1.61)

This inequality has to be understood in the weak sens:

$\int_{\mathbb{R\times R^{+}}}(\eta(u)_{t}\phi(x,t)+q_{x}\phi(x,t))dxdt\leq 0$
(1.62)

For a positive test function $\phi$. Then by Green’s formula:

$\int_{\mathbb{R\times
R^{+}}}(\eta(u)\phi_{t}(x,t)+q(u)\phi_{x}(x,t))dxdt-\int_{\mathbb{R^{+}}}\eta(u)\phi(0,t)dt\geq
0$ (1.63)

###### Remarks

  1. 1.

Suppose the previous inequality is verified for $(\eta_{1},q_{1})$ and for
$(\eta_{2},q_{2})$ then it is verified for $\eta=\eta_{1}+\eta_{2}$ associated
with the flux $q=q_{1}+q_{2}$. Note however that the sum might not be convex.

  2. 2.

One needs not to check the previous inequality for all convex continuous
$\eta$, it’s enough to check it for a this special family:
$\\{\eta_{k}(u)=|u-k|\quad k\in\mathbb{R}\\}$. The associated flux of this
family is $q_{k}(u)=(f(u)-f(k))sgn(u-k)$, with $sgn$ is the sign function. The
entropic inequality then becomes:

$\int_{\mathbb{R\times
R^{+}}}(|u-k|\phi_{t}+(f(u)-f(k))sgn(u-k)\phi_{x})dxdt-\int_{\mathbb{R^{+}}}|u_{0}(x)-k|\phi(x,0)dx\geq
0$ (1.64)

Proof: It’s actually possible to establish a sequence of piece-wise affine
functions that converge to any continuous convex function such that each
element of the sequence is of the form:

$\eta^{n}=a^{n}+b^{n}u+\sum_{k}c^{n}_{k}|u-k|$ (1.65)

Each term of the sequence verifies the inequality thanks to the previous
remark, so is the limit.

  3. 3.

If $f$ is convex then it’s enough to check for one $\eta$

  4. 4.

Oleinik’s entropy condition If we consider a solution that is smooth
everywhere except for a discontinuity at $s(t)$ at time $t$, then we can apply
on $\eta(u)$ the same calculations as for the Hugnoiot-Rankine condition, and
it would lead to the inequality:

$(\eta(u^{L})-\eta(u^{R}))n^{t}+(q(u^{L})-q(u^{R}))n^{x}\geq 0$ (1.66)

if we apply this to the family $\eta_{k}(u)=|u-k|=(u-k)sgn(u-k)$, we get:

$\underset{u_{R}}{\overset{u^{L}}{[}}|u-k|]n^{t}+\underset{u_{R}}{\overset{u^{L}}{[}}(f(u)-f(k))sgn(u-k)]n^{x}\geq
0$ (1.67)

Let’s choose $k$ in between $u^{L}$ and $u^{R}$: $k=\lambda
u^{L}+(1-\lambda)u^{R}$ with $\lambda\in[0,1]$

$((u^{L}+u^{R}-2k)n^{t}+(f(u^{L})+f(u^{R})-2f(k))n^{x})sgn(u^{L}-u^{R})\geq 0$
(1.68)

Now using the Hugoniot-Rankine relation:

$(\frac{f(u^{L})-f(u^{R})}{u^{L}-u^{R}}(u^{L}+u^{R}-2k)+(f(u^{L})+f(u^{R})-2f(k)))sgn(u^{L}-u^{R})\geq
0$ (1.69)

noticing that $u^{L}+u^{R}-2k=(2\lambda-1)(u^{L}-u^{R})$ we get:

$((\lambda f(u^{L})+(1-\lambda)f(u^{R})-f(k)))sgn(u^{L}-u^{R})\geq 0$ (1.70)

Which means:

     * •

If $u^{L}>u^{R}$, then $f$ restricted to the $[u^{R},u^{L}]$ is under its
chord, figure 2.7(b)

     * •

If $u^{L}<u^{R}$, then $f$ restricted to the $[u^{L},u^{R}]$ is above its
chord , figure 2.7(a)

This constitutes Oleinik’s entropy condition for the admissibility of
discontinuous shocks.

  5. 5.

The stability interpretation:

One can understand Oleinik’s entropy condition intuitively in terms of the
stability of the shock: assume $u^{R}<u^{*}<u^{L}$, then the condition can be
written as: (figure 1.5 )

$\frac{f(u^{*})-f(u^{R})}{u^{*}-u^{R}}\leq\frac{f(u^{L})-f(u^{*})}{u^{L}-u^{*}}$
(1.71)

This means that if the shock between $u^{L}$ and $u^{R}$ is split (due to a
small perturbation) into two shocks: one between $u^{L}$ and an intermediate
value $u^{*}$ followed by one between $u^{*}$ and $u^{R}$, in order for the
two shocks to unite again, the first shock has to be faster than the second,
which is given by (1.71).See figure 1.6 for illustration. This condition is
referred to sometimes as Liu entropy condition.

$u$$f(u)$$u^{L}$$u^{R}$$u^{*}$ (a)

$u$$f(u)$$u^{L}$$u^{R}$$u^{*}$ (b)

Figure 1.5: The situations where the Oleinik condition is verified
$x$$u$$u^{L}$$u^{*}$$u^{R}$ Figure 1.6: Illustration of the stability
condition of a decreasing shock: the upper sub shock should have a higher
speed than the lower sub shock for any intermediate split value $u^{*}$

  6. 6.

If $f$ is convex, then the Oleinik’s entropy condition becomes particularly
simple: only decreasing shocks are admissible:

$u^{L}>u^{R}$ (1.72)

If $f$ is concave, the admissible shocks are the increasing ones.

  7. 7.

We can rewrite 1.71 slightly differently:

$\frac{f(u^{L})-f(u^{R})}{u^{L}-u^{R}}\leq\frac{f(u^{L})-f(u^{*})}{u^{L}-u^{*}}$
(1.73)

One can understand the equivalence between the two inequality easily by
contemplating figure 1.5. the relevance of this form is its adaptability to a
generalization to the non-scalar case, as we will see later.

We saw that the entropy condition is a necessary condition for a viscous
solution. It is possible to show that an entropic solution is unique. so, if
we admit the existence of a viscous solution, the entropic condition is
sufficient for selecting it. This result is known as Kružkov uniqueness
theorem. The proof relies on an $L^{1}$ contraction property of entropic
solutions. It can be found in chapter 3 of [54]

####  1.1.7 Relation to Hamilton-Jacobi equation

The objective here is to state one of the most classical formulas for scalar
conservation laws: Lax-Okeinik formula that describes an entropic solution for
arbitrary bounded initial condition. We choose to arrive from a path familiar
to most physicist: the Hamilton-Jacobi equation.

Let’s start by recalling the context of the HJE. Consider a system endowed
with a Lagrangian: $\mathcal{L}(q,\dot{q},t)$. We define the action as:

$S_{q_{0},t_{0}}(q,t):=min\int_{t_{0},q_{0}}^{t,q}\mathcal{L}(q,\dot{q},\tau)d\tau$
(1.74)

Where the minimum is taken over all the trajectories $q(.)$ such that
$q(t_{0})=q_{0}$ and $q(t_{1})=q_{1}$. Note that this definition is slightly
technically different from the usual definition of the action that is a
functional of the trajectories and not a function of
$(q,t)\in\mathbb{R}^{n}\times\mathbb{R}^{+}$. For what follows we set $n=1$
for simplicity. Many properties can be generalized trivially to higher
dimensions.

Let’s notice that the action verifies the following property: consider an
infinite path $\gamma$ that divides the the $(q,t)$ space into two parts, one
containing $(q_{0},t_{0})$ and the other $(q,t)$ then:

$S_{q_{0},t_{0}}(q,t):=\min_{(\tilde{q},\tilde{t})\in\Gamma}(S_{q_{0},t_{0}}(\tilde{q},\tilde{t})+S_{\tilde{q},\tilde{t}}(q,t))$
(1.75)

This suggest a different way to initialize the action. Consider that we know
the action over a path (for simplicity $\mathbb{R}\times\\{0\\}$ ):

$S(q,0)=S_{0}(q)$ (1.76)

then we can define the action all over the space by:

$S(q,t)=\min_{\tilde{q}}(S_{0}(\tilde{q})+S_{\tilde{q},0}(q,t))$ (1.77)

We will see a bit later the relevance of this definition. Let’s continue first
our development of the HJE. elementary calculus of variations of (1.74) allows
to establish:

$\frac{\partial S}{\partial q}=\frac{\partial\mathcal{L}}{\partial\dot{q}}:=p$
(1.78)

And equally:

$\frac{\partial S}{\partial t}=-\dot{q}p+\mathcal{L}:=-H$ (1.79)

Note that $q$ and $t$ are only spectators regarding the Legendre
transformation between $H$ and $\mathcal{L}$. After writing $\dot{q}$ in terms
of $p$,$q$ and $t$ and then $H$ in terms of $(q,p,t)$, one can use the
previous two equations to have a first order PDE for $S$:

$\frac{\partial S}{\partial t}+H(q,\frac{\partial S}{\partial q},t)=0$ (1.80)

This is the Hamilton Jacobi equation. Once $S$ is found, one can find the
trajectories in a similar fashion to finding the light rays in geometrical
optics out of the wavefronts. ( When $q$ is scalar, the image of the
trajectories is trivial, one needs only to find $\dot{q}$ as a function of
$q$, which can be done directly from (1.78) ), for more details: [64].
Actually, the characteristics (1.6) of this equation are Hamilton’s equations.
At this stage, it’s easy to make sense of the famous Hopf-Lax formula

###### Hopf Lax formula:

If we consider the HJE with the initial condition:

$S(q,t)=S_{0}(q)$ (1.81)

From the previous discussion, we can write $S$ as:

$S(x,t)=\inf\\{\int_{0}^{t}L(\dot{q}(s))ds+S_{0}(y)\quad|\quad
q_{0}=y,q(t)=x\\}$ (1.82)

Let’s now assume that $H$ is convex and superlinear
$\lim_{|p|\rightarrow\infty}\frac{H(p)}{p}=\infty$, then the Hopfs-Lax formula
tells us that we don’t actually need to minimize over all the trajectories
starting from the path and reaching the point$(x,t)$, it’s enough to minimize
among straight lines with the constant speed $\frac{x-y}{t}$ :

$S(x,t)=\min_{y\in\mathbb{R}}\\{t\mathcal{L}(\frac{x-y}{t})+S_{0}(y)\\}$
(1.83)

###### Proof.

If we choose a straight trajectory from point $y\in\mathbb{R}$ to $x$ with a
constant speed $\dot{q}=\frac{x-y}{t}$ It’s obvious that

$S(x,t)\leq\min_{y\in\mathbb{R}}\\{t\mathcal{L}(\frac{x-y}{t})+S_{0}(y)\\}$

to show the other direction inequality, we use the convexity of $\mathcal{L}$
by applying Jensen’s inequality:

$\mathcal{L}(\int_{0}^{t}\frac{1}{t}\dot{q}(s)ds)\leq\frac{1}{t}\int_{0}^{t}\mathcal{L}(\dot{q}(s))ds$
(1.84)

So that makes:

$t\mathcal{L}(\frac{x-y}{t})\leq\int_{0}^{t}\mathcal{L}(\dot{q}(s))ds$ (1.85)

Which completes the proof. ∎

###### Regularity and uniqueness of the solution

The Hopf-Lax formula provides a weak form solution for the HJE as it doesn’t
require for $S$ to be differentiable. Let’s be more precise: assumes $S_{0}$
to be Lipschitz with $Lip(S_{0})\leq M$ then $S$ defined by Hopf-Lax solves
the HJE a.e., and it is Lipschitz with $Lip(S_{t})\leq M$ and differentiable
almost everywhere.

Now we reach the stage where we can show the link to the conservation law. If
we derive (1.80) with respect to $q$, and ignore the dependence of $H$ on $q$
and $t$ we get

$\frac{\partial}{\partial t}\frac{\partial S}{\partial q}+\frac{\partial
H}{\partial q}(\frac{\partial S}{\partial q})=0$ (1.86)

Now we can identify $u$ with $\frac{\partial S}{\partial q}$ and $H$ with $f$
and we find our conservation law.

###### Lax-Oleinik formula

The solution $S$ defined by Lax-Hopf is differentiable almost everywhere
(Rademacher’s theorem). We would like to work out its derivative:

$u(x,t):=\frac{\partial}{\partial
x}\min_{y\in\mathbb{R}}\\{t\mathcal{L}(\frac{x-y}{t})+S_{0}(y)\\}$ (1.87)

Under some assumptions:

  * •

$f(0)=0$

  * •

$f$ is uniformly convex

  * •

$f$ is smooth

  * •

$u_{0}$ is bounded

And let $G(x)=(f^{\prime})^{-1}(x)$, then

  * •

There exists for almost all values of $x$, a unique $y(x,t)$ such that the
minimum is attained:

$\min_{y\in\mathbb{R}}\\{t\mathcal{L}(\frac{x-y}{t})+S_{0}(y)\\}=t\mathcal{L}(\frac{x-y(x,t)}{t})+S_{0}(y(x,t))$
(1.88)

  * •

$x\rightarrow y(x,t)$ is non decreasing

  * •

Almost everywhere for $x$ the previous derivative is:

$u(x,t)=G(\frac{x-y(x,t)}{t})$ (1.89)

This generalizes the Hopfs treatment of Burgers equation previously
encountered.

####  1.1.8 Riemann problem

The Riemann problem is a conservation system with constant initial data except
at zero:

$u_{0}(x)=\begin{cases}u^{L}&\text{if $x<0$}\\\ u^{R}&\text{if
$x>0$}\end{cases}$ (1.90)

The advantage of this initial condition is that it allows for a solution which
is invariant under the resealing:

$u(x,t)=u(\lambda x,\lambda t)$ (1.91)

In other words, the solution is a function of $\frac{x}{t}:=\xi$

$u(x,t)=u(\frac{x}{t}):=u(\xi)$ (1.92)

We can write the conservation equation as:

$u^{{}^{\prime}}(\xi)(f^{{}^{\prime}}(u(\xi))-\xi)=0$ (1.93)

This equation in the strong sense can give us insight into the regular
solutions: they can be either constants or of the form:
$u(\xi)=(f^{{}^{\prime}})^{-1}(\xi)$. This requires $f^{{}^{\prime}}$ to be
invertible. We need to take into account the shocks and the initial condition
and to treat the case where $f^{{}^{\prime}}$ is not invertible. It’s
convenient to start the discussion with a convex (or concave), then move to
the general form of flux.

##### Convex flux

If $u^{L}>u^{R}$ then there is a simple solution which is a shock at a speed
$\frac{f(u^{L})-f(u^{R})}{u^{L}-u^{R}}$ and it is an entropic solution since
it verifies the Oleinik condition.

If $u^{L}<u^{R}$, then we can have a continuous entropic solution:

$u(\xi)=\begin{cases}u^{L}&\text{if $\xi<f^{{}^{\prime}}(u^{L})$}\\\
(f^{{}^{\prime}})^{-1}(\xi)&\text{if
$f^{{}^{\prime}}(u^{L})<\xi<f^{{}^{\prime}}(u^{R})$}\\\ u^{R}&\text{if
$\xi>f^{{}^{\prime}}(u^{R})$}\end{cases}$ (1.94)

Thanks to the convexity, $f^{{}^{\prime}}$ is increasing and thus invertible.
The part the solution on the interval
$[f^{{}^{\prime}}(u^{L}),f^{{}^{\prime}}(u^{R})]$ is called a rarefaction fan

Note that if $f$ is not differentiable (but only continuous) at some point
$\tilde{u}$ then the solution is constant on the interval
$[f^{{}^{\prime}}(\tilde{u}-0^{+}),f^{{}^{\prime}}(\tilde{u}+0^{+})]$, so we
can have multiples rarefaction fans. For a convex flux, one cannot observe at
the same time a shock and a rarefaction fan.

The case of a concave flux is treated in exactly similar fashion except that
the shock will appear now when $u^{L}>u^{R}$ while the rarefaction for
$u^{L}<u^{R}$. We have actually the symmetry $f(u)\leftrightarrow-f(-u)$ that
allows to passe from one case to the other.

##### Non convex flux

The general non-convex, non-concave flux case has been treated by S.Osheri in
1983 [65] According to him the solution is:

if $u^{L}<u^{R}$:

$\xi u(\xi)-f(u(\xi))=\max_{v\in[u^{L},u^{R}]}\\{\xi v-f(v)\\}$ (1.95)

if $u^{L}>u^{R}$:

$\xi u(\xi)-f(u(\xi))=\min_{v\in[u^{L},u^{R}]}\\{\xi v-f(v)\\}$ (1.96)

There is an equivalent very simple formation, figure 1.7 (that I astonishingly
haven’t encountered it in the literature):

if $u^{L}<u^{R}$, we replace $f$ by its convex hull:

$\check{f}=\max\\{g\in\mathcal{C}^{0};g\leq f\\}$ (1.97)

if $u^{L}>u^{R}$, we replace $f$ by its concave hull:

$\hat{f}=\min\\{g\in\mathcal{C}^{0};g\geq f\\}$ (1.98)

One can understand the equivalence of the two formulations with the help of
some elementary geometrical constructions.

###### Remarks

  * •

if part of the flux is linear then this part corresponds to a discontinuity
moving at a speed equal to the slope of straight line, which means that the
Oleinik condition is a particular case of this formulation since the shock for
convex $f$ and $u^{L}>u^{R}$ can be as well seen as the solution of a linear
flux in the interval $[u^{R},u^{L}]$. This flux is called ”contact flux” in
the literature.

  * •

If the initial condition is ”Riemann like”, in the sense that it’s uniform on
the left and on the right, except on some bounded interval, then the re-scaled
solution converges to the limit shape of the corresponding Riemann problem.

$u$$f(u)$$u^{L}$$u^{R}$ (a) The concave hull (in red) of flux (in green),
needed when $u^{R}<u^{L}$

$u$$f(u)$$u^{R}$$u^{L}$ (b) The convex hull (in red) of flux (in green),
needed when $u^{R}>u^{L}$

Figure 1.7: Non convex/concave flux treatment, equivalent of Osheri’s
solution.

###  1.2 Hyperbolic Systems of Conservation Laws

In the first part we treated the case of a single conserved quantity, the
problem becomes significantly more complex with $n$ conserved quantities
$\mathbf{u}=(u_{1},...,u_{n})^{t}$ with a flux for the quantity $i$ that
depends on all of the quantities:
$\mathbf{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. The conservation law
becomes a system of $n$ coupled PDE’s:

$\begin{split}&\mathbf{u}_{t}+(\mathbf{f}(\mathbf{u}))_{x}=0\\\
&\mathbf{u}(x,0)=\mathbf{u}_{0}(x)\end{split}$ (1.99)

This system is said to be strictly hyperbolic when the differential of the
flux $A(u)=D_{u}f$ is diagonalisable in $\mathbb{R}$ and its eigenvalues are
distinct for all $\mathbf{u}$:

$\lambda_{1}<\lambda_{2}<..<\lambda_{n}$ (1.100)

This allows to choose the left and the right eigenvectors ($\mathbf{l}_{i}$
and $\mathbf{r}_{i}$ respectively) such that:

$\mathbf{l}_{i}.\mathbf{r}_{j}=\delta_{i,j}$ (1.101)

The proof of this is elementary:

$\mathbf{l}_{i}^{t}A\mathbf{r}_{j}=\mathbf{l}_{i}^{t}\lambda_{i}\mathbf{r}_{j}=\mathbf{l}_{i}^{t}\lambda_{j}\mathbf{r}_{j}$
(1.102)

Finally:

$(\lambda_{i}-\lambda_{j})\mathbf{l}_{i}.\mathbf{r}_{j}=0$ (1.103)

Before treating the general non-linear system, let’s have a look at the simple
case of a linear one:

####  1.2.1 A linear system

A simple situation when is $\mathbf{f}$ is linear:
$\mathbf{f}(\mathbf{u})=A\mathbf{u}$ The conservation system is:

$\mathbf{u}_{t}+A\mathbf{u}_{x}=0$ (1.104)

We can write the vector $\mathbf{u}$ as :

$\mathbf{u}=\sum_{i}(\mathbf{u}.\mathbf{l_{i}})\mathbf{r_{i}}$ (1.105)

Lets define $n$ new quantities: $\mathbf{\tilde{u}}=(\tilde{u}_{i})_{1\leq
i\leq n}$

$\tilde{u}_{i}=\mathbf{u}.\mathbf{l_{i}}$ (1.106)

Which are the densities in the base of the right eigenvectors. Then we can
realize by multiplying both sides of equation (1.104) by $\mathbf{l_{i}}$ that
each of these quantities verifies a scalar conservation law:

$(\tilde{u}_{i})_{t}+\lambda_{i}(\tilde{u}_{i})_{x}=0$ (1.107)

Where $\lambda_{i}$ is the eigenvalue associated with $\mathbf{l_{i}}$. So,
the new quantities evolve independently, and the solution of the original
system is simply:

$\mathbf{u}(x,t)=\sum_{i}(\mathbf{u}(x-\lambda_{i}t,0).\mathbf{l_{i}})\mathbf{r_{i}}$
(1.108)

The situation becomes significantly more complex when the matrix $A$ is a
function of $\mathbf{u}$. The different waves can now interact with each
other. We will consider only the situation when the system can be written in a
conservative form, i.e. when $A$ is the differential of a flux function. We
assume as well the strict hyperbolic condition as previously.

####  1.2.2 Weak solutions, the Rankine-Hugoniot condition

Similarly to the scalar case, one has to interpret the conservation equation
in the sense of distributions. This allows for discontinuous solutions to
exist. A straightforward generalization of the Rankine-Hugoniot condition is
possible: The discontinuities should verify:

$(\mathbf{u}^{L}-\mathbf{u}^{R})\sigma=(\mathbf{f}(\mathbf{u}^{L})-\mathbf{f}(\mathbf{u}^{R}))$
(1.109)

Unlike its scalar counterpart, this condition doesn’t only provide the shock
speed, but also constraint the densities between which one can have a shock
solution, namely for all $1\leq i,j\leq n$

$Det\begin{pmatrix}u_{i}^{L}-u_{i}^{R}&f_{i}(\mathbf{u}^{L})-f_{i}(\mathbf{u}^{R})\\\
u_{j}^{L}-u_{j}^{R}&f_{j}(\mathbf{u}^{L})-f_{j}(\mathbf{u}^{R})\end{pmatrix}=0$
(1.110)

####  1.2.3 The shock curves

Let’s fix a point in the density space $\mathbf{u^{*}}$ and search for all the
points that can be connected to $\mathbf{u^{*}}$ through a shock. One can see
this set as parameterized by $\sigma$, so it is expected to form a 1d
manifold, i.e a curve, but actually, it is composed of $n$ curves that passes
by $\mathbf{u^{*}}$. To see this, one can linearize 1.109 in the neighborhood
of $\mathbf{u^{*}}$:

$\mathbf{u}=\mathbf{u^{*}}+\frac{1}{\sigma}A(\mathbf{u^{*}})(\mathbf{u}-\mathbf{u^{*}})$
(1.111)

Since $A$ has $n$ real distinct eigenvalues, we can conclude that this
equation admits $n$ independent 1d eigenspaces that would represent the
tangents of $n$ shock curves. We can parameterize each by its speed:
$S^{i}_{\mathbf{u^{*}}}(\sigma)$, figure 1.8. Obviously, the small
perturbations in the i-shocks propagate at a speed $\lambda_{i}$:

$\lim_{\sigma\rightarrow\lambda_{i}}S^{i}_{\mathbf{u^{*}}}(\sigma)=\mathbf{u^{*}}$
(1.112)

Note that the $i$-shock curve emerging from one point does not coincide in
general with the $i$-shock curve emerging from another point located at the
former.

$\mathbf{u^{*}}$$\mathbf{r_{0}}$$\mathbf{r_{1}}$$S^{1}_{\mathbf{u^{*}}}$$S^{0}_{\mathbf{u^{*}}}$
Figure 1.8: Shock curves in a 2d density space

####  1.2.4 Admissibility conditions

Since weak solutions are not necessarily unique, we need to select among them
the ”physical” ones. Conceptually, we can add an infinitesimal diffusion term
to the conservation equation:

$\mathbf{u}_{t}+(\mathbf{f}(\mathbf{u}))_{x}+\epsilon\mathbf{u}_{xx}=0$
(1.113)

And search for a solution that are a limit in $L^{{}^{1}}_{loc}$ when
$\epsilon\rightarrow 0$.

Although, it was possible for the Burgers equation to be treated in this
manner as Hopf did, this approach doesn’t provide a practical procedure
allowing to eliminate non-physical solutions. One has to look for alternative
approaches.

##### Entropy condition

The notion of Entropy can be extended to the multi-dimensional case; however,
its existence is no longer guaranteed. Let’s recall that entropy is a smooth
convex scalar function, associated with a flux that verifies:

$D_{\mathbf{u}}q=D_{\mathbf{u}}\eta A(\mathbf{u})$ (1.114)

This is a system of $n$ first-order PDE with two scalar variables. For $n>2$
the system is over-determined and doesn’t in general have a solution. If it
does, then it allows for classifying solutions within three categories:

  * •

regular solutions (in the sense $C^{1}$) conserve the entropy under time
evolution.

  * •

physical singular solutions consume entropy

  * •

non-physical solutions produce entropy.

The second category means that a $\mathbf{u^{L}}-\mathbf{u^{R}}$
discontinuity, traveling at a speed $\lambda$, should verify:

$(\eta(\mathbf{u^{L}})-\eta(\mathbf{u^{R}}))\lambda\leq
q(\mathbf{u^{L}})-q(\mathbf{u^{R}})$ (1.115)

It’s not possible to extend the Oleinik condition in a straightforward way
even if we restrict ourselves to a particular i-shock curve. For this to
happen, the stability condition has to be formulated in the sense of Liu.

##### Liu Condition

Consider an i-shock curve that originates at $\mathbf{u^{L}}$ and let
$\mathbf{u^{R}}$ be a point that belongs to that curve:
$\mathbf{u^{R}}=S^{i}_{\mathbf{u^{L}}}(\sigma)$ and let $\mathbf{u^{*}}$ be an
intermediate point on the curve between $\mathbf{u^{L}}$ and $\mathbf{u^{R}}$:
$\mathbf{u^{*}}=S^{i}_{\mathbf{u^{L}}}(s)$. The Liu stability condition states
that:

$\sigma(\mathbf{u^{L}},\mathbf{u^{R}})\leq\sigma(\mathbf{u^{L}},\mathbf{u^{*}})$
(1.116)

This obviously means that the intermediate perturbation shock will not form a
separate shock but rather join back with the mother shock. This follows the
same logic as the scalar case, figure 1.6, restricted on one i-shock curve.
This condition was developed by Liu in his paper: [66]. Another handy directly
applicable condition is the Lax condition, introduced in the next paragraph.

##### Lax condition

Let $\lambda_{i}(\mathbf{u^{L}})$,$\lambda_{i}(\mathbf{u^{R}})$ be the
i-eigenvalues at $\mathbf{u^{L}}$, $\mathbf{u^{R}}$ respectively, then Lax
stability condition can be expressed as:

$\lambda_{i}(\mathbf{u^{L}})\geq\sigma(\mathbf{u^{L}},\mathbf{u^{R}})\geq\lambda_{i}(\mathbf{u^{R}})$
(1.117)

This means that the small perturbations on the left and on the right of the
shock should move towards the shock. In terms of characteristics: In the
neighborhood of the i-shock, the neighboring i-characteristics should be
entering the shock and not leaving it, figure 1.9.

$x$$t$Admissible (a) Characteristics are joining the shock

$x$$t$Not admissible (b) Characteristics are leaving the shock

Figure 1.9: Illustration of Lax condition in terms of the behavior of of
characteristic around the shock represented by the brown segment

In a sense, this condition represents the time-irreversible character of the
singular solution. The information of the initial data is lost at the shocks.
mathematical solutions that inverse this process are not physical.

##### Conclusion

Concretely, the admissibility conditions will eliminate for each shock curve
originating from $\mathbf{u^{*}}$, one of the two parts of the shock curve
separated by $\mathbf{u^{*}}$. Let’s denote $S^{i+}_{\mathbf{u^{*}}}$ the non-
admissible part and $S^{i-}_{\mathbf{u^{*}}}$ the admissible one.

We orientate $\mathbf{r}_{i}$ to the non-admissible part, figure 1.10. This
orientation will be compatible with further developments.

$\mathbf{u^{*}}$$\mathbf{r_{0}}$$\mathbf{r_{1}}$$S^{1}_{\mathbf{u^{*}}}$$S^{0}_{\mathbf{u^{*}}}$
Figure 1.10: Admissible sections of shock curves(green). Non-admissible ones
are colored(red)

Recall that in the case of a scalar system, the convexity of the current
allowed a simplification of the analysis in particular because the mapping
from the density to the speed of perturbations becomes monotonous, and thus
one to one. In higher dimensions, we need this character to be conserved on
the integral curves the eigenvectors of the Jacobian matrix as we will see in
what follows.

##### A simplifying hypothesis

A typical hypothesis in the textbooks that has its origins to Lax 1957 [67] is
to assume that each of the eigenvectors’ fields falls into one of the two
categories:

  * •

$D_{\mathbf{u}}\lambda_{i}.\mathbf{r}_{i}(\mathbf{u})>0$ for all $\mathbf{u}$
and in this case it’s said to be genuinely non linear

  * •

$D_{\mathbf{u}}\lambda_{i}.\mathbf{r}_{i}(\mathbf{u})=0$ for all $\mathbf{u}$
and is said to be linearly degenerate.

The first case means that the directional derivative of $\lambda_{i}$ in the
direction of $\mathbf{r}_{i}$ is positive. Obviously, the sign is irrelevant
as far as it doesn’t change. However, for later convenience, we choose the
direction of $\mathbf{r}_{i}$ so that this sign is positive. This means that
$\lambda_{i}$ is an increasing function along the directed integral curves of
the i-field. The second case simply means that $\lambda_{i}$ is constant along
these curves.

####  1.2.5 Rarefaction Curves

For a genuinely non-linear field, we obviously can parameterize an integral
curve by the corresponding eigenvalue field. So i-curve passing by a density
point $\mathbf{u^{*}}$ can be described by:

$\lambda_{i}\rightarrow R_{\mathbf{u^{*}}}^{i}(\lambda_{i})$ (1.118)

In other words, this is the solution (and the flow for the parameter
$\mathbf{u^{*}}$) of the ODE:

$\frac{d\mathbf{u}}{ds}=\frac{\mathbf{r}_{i}(\mathbf{u})}{D_{\mathbf{u}}\lambda_{i}.\mathbf{r}_{i}(\mathbf{u})}$
(1.119)

This curve is called the i-rarefaction curve.

The i-rarefaction curve that passes by $\mathbf{u^{*}}$ is tangent to the
i-shock curve that originates at $\mathbf{u^{*}}$. It is even possible to show
that they have the same curvature, figure 1.11. In a class of conservation
laws known as the Temple class, the two curves are identical for all the
fields.

$\mathbf{u^{*}}$$\mathbf{r_{0}}$$S^{0}_{\mathbf{u^{*}}}$$R^{0}_{\mathbf{u^{*}}}$
Figure 1.11: A rarefaction curve with a shock curve

$\mathbf{u^{*}}$ divides the rarefaction curve passing by it into two parts.
We denote $R^{i+}_{\mathbf{u^{*}}}$ the part verifying
$\lambda_{i}(\mathbf{u})>\lambda_{i}(\mathbf{u^{*}})$, and by
$R^{i-}_{\mathbf{u^{*}}}$, the other part.

####  1.2.6 T-curves

It will soon become meaningful to define a new curve by sticking
$R^{i+}_{\mathbf{u^{*}}}$ to $S^{i-}_{\mathbf{u^{*}}}$. This curve is called a
T-curve.

$S^{i-}_{\mathbf{u^{*}}}$$\mathbf{r_{i}}$$\mathbf{u^{*}}$$R^{i+}_{\mathbf{u^{*}}}$
Figure 1.12: The structure of a T-curve

In the next paragraph, we will be considering the solutions of the system of
conservation laws for the particular Riemann initial condition.

####  1.2.7 The Riemann Problem

We consider the system of conservation laws associated with an initial
condition:

$\mathbf{u}(x,0)=\mathbf{u^{L}}\mathds{1}_{x<0}+\mathbf{u^{R}}\mathds{1}_{x>0}$
(1.120)

It is easy to verify that the system is invariant under the transformation
$(x,t)\rightarrow(\mu x,\mu t)$, so one can express the solution in terms of
the variable $\xi=\frac{x}{t}$ and we can convert the system into an equation
of $\mathbf{u}(\xi)$:

$\mathbf{u}^{{}^{\prime}}(\xi)(\xi-A)=0$ (1.121)

Besides trivial constant solutions we can identify the following ones:

##### Elementary solutions

We can distinguish the following simple solutions composed of a single type of
waves:

  * •

If $\mathbf{u^{L}}$ and $\mathbf{u^{R}}$ belong to the same i- shock curve and
verify: $\lambda_{i}(\mathbf{u^{L}})\geq\lambda_{i}(\mathbf{u^{R}})$ then the
solution of the Riemann problem is a simple shock:

$\mathbf{u}(\xi)=\begin{cases}\mathbf{u^{L}}&\text{if
}\xi<\sigma(\mathbf{u^{L}},\mathbf{u^{R}})\\\ \mathbf{u^{R}}&\text{if
}\xi>\sigma(\mathbf{u^{L}},\mathbf{u^{R}})\end{cases}$ (1.122)

  * •

If $\mathbf{u^{L}}$ and $\mathbf{u^{R}}$ belong to the same genuinely non-
linear i-rarefaction curve and verify:
$\lambda_{i}(\mathbf{u^{L}})<\lambda_{i}(\mathbf{u^{R}})$ then the solution of
the Riemann problem is a simple Rarefaction wave:

$\mathbf{u}(\xi)=\begin{cases}\mathbf{u^{L}}&\text{if
}\xi<\lambda_{i}(\mathbf{u^{L}})\\\ R_{\mathbf{u^{L}}}^{i}(\xi)&\text{if
}\lambda_{i}(\mathbf{u^{L}})<\xi<\lambda_{i}(\mathbf{u^{R}})\\\
\mathbf{u^{R}}&\text{if }\xi>\lambda_{i}(\mathbf{u^{R}})\end{cases}$ (1.123)

To verify why the branch in the middle holds, it’s enough to multiply the
equation 1.121 from the right by $\mathbf{r_{i}}$, we get a factor with the
eigenvalue equation that is solved by this branch.

  * •

If $\mathbf{u^{L}}$ and $\mathbf{u^{R}}$ belong to the same linearly
degenerate i- rarefaction curve with a constant eigenvalue $\lambda_{i}$ then
the solution of the Riemann problem is a shock:

$\mathbf{u}(\xi)=\begin{cases}\mathbf{u^{L}}&\text{if }\xi<\lambda_{i}\\\
\mathbf{u^{R}}&\text{if }\xi>\lambda_{i}\end{cases}$ (1.124)

This type of shock is sometimes referred to as contact discontinuity. Unlike
the two previous cases, the linearly degenerate curve is bi-directional, it’s
at the same time a rarefaction and a shock curve.

In conclusion, to have a physical elementary solution, $\mathbf{u^{R}}$ has to
belong to one of the T-curves generated by $\mathbf{u^{L}}$.

##### Combined solutions

If $\mathbf{u^{L}}$ and $\mathbf{u^{R}}$ don’t belong to the same curve, one
has to combine different T-curves to connect the two. It’s possible to show
the existence and unicity of this combination of curves for $\mathbf{u^{L}}$
and $\mathbf{u^{R}}$ sufficiently close.

####  1.2.8 Riemann Variables

Let’s contemplate the vector field of the right eigenvectors:

$\mathbf{l_{i}}A=\lambda_{i}\mathbf{l_{i}}$ (1.125)

Let’s assume that up to a multiplicative scalar field $\mathbf{l_{i}}$ can be
derived from a scalar potential $z_{i}(\mathbf{u})$. i.e:
$\mathbf{l_{i}}\propto\nabla z_{i}$. We can choose the norm of
$\mathbf{l_{i}}$ so that we have equality:

$\mathbf{l_{i}}=\nabla z_{i}$ (1.126)

We call these scalar fields, the Riemann variables. They always exist for
$n=2$. For $n>2$ they don’t exist in general. Riemann variables simplify the
analysis of the conservation laws. They allow to partially decouple the
system. More precisely We have:

$\begin{split}\partial_{t}z_{i}(\mathbf{u})+\lambda_{i}\partial_{x}z_{i}(\mathbf{u})&=\nabla
z_{i}\partial_{t}(\mathbf{u})+\lambda_{i}\nabla
z_{i}\partial_{x}(\mathbf{u})\\\ &=\nabla
z_{i}(\partial_{t}(\mathbf{u})+\lambda_{i}\partial_{x}(\mathbf{u}))\\\
&=\mathbf{l_{i}}(\partial_{t}(\mathbf{u})+\lambda_{i}\partial_{x}(\mathbf{u}))\\\
&=\mathbf{l_{i}}(\partial_{t}(\mathbf{u})+A\partial_{x}(\mathbf{u}))=0\end{split}$
(1.127)

This means that if we express the system in terms of the Riemann variables,
the coupling between the equations appears only in the velocity coefficient
$\lambda_{i}(\mathbf{z})$.

The sets where $z_{i}$ is constant form a foliation of manifolds of dimension
$n-1$ in the density space, that are perpendicular to the integral curves of
$\mathbf{l_{i}}$ and since we have,
$\mathbf{l_{i}}.\mathbf{r_{j}}=\delta_{i,j}$, This means that $z_{i}$ is
constant over all the j-rarefaction curves such that $j\neq i$. For this
reason, they are sometimes called the Riemann invariants.

One can show that these speeds can be obtained by:

$\lambda_{i}(\bm{z})=\frac{\partial J_{k}}{\partial z_{i}}/\frac{\partial
u_{k}}{\partial z_{i}}$ (1.128)

This is true for any $k$ and it implies that:

$\frac{\partial J_{n}}{\partial z_{i}}\frac{\partial u_{m}}{\partial
z_{i}}=\frac{\partial J_{m}}{\partial z_{i}}\frac{\partial u_{n}}{\partial
z_{i}}$ (1.129)

##### Remark

The condition for the existence of Riemann variables is provided Frobenius
theorem. It’s more convenient to express it in the language of differentiable
forms. Let’s see $\mathbf{l_{i}}$ as a 1-form. If they exist, then we have a
scalar field $f$ such that $\mathbf{l_{i}}=fdz_{i}$. We have:
$d\mathbf{l_{i}}=df\wedge dz_{i}$. Now since $fdz_{i}\wedge df\wedge
dz_{i}=0$, we have the Frobenius condition:

$\mathbf{l_{i}}\wedge d\mathbf{l_{i}}=0$

It’s easy to see that in 2D, this is always verified

####  1.2.9 Temple Class Systems

A particular case of conservation laws for which explicit calculations are
typically possible and simpler is the Temple class. It was first noticed and
studied by Blake Temple 1982 [68]. I will be simplifying the main ideas of
this paper.

##### Definition

We say that a system is of temple class if for all $i$ the i-rarefaction curve
coincides with the i-shock curve. This of course happens trivially for an
i-field in the case of a contact discontinuity, which means that $\lambda_{i}$
is constant over each of the integral curves. i.e the i-field is linearly
degenerate. We will assume that none of the fields in this situation for what
follows.

##### Important Property

A system is of Temple class if and only if all the rarefaction curves and the
shock. curves are affine.

One of the directions of the implications is trivial. Since we know that the
i-rarefaction and the i-shock curve originated from a point are tangent, being
affine implies coinciding. For the other direction, I will not provide
rigorous proof, it can be found in the paper of Temple, but I can provide an
intuitive understanding: Since the shock curve coincides with the rarefaction
curve, it means that the shock curve in this case does not depend on the point
where it is originated from along the rarefaction curve. so if we take three
distinct points on this curve: $\mathbf{u_{1}}$, $\mathbf{u_{2}}$,
$\mathbf{u_{3}}$. we can write three Hugoniot conditions:

$f(\mathbf{u_{1}})-f(\mathbf{u_{2}})=\sigma_{1}(\mathbf{u_{1}}-\mathbf{u_{2}})$
(1.130)
$f(\mathbf{u_{2}})-f(\mathbf{u_{3}})=\sigma_{2}(\mathbf{u_{2}}-\mathbf{u_{3}})$
(1.131)
$f(\mathbf{u_{1}})-f(\mathbf{u_{3}})=\sigma_{3}(\mathbf{u_{1}}-\mathbf{u_{3}})$
(1.132)

If the curve is linearly degenerate, then $\sigma_{1}=\sigma_{2}=\sigma_{3}$.
However, we excluded this possibility. so the sigmas are not all identical
(except potentially at a subset of the curve of Lebesgue measure zero, at
which we can prolong the arguments by continuity). Obviously, the left side of
the third equation is the sum of the first two, so the same must hold for the
right side:

$\sigma_{1}(\mathbf{u_{1}}-\mathbf{u_{2}})+\sigma_{2}(\mathbf{u_{2}}-\mathbf{u_{3}})=\sigma_{3}(\mathbf{u_{1}}-\mathbf{u_{3}})$
(1.133)

It is easy to see now that if the sigmas are not all equal, then they are all
different. This implies a linear relationship between the three points and
means that they are aligned.

##### Remarks

  * •

If the initial data of our problem belong to a single rarefaction-shock curve,
then it will evolve staying on this curve for any time, so this curve is
sometimes called an invariant manifold of dimension one. In other words, it is
possible to restrict the system of conservation laws to the manifold, and the
restriction will be a scalar conservation law:

$\partial_{t}\rho+\lambda(\phi(\rho))\partial_{x}\rho=0$ (1.134)

where $\phi:\mathbb{R}\rightarrow\mathbb{R}^{n}$ is a parameterization of the
curve. This similarity between Temple class and scalar conservation laws
allows extending some of the general theories that were proven only for scalar
systems to Temple classes. For instance, the existence and uniqueness of
physical solutions were proven in [41]

  * •

Riemann variables are build-in within Temple systems: Let
$\mathbf{l}_{i}(\mathbf{u}^{*})$ be the left eigenvector at $\mathbf{u}^{*}$.
The j-rarefaction lines for all $j\neq i$ form a hyperplane perpendicular to
$\mathbf{l}_{i}(\mathbf{u})$ will be perpendicular to this hyperplane for all
$\mathbf{u}$ belonging to it. If we parameterize the family of hyperplanes
associated with the i-left field by $z_{i}$ then this will constitute
obviously a Riemann variable: $\mathbf{l_{i}}\propto\nabla z_{i}$. For a more
rigorous treatment of the existence of Riemann variables, one has to make use
of the Frobenius theorem.

  * •

A conservation system can of course be partially of Temple class, in the sense
that the properties of Temple apply only to some of the characteristic fields.
These fields would form an invariant manifold where initial data stay on it
and do not leave it. The restriction of the conservation system to this
manifold would be a Temple system.

##### Temple class for a system of two conservation laws.

Consider two conservation laws corresponding to non-degenerate fields.

$\begin{split}\partial_{t}u+\partial_{x}f(u,v)=0\\\
\partial_{t}u+\partial_{x}g(u,v)=0\end{split}$ (1.135)

We assume that the rarefaction lines are not parallel for either of the
fields. If one of them is, then this one decouples trivially from the system.
We can parameterize the i-rarefaction field by the slope of its lines, say
$z_{i}(u,v)$. These variables would obviously be Riemann invariant. We are
interested here in determining the class of currents that can generate such
systems. One can show that:

$\begin{split}g=fz_{0}+H_{1}(z_{0})\\\ g=fz_{1}+H_{2}(z_{1})\end{split}$
(1.136)

This means that if we know the currents on the boundaries of a temple class
system, we can determine the currents on the bulk.

##  Chapter 2 Hydrodynamic behavior of the two–TASEP

The content of this chapter is based on the article Cantini, Zahra (2022) that
is published in Journal of Physics A: Mathematical and Theoretical 55.30
(2022) https://iopscience.iop.org/article/10.1088/1751-8121/ac79e3/meta

### Abstract

We address the question of the hydrodynamic behavior of a 2-species
generalization of the TASEP, called 2–TASEP, introduced by Derrida [32] and
Mallick [33]. We find that the auxiliary variables, introduced previously in
the literature to express the density dependence of particle currents, turn
out to be the Riemann variables of the conservation equations. This allows us
to work out quite explicitly the rarefaction and shock solutions and to
completely solve the associated Riemann problem. Our theoretical results are
confirmed by Monte Carlo simulations.

###  2.1 Introduction

The asymmetric simple exclusion process (ASEP) is a minimal model of transport
in (quasi) one–dimensional systems. It consists of particles which occupy the
sites of a one dimensional lattice with only one particle allowed on each
lattice site. These particles hop under the effect of an external driving
force which breaks detailed balance and creates a stationary current. This
model was introduced in the late 60s in biology to model translation in
protein synthesis [25] and independently in probability [19] and afterwards it
has found a wide spectrum of applications, ranging from theoretical and
experimental studies of biophysical transport [17] to modeling traffic flow
[42, 69]. As soon one considers models which are more suited for
physical/biological systems, one will encounter variants of ASEP containing
localized or mobile defects and several species of particles, which have
different behaviors. As a result, typically these models are not exactly
solvable and even for some of the most basic questions, like the study of
large scale behavior of the system (which in the case of ASEP is known to be
described by the Burgers equation [18, 70, 71, 72]), approximations schemes
like mean–field are necessary.

In this paper we address the question of the large scale or hydrodynamic
behavior of an exactly solvable multispecies generalization of ASEP,
consisting of two kinds of particles, $\bullet$–particles and
$\circ$–particles, moving in opposite directions. One can think of them as
opposite charged particles moving under the influence of an external electric
field or as cars moving on two opposite lanes. Each site of a one–dimensional
lattice is either empty or occupied by one of the two kinds of particles. For
convenience, empty sites can be treated as a third species of particles, the
$\ast$–particles. In continuous time, a $\bullet$–particle jumps forward on
empty sites with rate $\beta$, while a white $\circ$–particle jumps backward
on empty sites with rate $\alpha$. On top of this, an adjacent pair
$\bullet\circ$ swaps to $\circ\bullet$ with rate $1$.

$\begin{split}\bullet\,\ast&\leavevmode\nobreak\
\rightarrow\leavevmode\nobreak\ \ast\,\bullet\quad\text{rate}\quad\beta\\\
\ast\,\circ&\leavevmode\nobreak\ \rightarrow\leavevmode\nobreak\
\circ\,\ast\quad\text{rate}\quad\alpha\\\ \bullet\,\circ&\leavevmode\nobreak\
\rightarrow\leavevmode\nobreak\ \circ\,\bullet\quad\text{rate}\quad
1\end{split}$ (2.1)

This model has appeared in the literature under different names. It has been
first considered in [32, 33], where the stationary measure on a finite
periodic lattice was written in a matrix product ansatz form [38, 16]. It is
also a particular case ($q=0$) of the so called AHR model [73, 74, 75], in
which the swap $\circ\bullet\rightarrow\bullet\circ$ is allowed with rate $q$.
Being a natural 2–species generalization of TASEP we shall call this model
$2$–TASEP. It turns out that the 2–TASEP is Yang–Baxter integrable, this was
first proven in [76] in a particular case with the constraint
$\alpha+\beta=1$, which happens to be the same condition for the system to
have a product invariant measure. For arbitrary values of $\alpha$ and
$\beta$, the Yang-Baxter integrability was proven in [46]. It belongs indeed
to a larger family of integrable multispecies exclusion processes introduced
in [77]. Bethe ansatz techniques can be used to solve exactly for the long
time limit behavior of the generating function of the currents [78, 46]. More
recently, in the case $\alpha+\beta=1$, the transition probabilities as well
as the joint current distribution for some specific initial distribution of a
finite number of $\bullet$ and $\circ$–particles have been obtained [79, 80],
and an asymptotic analysis of these results has allowed to prove that the
joint current distribution is given by a product of a Gaussian and a GUE
Tracy-Widom distribution in the long time limit, as predicted by non–linear
fluctuating hydrodynamics [81, 82, 83].

When $\alpha+\beta=1$ the stationary measure factorizes and the currents have
a simple expression as function of the densities. In [84, 85] the hydrodynamic
limit of the 2–TASEP for $\alpha=\beta=\frac{1}{2}$ has been studied and
proven to converge to the classical Leroux system of conservation laws [86,
87]. The Leroux system is a notable example of a _Temple class system_ i.e. a
2–components conservation law whose shock and rarefaction curves coincide
[68]. The theory of Temple class systems shares several common features with
the theory of single component conservation laws [88], in particular well-
posedness results for Temple systems are available for a much larger class of
initial data compared to general systems of conservation laws.

For arbitrary $\alpha$ and $\beta$ only numerical results based on mean field
approximation are available [89]. In the present paper we study the exact
hydrodynamic equations of the 2–TASEP and show that they are Temple class for
arbitrary $\alpha$ and $\beta$. This allows us to compute their rarefaction
and shock solution, as well as to solve completely the Riemann problem, which
consists in determining the density profile starting from a domain-wall
initial data.

The paper is organized as follows. In Section 2.2 we review and expand on
results about the 2–TASEP currents obtained in [46]. The core of the paper is
Section 2.3 where the conservation laws are studied. We derive the rarefaction
waves as well as the shock solutions and finally we solve the full Riemann
problem. In Section 2.4 we compare the prediction of the hydrodynamic
equations with Monte Carlo simulations. Conclusions and some outlooks for
further works are discussed in Section 2.5.

###  2.2 Currents

In this section we reproduce and expand the results of the analysis in [46] in
a convenient way, which makes manifest the symmetries of the model. In order
to compute the particle currents as functions of the local densities, we
consider our model on a periodic ring with a fixed number of particles of each
species. Call $M_{i}$ the number of particles of species $i$, they are related
to $N$, the length of the ring, by $N=M_{\bullet}+M_{\circ}+M_{\ast}$. Let the
system evolve starting at time $t=0$ from an arbitrary fixed configuration and
call $n_{i,j}(t)$ the number of swaps of consecutive ordered pairs of
particles of type $i,j$ up to time $t$. This number increases by $+1$ each
time two consecutive ordered particles of species $i,j$ exchange their
position $i,j\rightarrow j,i$. The average rate of swaps $i,j\rightarrow j,i$
in the steady state is just given by
$\lim_{t\rightarrow+\infty}\frac{1}{t}\mathbb{E}\left[n_{i,j}(t)\right]$,
irrespectively of the initial state. The particle currents in the steady state
are hence given by

$J_{i}=\lim_{t\rightarrow+\infty}\frac{1}{Nt}\mathbb{E}\left[\sum_{j\neq
i}n_{i,j}(t)-n_{j,i}(t)\right]$ (2.2)

In our case, it is convenient to introduce the following quantity

$\Phi(\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})=\lim_{t\rightarrow+\infty}\frac{1}{t}\mathbb{E}\left[\nu_{\bullet,\circ}\,n_{\bullet,\circ}(t)+\nu_{\bullet,\ast}\,n_{\bullet,\ast}(t)+\nu_{\ast,\circ}\,n_{\ast,\circ}(t)\right].$
(2.3)

The currents are obtained as specialization of
$\Phi(\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})$

$J_{\bullet}=\frac{\Phi(1,1,0)}{N},\quad
J_{\circ}=\frac{\Phi(-1,0,-1)}{N},\quad J_{\ast}=\frac{\Phi(0,-1,1)}{N}.$
(2.4)

In [46] the function
$\Phi(\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})$ was shown to
be given by the solution of the following equation

$\det
G(\Phi(\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ}),\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})=0.$
(2.5)

where the matrix
$G(\Phi,\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})$ is given by

$G(\Phi,\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})=\scalebox{0.75}{$\left(\begin{array}[]{ccc}\Phi&F_{\alpha}[M_{\circ},M_{\bullet},M_{\ast}]&F_{\beta}[M_{\bullet},M_{\circ},M_{\ast}]\\\
\nu_{\bullet,\circ}M_{\bullet}+\nu_{\ast,\circ}M_{\ast}&F_{\alpha}[M_{\circ}+1,M_{\bullet},M_{\ast}]&-F_{\beta}[M_{\bullet},M_{\circ}+1,M_{\ast}]\\\
\nu_{\bullet,\circ}M_{\circ}+\nu_{\bullet,\ast}M_{\ast}&-F_{\alpha}[M_{\circ},M_{\bullet}+1,M_{\ast}]&F_{\beta}[M_{\bullet}+1,M_{\circ},M_{\ast}]\end{array}\right)$}$
(2.6)

with

$F_{\gamma}[a,b,c]:=\oint_{0}\frac{dz}{2\pi
i}\frac{1}{z^{a}(z-1)^{b}(z-\gamma)^{c}}.$ (2.7)

When one of the particle species is _strictly_ absent (i.e. when one among
$M_{\bullet},M_{\circ},M_{\ast}$ vanishes) the model reduces to a single
species TASEP and it is not difficult to see that one of the currents
vanishes, while the others boil down to the usual TASEP current. On the other
hand in the following we shall assume that at least one particle per species
is present ($M_{i}\neq 0$) and we shall be mainly interested in the
thermodynamic limit of these quantities as $N\rightarrow\infty$, with
$\lim_{N\rightarrow\infty}\frac{M_{i}}{N}=\rho_{i}$. We shall see, as already
found in [33, 32], that the presence of even a single particle of a given
species (i.e. an infinitesimally vanishing but not strictly zero density) can
affect the macroscopic behavior of the system. With this in mind, we consider
the limit $a\longrightarrow\infty$, with $b/a$ and $c/a$ fixed, of the
function $F_{\gamma}[a,b,c]$, that behaves like111Here we are supposing
$a,b,c,\gamma>0$.

$F_{\gamma}[a,b,c]\sim\frac{1}{z_{\gamma}^{a}(z_{\gamma}-1)^{b}(z_{\gamma}-\gamma)^{c}},$

where $z_{\gamma}$ is the zero of the saddle point equation
$\frac{a}{z}+\frac{b}{z-1}+\frac{c}{z-\gamma}=0$, belonging to the interval
$[0,\min[1,\gamma]]$. Applying this expression in eq.(2.5) we get in the
thermodynamic limit

$\lim_{N\rightarrow\infty}\frac{\Phi(\nu_{\bullet,\circ},\nu_{\bullet,\ast},\nu_{\ast,\circ})}{N}=(\nu_{\bullet,\circ}\rho_{\bullet}+\nu_{\ast,\circ}\rho_{\ast})z_{\alpha}(1-z_{\beta})+(\nu_{\bullet,\circ}\rho_{\circ}+\nu_{\bullet,\ast}\rho_{\ast})z_{\beta}(1-z_{\alpha})$
(2.8)

where with $z_{\alpha}\in[0,\min(1,\alpha)]$ and
$z_{\beta}\in[0,\min(1,\beta)]$ are the 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$
(2.9)
$\displaystyle\frac{\rho_{\bullet}}{z_{\beta}}+\frac{\rho_{\circ}}{z_{\beta}-1}+\frac{1-\rho_{\circ}-\rho_{\bullet}}{z_{\beta}-\beta}=0.$
(2.10)

The result for the currents then reads

$\displaystyle
J_{\circ}=z_{\alpha}(z_{\beta}-1)+\rho_{\circ}(z_{\alpha}-z_{\beta})$ (2.11)
$\displaystyle
J_{\bullet}=z_{\beta}(1-z_{\alpha})+\rho_{\bullet}(z_{\alpha}-z_{\beta})$
(2.12) $\displaystyle J_{\ast}=\rho_{\ast}(z_{\alpha}-z_{\beta})$ (2.13)

Notice that eqs.(2.9,4.29) are invariant under exchange
$\rho_{\circ}\leftrightarrow\rho_{\bullet}$, $\alpha\leftrightarrow\beta$ and
$z_{\alpha}\leftrightarrow z_{\beta}$. This implies as expected, that under
exchange $\rho_{\circ}\leftrightarrow\rho_{\bullet}$ and
$\alpha\leftrightarrow\beta$ we have $J_{\circ}\leftrightarrow-J_{\bullet}$.
Let us finish this section by showing how some known results fit in the
analysis above.

  * $\Diamond$

$\beta=1$. In this case $\bullet$–particles do not distinguish
$\circ$–particles from $\ast$–particles and so they behave just as particles
in a single species TASEP. This is reflected in eq.(4.29), where
$\rho_{\circ}$ disappears and one finds $z_{\beta}=\rho_{\bullet}$, which
replaced in eq.(4.27) gives $J_{\bullet}=\rho_{\bullet}(1-\rho_{\bullet})$.
The case $\alpha=1$ is completely analogous: $z_{\alpha}=\rho_{\circ}$ and
$J_{\circ}=\rho_{\circ}(\rho_{\circ}-1)$.

  * $\Diamond$

$\alpha+\beta=1$. In this case it is known that the stationary measure takes a
factorized form [75]. At the level of the currents, we have indeed
$J_{\circ}=-\rho_{\circ}(\rho_{\bullet}+\alpha\rho_{\ast})$ and
$J_{\bullet}=\rho_{\bullet}(\rho_{\circ}+\beta\rho_{\ast})$.

####  2.2.1 The $z$ variables

In our analysis the variables ${\bf z}=(z_{\alpha},z_{\beta})$ will play a
prominent role, it is therefore important to work out their domain of
definition $\mathcal{D}_{z}(\alpha,\beta)$ corresponding to the physical
domain $\mathcal{D}$ in the variables
$\bm{\rho}=(\rho_{\circ},\rho_{\bullet})$, $\rho_{\circ},\rho_{\bullet}\geq
0,\rho_{\circ}+\rho_{\bullet}\leq 1$.

First of all we have already seen above that ${\bf z}$ has to satisfy
$z_{\alpha}\in[0,\min(1,\alpha)]$ and $z_{\beta}\in[0,\min(1,\beta)]$. At
fixed ${\bf z}$, the system of equations (2.9,4.29) is just the crossing of
two lines in the $\bm{\rho}$ plane: $\ell_{\alpha}$ coming from eq.(2.9) and
$\ell_{\beta}$ coming from eq.(4.29). In Fig. 2.1 we show by a simple
geometrical argument that these lines cross inside $\mathcal{D}$ whenever
$z_{\alpha}+z_{\beta}\leq 1$. So in conclusion the domain
$\mathcal{D}_{z}(\alpha,\beta)$ is given by $z_{\alpha}\in[0,\min(1,\alpha)]$,
$z_{\beta}\in[0,\min(1,\beta)]$ and $z_{\alpha}+z_{\beta}\leq 1$. The
geometrical reasoning explained in Fig. 2.1 allows also to conclude that at
fixed $z_{\alpha}$, $\rho_{\bullet}$ is an increasing function of $z_{\beta}$,
while at fixed $z_{\beta}$, $\rho_{\circ}$ is an increasing function of
$z_{\alpha}$.

$1$$1$$\rho_{\circ}$$\rho_{\bullet}$$\ell_{\beta}$$\frac{z_{\beta}}{\beta}$$(1-z_{\beta},z_{\beta})$$\ell_{\alpha}$$\frac{z_{\alpha}}{\alpha}$$(z_{\alpha},1-z_{\alpha})$
Figure 2.1: The shaded triangle is the physical region $\mathcal{D}$. Given
$z_{\alpha}\in[0,\min(1,\alpha)]$, the line $\ell_{\alpha}$ (red) intersects
the boundary of $\mathcal{D}$ at the bottom side and at the diagonal side at
the point $C_{\alpha}$ of coordinates
$\rho_{\circ}=z_{\alpha},\rho_{\bullet}=1-z_{\alpha}$. Similarly given
$z_{\beta}\in[0,\min(1,\beta)]$, the line $\ell_{\beta}$ (blue) intersects the
boundary of $\mathcal{D}$ at the left side and at the diagonal side at the
point $C_{\beta}$ of coordinates
$\rho_{\circ}=1-z_{\beta},\rho_{\bullet}=z_{\beta}$. For the two lines to
cross we need the point $C_{\alpha}$ to lie above the point $C_{\beta}$, which
is the case iff $z_{\alpha}+z_{\beta}\leq 1$.

$\rho_{\circ}$$\rho_{\bullet}$$\frac{1}{\beta}$$\frac{1}{\alpha}$$1$$1$$z_{\alpha}$$z_{\beta}$$1$$1$
(a) $\alpha,\beta>1$

$\rho_{\circ}$$\rho_{\bullet}$$(1-\beta,\beta)$$(\alpha,1-\alpha)$$1$$1$$z_{\alpha}$$z_{\beta}$$\alpha$$\beta$$1$$1$
(b) $\alpha,\beta<1$, $\alpha+\beta>1$

$\rho_{\circ}$$\rho_{\bullet}$$1$$1$$z_{\alpha}$$z_{\beta}$$\alpha$$\beta$$1$$1$
(c) $\alpha+\beta<1$

Figure 2.2: On the left the physical domain $\mathcal{D}$ in the densities
plane, on the right the corresponding domain $\mathcal{D}_{z}$ in the ${\bf
z}$ variables plane. On both sides we have drawn in red the lines
$\ell_{\alpha}$ at constant $z_{\alpha}$ and in blue the lines $\ell_{\beta}$
at constant $z_{\beta}$. The lines $\ell_{\alpha}$ (red) have slope
$\frac{\alpha(1-z_{\alpha})}{(\alpha-1)z_{\alpha}}$ and in particular they
have positive slope for $\alpha>1$ as in figure (a) and negative slope for
$\alpha<1$ as in figures (b) and (c). Similarly, the lines $\ell_{\beta}$
(blue) have slope $\frac{(\beta-1)z_{\beta}}{\beta(1-z_{\beta})}$, they have
positive slope for $\beta>1$ as in figure (a) and negative slope for $\beta<1$
as in figures (b) and (c).

In Figs. 2.2 we have reported on the left the physical domain $\mathcal{D}$ in
the densities plane and on the right the corresponding domain
$\mathcal{D}_{z}$ in the ${\bf z}$ variables plane for the cases
$\alpha,\beta>1$ (2.2(a)) , $\alpha,\beta<1$ and $\alpha+\beta>1$ (2.2(b)) and
$\alpha+\beta<1$ (2.2(c)). Notice that in Fig. 2.2(a) the thick red segment on
the $\rho_{\circ}$ axis is mapped to the point $(z_{\alpha}=1,z_{\beta}=0)$
and the thick blue segment on the $\rho_{\bullet}$ axis is mapped to the point
$(z_{\alpha}=0,z_{\beta}=1)$. In Fig. 2.2(c): the overlap of the red and blue
segments on the boundary $\rho_{\ast}=0$ is mapped to the point
$(z_{\alpha}=\alpha,z_{\beta}=\beta)$. This indicates that the mapping
$\bm{\rho}\rightarrow{\bf z}$ can be singular on the boundary of the physical
domain $\mathcal{D}$, where at least one of the densities vanishes. Let’s
analyze the different possibilities and work out the portion of the boundary
where the mapping is singular.

  * $\rho_{\circ}\rightarrow 0$ In this case from eq.(2.9) we deduce that $z_{\alpha}=0$ while the two solutions of eq.(4.29) are $z_{\beta}=1,z_{\beta}=\rho_{\bullet}\beta$ and we have to retain the smallest one. If $\beta\leq 1$ then on the $\rho_{\circ}=0$ axis the map is $1$-to-$1$ and there are no singularities, on the other hand, if $\beta>1$ then all the points $\rho_{\bullet}\geq\beta^{-1}$ are mapped to the same point $(z_{\alpha}=0,z_{\beta}=1)$.

  * $\rho_{\bullet}\rightarrow 0$ This case is treated similarly to the previous one: we have $z_{\beta}=0$ and the two solutions of eq.(2.9) are $z_{\alpha}=1,z_{\alpha}=\rho_{\circ}\alpha$. If $\alpha\leq 1$ then on the $\rho_{\bullet}=0$ axis the map is $1$-to-$1$ and there are no singularities, on the other hand, if $\alpha>1$ then all the points $\rho_{\circ}\geq\alpha^{-1}$ are mapped to the same point $(z_{\alpha}=1,z_{\beta}=0)$.

  * $\rho_{\ast}\rightarrow 0$ Eqs.(2.9, 4.29) have solutions, $z_{\alpha}=\alpha,z_{\alpha}=\rho_{\circ}$, $z_{\beta}=\beta,z_{\beta}=1-\rho_{\circ}$. Whenever $\alpha+\beta\leq 1$, all the points on the line $\rho_{\circ}+\rho_{\bullet}=1$ such that $\alpha\leq\rho_{\circ}\leq 1-\beta$ are mapped to the single point $z_{\alpha}=\alpha,z_{\beta}=\beta$.

####  2.2.2 Behavior at the boundary of the physical domain

The singularities of the mapping $\bm{\rho}\rightarrow{\bf z}$ reflect some
important features of the model. Let’s consider the currents of non zero
density particles at the boundary of $\mathcal{D}$

$\rho_{\circ}\rightarrow 0$ We have for the current

$J_{\bullet}(\rho_{\bullet})=\left\\{\begin{array}[]{ll}\beta\rho_{\bullet}(1-\rho_{\bullet})&0\leq\rho_{\bullet}\leq\beta^{-1}\\\
(1-\rho_{\bullet})&\beta^{-1}\leq\rho_{\bullet}\leq 1.\end{array}\right.$
(2.14)

$\rho_{\bullet}\rightarrow 0$ We have for the current

$J_{\circ}(\rho_{\circ})=\left\\{\begin{array}[]{ll}-\alpha\rho_{\circ}(1-\rho_{\circ})&0\leq\rho_{\circ}\leq\alpha^{-1}\\\
-(1-\rho_{\circ})&\alpha^{-1}\leq\rho_{\circ}\leq 1.\end{array}\right.$ (2.15)

$\rho_{\ast}\rightarrow 0$ We have for the current

$J_{\bullet}(\rho_{\bullet})=\left\\{\begin{array}[]{ll}\rho_{\bullet}(1-\rho_{\bullet})&0\leq\rho_{\bullet}\leq\beta,1-\alpha\leq\rho_{\bullet}\leq
1\\\ \beta(1-\alpha)+(\alpha-\beta)\rho_{\bullet}&\beta\leq\rho_{\bullet}\leq
1-\alpha.\end{array}\right.$ (2.16)

$\rho_{\bullet}$$\frac{J_{\bullet}(\rho_{\bullet})}{\beta}$$1$$\beta^{-1}$ (a)

$\rho_{\circ}$$-\frac{J_{\circ}(\rho_{\circ})}{\alpha}$$1$$\alpha^{-1}$ (b)

$\rho_{\bullet}$$J_{\bullet}(\rho_{\bullet})$$1$$\beta$$1-\alpha$ (c)

Figure 2.3: (a) Current $J_{\bullet}$ for $\beta>1$ and $\rho_{\circ}=0$. (b)
Current $J_{\circ}$ for $\alpha>1$ and $\rho_{\bullet}=0$. (c) Current
$J_{\bullet}$ for $\alpha+\beta<1$ and $\rho_{\ast}=0$. The dashed lines
correspond to the current for a strict absence of $\circ$–particles (a),
$\bullet$–particles (b), $\ast$–particles (c).

These results have to be compared to the situation in which we have _strict_
absence of a species of particles (and not just vanishing density). Consider
for example a system without $\circ$–particles. Such a system is effectively a
single species TASEP with jump rates equal to $\beta$ and hence with current
just $J_{\bullet}(\rho_{\bullet})=\beta\rho_{\bullet}(1-\rho_{\bullet})$.
Comparing this with eq.(2.14) we see that for $\beta>1$ this behavior holds
only for $0\leq\rho_{\bullet}\leq\beta^{-1}$, while for
$\rho_{\bullet}>\beta^{-1}$ the presence of even a single $\circ$–particle
affects the macroscopic behavior of the system, giving rise to a modified
current.

At the boundary of $\mathcal{D}$, the average speed of the zero density
particles displays also an interesting behavior.

  * $\rho_{\circ}\rightarrow 0$ Speed of $\circ$–particles

$v_{\circ}(\rho_{\bullet})=\left\\{\begin{array}[]{ll}-\frac{\alpha+\beta(1-\alpha)\rho_{\bullet}(1-\rho_{\bullet})}{1+(\alpha-1)\rho_{\bullet}}-\beta\rho_{\bullet}&0\leq\rho_{\bullet}\leq\beta^{-1}\\\
-1&\beta^{-1}\leq\rho_{\bullet}\leq 1,\end{array}\right.$ (2.17)

  * $\rho_{\bullet}\rightarrow 0$ Speed of $\bullet$–particles

$v_{\bullet}(\rho_{\circ})=\left\\{\begin{array}[]{ll}\frac{\beta+\alpha(1-\beta)\rho_{\circ}(1-\rho_{\circ})}{1+(\beta-1)\rho_{\circ}}-\alpha\rho_{\circ}&0\leq\rho_{\circ}\leq\alpha^{-1}\\\
1&\alpha^{-1}\leq\rho_{\circ}\leq 1,\end{array}\right.$ (2.18)

  * $\rho_{\ast}\rightarrow 0$ Speed of $\ast$–particles

$v_{\ast}(\rho_{\bullet})=z_{\alpha}-z_{\beta}$ (2.19)

with

$z_{\alpha}=\left\\{\begin{array}[]{cc}\alpha&\rho_{\bullet}\leq 1-\alpha\\\
1-\rho_{\bullet}&\rho_{\bullet}\geq 1-\alpha\end{array}\right.,\qquad
z_{\beta}=\left\\{\begin{array}[]{cc}\beta&\rho_{\bullet}\geq\beta\\\
\rho_{\bullet}&\rho_{\bullet}\leq\beta\end{array}\right.$ (2.20)

The results in eqs.(2.17–2.19) have been obtained in the literature by
considering systems with a single particle of either species: the cases
$\rho_{\ast}\rightarrow 0$, eqs.(2.16,2.19) first appeared in [33, 78], the
cases $\rho_{\circ}\rightarrow 0$ or $\rho_{\bullet}\rightarrow 0$,
eqs.(2.14,2.17) and eqs.(2.15,2.18) first appeared in [90].

###  2.3 Conservations laws

Under Euler–scaling (where site position and time scale as
$\epsilon^{-1}n,\epsilon^{-1}t$ for $\epsilon\rightarrow 0$) the density
profiles are expected to evolve deterministically as solutions of a system of
conservation laws. Consider initial data $\rho_{i}^{(0)}(x)$ and to such data
associate a family of initial conditions of the 2-TASEP of product Bernoulli
form, with local probability at site $n$ given by

$\mathbb{E}[\chi_{i}^{\epsilon}(n,t=0)]=\rho_{i}^{(0)}(\epsilon n),$

where $\chi_{i}^{\epsilon}(n,t)$ is the $i$–th species indicator function at
time $t$ and site $n$. We expect that the random variable
$\chi_{i}^{\epsilon}(\lfloor\epsilon^{-1}x\rfloor,\epsilon^{-1}t,)$ converges
for $\epsilon\rightarrow 0$ to a deterministic density profile. More precisely
we expect that

$\lim_{\epsilon\longrightarrow 0}\sum_{n:a\leq\epsilon n\leq
b}\epsilon\,\chi_{i}^{\epsilon}(n,\epsilon^{-1}t)=\int_{a}^{b}\rho_{i}(x,t)dx,\qquad
a.s.$ (2.21)

where $\bm{\rho}=(\rho_{\circ},\rho_{\bullet})$ is the solutions of a system
of conservation laws

$\partial_{t}\bm{\rho}+\partial_{x}{\bf J}=0.$ (2.22)

with initial condition
$\bm{\rho}(t=0)=\bm{\rho}^{(0)}=(\rho^{(0)}_{\circ},\rho^{(0)}_{\bullet})$. By
making the usual hypothesis of _local stationarity_ we identify the local
currents with the stationary currents at density ${\bm{\rho}}$, given by
eqs.(2.11,4.27). A more precise statement and proof of this result for the
case $\alpha=\beta=\frac{1}{2}$ can be found in [84]. While the approach
developed in [84] can be extended to the full $\alpha+\beta=1$ line (for which
the stationary measure is product), it is not expected to work for arbitrary
values of $\alpha$ and $\beta$ [91]. In the present paper we take
eqs.(2.21,2.22) as a working hypothesis. Eqs.(2.22) form a system of coupled
conservation laws, whose non-linearity is known to be at the origin of
characteristic phenomena such as shocks formation in finite time, and
rarefaction waves, i.e. self-similar solutions, which present regions
expanding in time at constant speed where the densities interpolate between
two boundary values. In the following we shall analyze in detail eqs.(2.22).
We shall show that the variables ${\bf z}$ are Riemann variables for this
system. On general grounds we know that the rarefaction fans can be expressed
in an implicit form involving the Riemann variables. What is more surprising
is that for our system also the shock solutions are explicitly written in
terms of the Riemann variables: they correspond to a discontinuity of only one
of the two Riemann variables (the other one being continuous). Putting
together fans and shock we can explicitly solve the Riemann problem. In the
Section 2.4 these theoretical results are compared to Monte Carlo simulations
of the $2$–TASEP.

####  2.3.1 The cases $\alpha=\beta=1$ and $\alpha+\beta=1$

Before discussing the system of equation (2.22) in full generality let us
start with some comments on two particular cases: $\alpha=\beta=1$ and
$\alpha+\beta=1$

_$\bullet\quad\alpha=\beta=1$_.

When $\beta=1$, the $\bullet$–particles don’t distinguish $\circ$–particles
from $\ast$–particles. This means that the $\bullet$–particles evolve as in a
single species TASEP. At the level of currents, for $\beta=1$ we have indeed
$J_{\bullet}=\rho_{\bullet}(1-\rho_{\bullet})$. In this case the conservation
law for $\rho_{\bullet}$ completely decouples from that of $\rho_{\circ}$ and
takes the usual form of the non-viscous Burgers equation

$\partial_{t}\rho_{\bullet}+(1-2\rho_{\bullet})\partial_{x}\rho_{\bullet}=0.$

Analogously, for $\alpha=1$, the conservation law for $\rho_{\circ}$
completely decouples from that of $\rho_{\bullet}$ and takes the form

$\partial_{t}\rho_{\circ}-(1-2\rho_{\circ})\partial_{x}\rho_{\circ}=0.$

So for $\alpha=\beta=1$, system (2.22) just decouples completely into two
Burgers equations.

_$\bullet\quad\alpha+\beta=1$_.

In this case, thanks to the factorization of the stationary measure, the
currents can be explicitly written as functions of the densities

$J_{\circ}=-\rho_{\circ}(\rho_{\bullet}+\beta(1-\rho_{\circ}-\rho_{\bullet})),\qquad
J_{\bullet}=\rho_{\bullet}(\rho_{\circ}+\alpha(1-\rho_{\circ}-\rho_{\bullet})).$

One can consider conserved quantities $\rho$ and $v$, defined by

$\left(\begin{array}[]{c}\rho\\\
v\end{array}\right)=\left(\begin{array}[]{cc}-\frac{\alpha(\alpha+2\beta)}{3}&-\frac{\beta(\alpha+2\beta)}{3}\\\
\alpha&-\beta\end{array}\right)\left(\begin{array}[]{c}\rho_{\circ}\\\
\rho_{\bullet}\end{array}\right)+\left(\begin{array}[]{c}\frac{(2\alpha+\beta)(2\beta+\alpha)}{9}\\\
\frac{\beta-\alpha}{3}\end{array}\right)$ (2.23)

The associated currents are (up to irrelevant additive constants)