Some Exact Results for the Exclusion Process

# Some Exact Results for the Exclusion Process

## Abstract

The asymmetric simple exclusion process (ASEP) is a paradigm for non-equilibrium physics that appears as a building block to model various low-dimensional transport phenomena, ranging from intracellular traffic to quantum dots. We review some recent results obtained for the system on a periodic ring by using the Bethe Ansatz. We show that this method allows to derive analytically many properties of the dynamics of the model such as the spectral gap and the generating function of the current. We also discuss the solution of a generalized exclusion process with -species of particles and explain how a geometric construction inspired from queuing theory sheds light on the Matrix Product Representation technique that has been very fruitful to derive exact results for the ASEP.

Non-equilibrium statistical physics; ASEP; Bethe Ansatz; large deviations, Matrix Ansatz.
###### pacs:
05-40.-a;05-60.-a

## I Introduction

Equilibrium Statistical Mechanics is a well-established field that embodies Classical Thermodynamics (see e.g. F. Reif, 1965 (82)). In an extremely terse formulation, this theory could be summarized by saying that the relevant features of the microscopic physics of any system are encoded in a probability measure for which, at thermal equilibrium, a mathematical formula is known. Indeed, a system in contact with a heat bath at temperature occupies all accessible microscopic configurations with a probability , given by the Boltzmann-Gibbs canonical law:

 Peq(C)=e−E(C)/kTZ, (1)

being the energy of the configuration . The Partition Function (or ‘sum over states’) is related to the thermodynamic Free Energy F via the relation

 F=−kTLogZ.

This relation is nothing but an avatar of the celebrated Boltzmann law, . Macroscopic observables are obtained by taking the expectation values of the corresponding microscopic observables with respect to the canonical measure (1). These postulates provide us with a well-defined prescription to analyze systems at equilibrium. In particular, equilibrium statistical mechanics predicts macroscopic fluctuations (typically Gaussian) that are out of reach of Classical Thermodynamics: the paradigm of such fluctuations is the Brownian Motion.

For systems far from equilibrium, a fundamental theory that would generalize the formalism of equilibrium Statistical Mechanics to time-dependent processes is not yet available. A schematic non-equilibrium process can be represented as a rod in contact with two reservoirs at different temperatures, or at different electrical (chemical) potentials (see Figure 1). In the stationary regime, a constant current flows through the system. Many fundamental questions remain to be answered in order to understand the physics of such a simple system: What are the relevant parameters that would fully characterize the macroscopic state of a non-equilibrium system? Can the stationary state be defined as the optimum of some (unknown) macroscopic functions of some (unknown) parameters? Does a general equation of state exist? Can one classify non-equilibrium processes into ‘Universality classes’? Can one postulate a general form for some non-equilibrium microscopic measures? What do the fluctuations in the vicinity of a stationary state look like?

A large amount of research has been devoted to these questions and to some related ones. Although the theory is far from being complete, substantial progress has been made, particularly during the last twenty years. There are different ways to try to tackle the profound questions raised above. One line of research consists in exploring structural properties of non-equilibrium systems: this endeavour has led to celebrated results such as Fluctuation Theorems (Gallavotti and Cohen, 1995 (42); see also (63)), Non-equilibrium Work Relations (Jarzynski, 1997 (52)) or Macroscopic Fluctuation Theory (G. Jona-Lasinio and co-workers, see e.g. L. Bertini et al. (16)).

Another strategy to gain insight into non-equilibrium physics is to extract as much information as possible from analytical studies and from exact solutions of some special models. The Ising model, that has played a fundamental role in the theory of phase transitions, critical exponents and renormalization group, is a landmark of this style of research. In the field of non-equilibrium statistical mechanics, the Asymmetric Simple Exclusion Process (ASEP) is reaching the status of such a paradigm. The ASEP consists of particles on a lattice, that hop from a site to its immediate neighbours, and satisfy the exclusion condition: a given location must be occupied by at most one particle. Therefore, a jump is allowed only if the target site is empty. Physically, the exclusion constraint mimics short-range interactions amongst particles. Besides, in order to drive this lattice gas out of equilibrium, non-vanishing currents must be established in the system. This can be achieved by various means: by starting from non-uniform initial conditions, by coupling the system to external reservoirs that drive currents (57) through the system (transport of particles, energy, heat) or by introducing some intrinsic bias in the dynamics that favors motion in a privileged direction. Then, each particle is an asymmetric random walker that drifts steadily along the direction of an external driving force. Due to its simplicity, this model has appeared in different contexts. It was first proposed as a prototype to describe the dynamics of ribosomes along RNA (67); (68). In the mathematical literature, Brownian processes with hard-core interactions were defined by Spitzer (88) who coined the name exclusion process (see also (50); (64); (65)). The ASEP also describes transport in low-dimensional systems with strong geometrical constraints (33) such as macromolecules transiting through capillary vessels (14), anisotropic conductors, or quantum dots where electrons hop to vacant locations and repel each other via Coulomb interaction (103). Very popular modern applications of the exclusion process include molecular motors that transport proteins along filaments inside the cells and, of course, ASEP and its variants are ubiquitous in discrete models of traffic flow (34); (83). More generally, the ASEP belongs to the class of driven diffusive systems defined by Katz, Lebowitz and Spohn in 1984 (55). For a general discussion, we refer to e.g., the book of H. Spohn (91) the review of B. Schmittmann and R. K. P. Zia (81) and that of G. M. Schütz (86). An early review of the properties of the ASEP can be found in Derrida, 1998 (24). We emphasize that the ASEP is defined through dynamical rules: there is no energy associated with a microscopic configuration and thus there is no possibility of writing the stationary measure in the canonical form (1). More generally, the kinetic point of view seems to be a promising and fruitful approach to non-equilibrium systems: for such a presentation of statistical mechanics, we highly recommend to the reader the recent book by P. L. Krapivsky et al. (56).

To summarize, the ASEP is a minimal model to study non-equilibrium behaviour (see Figure 2). It is simple enough to allow analytical studies, however it contains the necessary ingredients for the emergence of a non-trivial phenomenology (106):

• ASYMMETRIC: The external driving breaks detailed-balance and creates a stationary current in the system (107). The model exhibits a non-equilibrium stationary state.

• EXCLUSION: The hard core-interaction implies that there is at most 1 particle per site. The ASEP is a genuine N-body problem.

• PROCESS: The dynamics is stochastic and Markovian: there is no underlying Hamiltonian.

The outline of this article is as follows. In section II, we study spectral properties of the ASEP on a ring by using the coordinate Bethe Ansatz. We review the general technique and focus on the TASEP case, which, in our opinion, is one of the simplest models that allows to understand how the Bethe Ansatz method works. In section III, we explain how the fluctuations of the total current in the ASEP can be calculated using a functional formulation of the Bethe Ansatz; we review the exact results obtained and the different scaling regimes that appear in the limit of systems of large sizes. In section IV, we study a generalization of the ASEP in which different classes of particles interact through hierarchical dynamical rules. We show how the steady state of such systems can be determined by using a Matrix Product Representation that involves tensor products of quadratic algebras.

## Ii Spectral Properties of the Markov Matrix

### ii.1 The model

We consider the exclusion process on a periodic one dimensional lattice with sites (sites and are identical) and particles. Because a lattice site cannot be occupied by more than one particle, the state of a site () can be characterized by the Boolean number according as the site is empty or occupied.

The system evolves in continuous time according to the following stochastic rule: a particle on a site at time jumps, in the interval between and , with probability to the neighbouring site if this site is empty (exclusion rule) and with probability to the site if this site is empty (see Figure 3). In the totally asymmetric exclusion process (TASEP) the jumps are totally biased in one direction (). On the other hand, the symmetric exclusion process (SEP) corresponds to the choice .

The number of particles is conserved by the dynamics. The total number of configurations for particles on a ring with sites is given by .

A configuration can be represented by the sequence We call the probability of configuration at time . As the exclusion process is a continuous-time Markov process, the time evolution of is determined by the master equation

 ddtPt(C)=∑C′M(C,C′)Pt(C′). (2)

The Markov matrix encodes the dynamics of the exclusion process: the element is the transition rate from configuration to and the diagonal term represents the exit rate from configuration .

A right eigenvector is associated with the eigenvalue of if

 Mψ=Eψ. (3)

The matrix is a real non-symmetric matrix and, therefore, its eigenvalues (and eigenvectors) are either real numbers or complex conjugate pairs. The spectrum of contains the eigenvalue and the associated right eigenvector is the stationary state. For the ASEP on a ring the steady state is uniform and all configurations have the same probability (24).

Because the dynamics is ergodic (i.e., is irreducible and aperiodic), the Perron-Frobenius theorem (see, for example, Gantmacher 1959 (43)) implies that 0 is a non-degenerate eigenvalue and that all other eigenvalues have a strictly negative real part; the relaxation time of the corresponding eigenmode is . The imaginary part of gives rise to an oscillatory behaviour.

In Figure 4, we display the example of a spectrum of a Markov matrix for the TASEP with distinguishable particles and a ring of sites. In this case, the dimension of the phase space is 1260. However, when the particles are indistinguishable, the dimension of the Markov matrix is reduced by a factor of and is given by 252. The spectrum of the indistinguishable case is included in the larger spectrum of the distinguishable problem.

### ii.2 The Bethe Ansatz

Another way to characterize a configuration is to specify the positions of the particles on the ring, with . In this representation, the eigenvalue equation (3) becomes

 Eψ(ξ1,…,ξN)= ∑i[ψ(ξ1,…,ξi−1, ξi−1, ξi+1,…,ξN)−ψ(ξ1,…,ξN)] + ∑jx[ψ(ξ1,…,ξj−1, ξj+1, ξj+1,…,ξN)−ψ(ξ1,…,ξN)], (4)

where the sums run over the indexes such that and over the indexes such that . In other words, the sums are restricted to jumps that respect the exclusion condition (24).

Since the works of Dhar 1987 (31), and Gwa and Spohn 1992 (48), it is known that the Bethe Ansatz can be applied to the ASEP. The idea of the Bethe Ansatz (H. Bethe 1931 (11)) consists in writing the eigenvectors of the Markov matrix as linear combinations of plane waves:

 ψ(ξ1,…,ξN)=∑σ∈ΣNAσzξ1σ(1)zξ2σ(2)…zξNσ(N), (5)

where is the group of the permutations of indexes. The coefficients and the wave-numbers are complex numbers to be determined. We observe that in the case where all particles are far from each other (i.e. no particles are located on adjacent sites so that for all ), each monomial that appears on the right hand side of the expression of the Bethe wave function is a solution of the eigen-equation (4), with an eigenvalue given by

 E(z1,z2…zN)=N∑i=11zi+xN∑i=1zi−N(1+x). (6)

However, in order for the trial wave function to be a genuine eigenfunction, equation (4) must be fulfilled for all configurations: therefore, one has to insure that when two or more particles are adjacent equation (4) is still satisfied. The case where exactly two particles are adjacent and all other particles are well separated from one another is enough to fix the value of all the ’s (up to an overall multiplicative constant). In other words, the form of the Bethe wave-function (5) is fully determined by the ‘two-body collisions’. However, more particles can form larger clusters: this corresponds to -body collisions with : such collisions impose further constraints that the Bethe wave-function has no reasons, a priori, to satisfy. Fortunately, in the present problem, the constraints imposed by the -body collisions can be written as linear combinations of the two-body constraints and are therefore automatically satisfied by the Bethe wave-function. In the Bethe Ansatz jargon, this fact is formulated by saying that the -body collisions factorize into -body collisions etc… This remarkable property, that can be verified explicitly in the case of equation (4), lies at the heart of the integrability of the exclusion process, i.e., it implies that the ASEP can be solved by Bethe Ansatz. In fact, the ASEP can be mapped exactly into well-studied systems such as quantum spin chains (2); (1), vertex models (9); (54); (101) or solid-on-solid models (80) which are well-known to be integrable.

Because we are studying the system on a homogeneous ring, the function must also satisfy the following periodic boundary conditions

 ψ(ξ1,ξ2,…,ξn)=ψ(ξ2,…,ξn,ξ1+L). (7)

These periodic conditions quantify the eigenvalues by imposing a set of equations satisfied by the ’s, the Bethe Equations:

 zLi=(−1)N−1N∏j=1xzizj−(1+x)zi+1xzizj−(1+x)zj+1. (8)

This set of non-linear algebraic equations must be satisfied by the fugacities. Therefore, in order to obtain the spectrum of the ASEP Markov matrix one has first to find all -tuplets that solve the Bethe equations (8). Then, given a solution, one can calculate the corresponding eigenvalue (6) and eigenfunction (5). Of course, in practice this program can be carried out only in some special situations (for example in a limited portion of the spectrum) and usually one has to take the ‘thermodynamic’ limit with a fixed value of the density . Besides, the completeness issue (i.e. whether the Bethe Ansatz does provide the full spectrum) is a difficult problem in algebraic geometry (10); (62).

### ii.3 Bethe Equations for the TASEP

The totally asymmetric exclusion process (TASEP), that corresponds to , provides a particularly instructive illustration of the Bethe Ansatz. The Bethe equations take a simpler form on which analytical calculations can be carried out even for finite values of and . The TASEP is one of the simplest non-trivial models that allows to understand how the Bethe Ansatz technique works. In this subsection, we present a complete and self-contained derivation of the Bethe Ansatz for the TASEP.

First, we show that the Bethe wave function of the TASEP can be written as a determinant (45). The form of this determinant can be guessed heuristically by working out explicitly examples for small systems. Here, we reverse the logic and define as

 ψ(ξ1,…,ξN)=det(R), (9)

where is a matrix with elements

 R(i,j)=zξji(1−zi)j  for 1≤i,j≤N, (10)

being complex numbers. By expanding the determinant, one recovers the generic form (5) for the Bethe wave function . We now show that , thus defined, is a solution of the eigenvalue equation (4) with . First, we need to prove that satisfies two identities which are valid for any values of and of . The first identity is

 Eψ(ξ1,…,ξN)=N∑k=1[ψ(ξ1,…,ξk−1,…,ξN)−ψ(ξ1,…,ξk,…,ξN)], (11)

where is given by We emphasize that this equality is true for any and without imposing the ordering . Besides, runs from 1 to without any restriction, i.e. we do not impose that . Equation (11) is proved by writing

 ψ(ξ1,…,ξk−1,…,ξN)−ψ(ξ1,…,ξk,…,ξN)= det(R(i,1),…,(1zi−1)R(i,k),…,R(i,N) ). (12)

This determinant is similar to except for the -th column. Expanding this determinant over all permutations of and performing the sum over leads to the desired equation. The second identity takes care of the two particles collision case: it is valid for any and any with , and is given by

 ψ(ξ1,…,ξk,ξk,…,ξn)−ψ(ξ1,…,ξk,ξk+1,…,ξn)=0. (13)

This identity is proved as follows: we rewrite the left hand side of equation (13) as the determinant where is a matrix that is identical to except for its -th column that is given by

 ~R(i,k)=zξki−zξk+1i(1−zi)k=zξki(1−zi)k−1=R(i,k−1)=~R(i,k−1). (14)

We remark that the -th and the -th columns of are equal and, therefore, , proving equation (13).

To conclude, we note that the eigenvalue equation (4) is similar to equation (11) except that in (4) the sum is restricted to the allowed jumps of particles, i.e., to the values of such that . However, in equation (11), the terms with vanish thanks to equation (13). Hence, equation (11) is in fact exactly the same as the eigenvalue equation when the eigenvector has the form assumed in equations (9, 10).

Finally, because of the periodic boundary conditions (7), the ’s are quantified by the Bethe equations, that we now rederive. Denoting by and the generic line and column of the matrix , we can write

 ψ(ξ2,…,ξN,ξ1+L)=det⎛⎝zξ2i1−zi,…,zξj+1i(1−zi)j,…,zξ1+Li(1−zi)N⎞⎠. (15)

By cyclic permutation of the columns, we obtain

 ψ(ξ2,…,ξN,ξ1+L) (16) = (−1)N−1det⎛⎝zξ1+Li(1−zi)N,zξ2i1−zi,…,zξji(1−zi)j−1,…⎞⎠ = (−1)N−1det(zLi(1−zi)N−1R(i,1),…,(1−zi)R(i,j),…) = (−1)N−1N∏k=1(1−zk) det(zLi(1−zi)NR(i,1),…,R(i,j),…).

The last expression will be equal to if are solutions of the TASEP Bethe equations:

 (zi−1)Nz−Li=−N∏k=1(1−zk)   for   i=1,…,N. (17)

These equations can of course be obtained by substituting in the general Bethe equations (8) for ASEP. We note that determinantal representations of the eigenvectors and of the exact probability distribution at finite time play an important role in the study of the TASEP since the seminal work of G. M. Schütz (85), that has been generalized further by V. B. Priezzhev (72) and N. M. Bogoliubov (13) (see e.g. the review of V. B. Priezzhev 2005 (73)).

### ii.4 Analysis of the TASEP Bethe Equations

The Bethe equations for the TASEP exhibit a remarkable ‘decoupling’ property that allows to perform analytical calculations even for finite values of and . Using the auxiliary variables , the Bethe equations (17) become

 (1−Zi)N (1+Zi)L−N=−2LN∏j=1Zj−1Zj+1   for   i=1,…,N. (18)

We note that the right-hand side of these equations is independent of the index : this property is true only for the TASEP where the Bethe equations decouple and can be reduced to a polynomial in one-variable plus a self-consistency equation as shall be explained below. The corresponding eigenvalue of the Markov matrix is given by

 2E=−N+N∑j=1Zj. (19)

The solutions of the Bethe equations (18) are the roots of the polynomial equation of degree

 (1−Z)N(1+Z)L−N=Y, (20)

where is determined self-consistently by the r.h.s. of equation (18). Given an arbitrary value of the complex number , the roots of the polynomial (20) display a simple geometrical layout. If we write , being a positive real number, we observe that equation (20) implies that

 |Z−1|ρ|Z+1|1−ρ=r (21)

where is the particle density in the system. The curves defined in the complex plane for are known as Cassini ovals (see Figure 5). As can be seen in Fig. 5, the topology of the Cassini ovals depends on the value of with a critical value:

 rc=2ρρ(1−ρ)1−ρ. (22)

For , the curve consists of two disjoint ovals with solutions on the oval surrounding and solutions on the oval surrounding . For , the curve is a deformed Bernoulli lemniscate with a double point at . For , the curve is a single loop with solutions. Note that the Cassini ovals are symmetrical only if ; in this case

This geometric layout leads to the following procedure for solving the Bethe Equations for the TASEP (for a review see e.g. (46)):

• SOLVE, for any given value of , the equation The roots are located on Cassini Ovals

• CHOOSE roots amongst the available roots, with a choice set

• SOLVE the self-consistent equation where

 Ac(Y)=−2LN∏j=1zc(j)−1zc(j)+1.
• DEDUCE from the value of , the ’s and the energy corresponding to the choice set  :

 2Ec(Y)=−N+N∑j=1zc(j).

Eigenvalues of the Markov matrix are thus related to the choice sets of roots on the Cassini ovals; there is a strong evidence (46) that choice sets are in one-to-one correspondence with the eigenvalues of the Markov matrix. A simple argument is that the total number of choice sets is given by which is precisely the size of the Markov matrix.

The choice function that selects the fugacities with the largest real parts (see Figure 5) provides the ground state of the Markov matrix. The spectral gap, given by the first excited eigenvalue, corresponds to the choice for and (48). For this choice set, the calculation can be carried out explicitly and one can show that, in the large limit, the first excited state is given by a solution of a transcendental equation. For a density , one obtains

 E1= −2√ρ(1−ρ)6.509189337…L3/2±2iπ(2ρ−1)L. (RELAXATION)(OSCILLATIONS)

The first excited state consists of a pair of conjugate complex numbers when is different from 1/2. The real part of describes the relaxation towards the stationary state: we find that the largest relaxation time scales as with the dynamical exponent (100); (31); (48); (60); (44). This value agrees with the dynamical exponent of the one-dimensional Kardar-Parisi-Zhang equation that belongs to the same universality class as ASEP (see the review of Halpin-Healy and Zhang 1995 (49)). The imaginary part of represents the relaxation oscillations and scales as ; these oscillations correspond to a kinematic wave that propagates with the group velocity : such traveling waves can be probed by dynamical correlations (66); (8).

The same procedure also allows to classify the higher excitations in the spectrum (22). For the partially asymmetric case (), the Bethe equations do not decouple and analytical results are much harder to obtain. A systematic procedure for calculating the finite size corrections of the upper region of the ASEP spectrum was developed by Doochul Kim (60); (61).

We emphasize that the results presented here are valid for the ASEP on a finite periodic lattice. Problems with open boundary conditions (21); (22); (23) or on infinite lattices (92) require different approaches. In particular, the exclusion process on the infinite line has fascinating connections with random matrix theory and with the Robinson-Schensted-Knuth correspondence. More precisely, consider the totally asymmetric case and suppose that the initial configuration of the TASEP is such that all sites are occupied and all sites with are empty. Then, the probability that a particle initially at moves at least steps to the right in time equals the probability distribution of the largest eigenvalue in a unitary Laguerre random matrix ensemble. This crucial fact, first understood by Johannson in 2000 (53), allows to relate the statistical properties of the totally asymmetric exclusion process to the classical Tracy-Widom laws (94) for the largest eigenvalue in ensembles of random matrices. The study of such relations has become a subfield per se and has stimulated many works (see e.g., (89); (90); (40)). In particular, the determinantal representation of the transition probabilities of the TASEP allows to retrieve Johannson’s result in a very appealing manner (84). More recently, in a series of articles (95); (96); (97); (98); (99), C. A. Tracy and H. Widom have found some integral formulas for the probability distribution of an individual particle starting from step initial conditions in the asymmetric exclusion process with general parameter values (i.e. allowing forward and backward jumps). These expressions can be rewritten as a single integral involving a Fredholm determinant, that is amenable to asymptotic analysis. In particular, a limit theorem is proven for the total current distribution. These breakthroughs extend the results of Johansson on TASEP to ASEP: this is a crucial progress because the weakly asymmetric process leads to a well-defined continuous limit. In particular, a very important outgrow of these studies are the recent papers of H. Spohn, T. Sasamoto and S. Prolhac, in which the height fluctuations of the Kardar-Parisi-Zhang equation and -point correlation functions are exactly derived, solving a problem that remained open for 25 years (see the articles of H. Spohn and T. Sasamoto, and of P. L. Ferrari in the present special issue). For an overview of all these fascinating problems, we refer the reader to the recent review article by Kriecherbauer and Krug 2010 (58).

Finally, we note that the Bethe Ansatz can be used in models that are closely related to the ASEP, such as the zero-range process (71), the raise and peel model (4), vicious walkers (32) or the Bernoulli matching model of sequence alignment (74).

## Iii Fluctuations of the total current

In this section, we explain how the statistics of the total current (38); (102) in the ASEP on a ring can be determined by the Bethe Ansatz. In particular, we show that the fluctuations of the current can display a non-Gaussian behaviour in contrast with the equilibrium case.

### iii.1 Current statistics and Bethe Ansatz

We call the total distance covered by all the particles between time 0 and time and define the joint probability of being in the configuration at time and having . The evolution equation of is:

 ddtPt(C,Y)=∑C′(M0(C,C′)Pt(C′,Y)+M1(C,C′)Pt(C′,Y−1)+M−1(C,C′)Pt(C′,Y+1)). (23)

Using the generating function

 Ft(C)=∞∑Y=0eγYPt(C,Y), (24)

the evolution equation becomes

 ddtFt(C)=∑C′(M0(C,C′)+eγM1(C,C′)+e−γM−1(C,C′))Ft(C′)=∑C′M(γ)(C,C′)Ft(C′). (25)

This equation is similar to the original Markov equation (2) for the probability distribution but now the original Markov matrix is deformed by a jump-counting fugacity into (which is not a Markov matrix in general), given by

 M(γ)=M0+eγM1+e−γM−1. (26)

In the long time limit, , the behaviour of is dominated by the largest eigenvalue and one can write

 ⟨eγYt⟩≃eE(γ)t. (27)

Thus, in the long time limit, the function is the generating function of the cumulants of the total current . But is also the dominant eigenvalue of the matrix . Therefore, the current statistics has been traded into an eigenvalue problem. Fortunately, the deformed matrix can still be diagonalized by the Bethe Ansatz. In fact, a small modification of the calculations described in Section II leads to the following Bethe Ansatz equations

 zLi=(−1)N−1N∏j=1xe−γzizj−(1+x)zi+eγxe−γzizj−(1+x)zj+eγ. (28)

The eigenvalues of are given by

 E(γ;z1,z2…zN)=eγN∑i=11zi+xe−γN∑i=1zi−N(1+x). (29)

The cumulant generating function corresponds to the largest eigenvalue.

Remark: One can define , the large-deviation function of the current, as follows

 Prob(Ytt=j)∼e−tG(j). (30)

It can be shown using (27) that and are the Legendre transforms of each other. Large deviation functions play an increasingly important role in non-equilibrium statistical physics (see H. Touchette 2009 (93)), in particular in relation to the Fluctuation Theorem (63). Thus, determining exact expression for these large deviation functions for interacting particle processes, either analytically or numerically is a major challenge in the field (we refer the reader to the review of B. Derrida, 2007 (25)). Besides, higher cumulants of the current and large deviations are also of experimental interest in relation to counting statistics in quantum systems (41).

### iii.2 The TASEP Case

For the TASEP, the Bethe equations (28) decouple and can be studied by using the procedure outlined in Section II.4. The case was completely solved by B. Derrida and J. L. Lebowitz in 1998 (29). These authors calculated by Bethe Ansatz to all orders in . More precisely, they obtained the following representation of the function in terms of an auxiliary parameter :

 E(γ) = −N∞∑k=1(kL−1kN)BkkL−1, (31) γ = −∞∑k=1(kLkN)BkkL. (32)

These expressions allow to calculate the cumulants of , for example the mean-current and the diffusion constant :

 J=limt→∞⟨Yt⟩t = dE(γ)dγ∣∣γ=0=N(L−N)L−1, (33) D=limt→∞⟨Y2t⟩−⟨Yt⟩2t = d2E(γ)dγ2∣∣γ=0=N2(2L−3)!(N−1)!2(L−N)!2(L−1)!2(2N−1)!(2L−2N−1)!. (34)

When , with a fixed density and , the large deviation function , defined in (30) can be written in the following scaling form:

 G(j)=√ρ(1−ρ)πN3H(j−Lρ(1−ρ)ρ(1−ρ)) (35)

with

 H(y)≃−2√35√πy5/2 for y→+∞, (36) H(y)≃−4√π3|y|3/2 for y→−∞. (37)

This large deviation function is not a quadratic polynomial, even in the vicinity of the steady state. Moreover, the shape of this function is skew: it decays as the exponential of a power law with an exponent for and with an exponent for .

### iii.3 The general case: Functional Bethe Ansatz

In the general case , the Bethe Ansatz equations do not decouple and a procedure for solving them was lacking. For example, it did not even seem possible to extract from the Bethe equations (28) a formula for the mean stationary current (which can be obtained very easily by other means from the fact that the stationary measure is uniform). Finally, the solution was found by rewriting the Bethe Ansatz as a functional equation and restating it as a purely algebraic problem. We explain here the method we followed (75); (77) and describe some of the results obtained. First, we perform the following change of variables,

 yi=1−e−γzi1−xe−γzi. (38)

In terms of the variables the Bethe equations read

 eLγ(1−yi1−xyi)L=−N∏j=1yi−xyjxyi−yjfori=1…N. (39)

Here again the equations do not decouple as soon as . However, these equations are now built from first order monomials in the ’s and they are symmetrical in these variables. This observation suggests to introduce an auxiliary variable that plays the same role with respect to all the ’s and allows to define the auxiliary equation:

 eLγ(1−T1−xT)L=−N∏j=1T−xyjxT−yjfori=1…N. (40)

This equation, in which is the unknown, and the ’s are parameters, can be rewritten as a polynomial equation:

 P(T)=eLγ(1−T)LN∏i=1(xT−yi)+(1−xT)LN∏i=1(T−xyi)=0. (41)

Because the Bethe equations (39) imply that for , the polynomial , defined as

 Q(T)=N∏i=1(T−yi), (42)

must divide the polynomial . Now, if we examine closely the expression of , we observe that the factors that contain the ’s inside the products over can be written in terms of . Therefore, we conclude that DIVIDES Equivalently, there exists a polynomial such that

 Q(T)R(T)=eLγ(1−T)LQ(xT)+xN(1−xT)LQ(T/x). (43)

This functional equation is equivalent to the Bethe Ansatz equations (it is also known as Baxter’s TQ equation). It can be used to determine the polynomial of degree that vanishes at the Bethe roots. In the present case, equation (43) can be solved perturbatively w.r.t. to any desired order. Knowing perturbatively an expansion of is derived, leading to the cumulants of the current and to the large deviation function. For example, this method allows to calculate the following cumulants of the total current:

Mean Current :

Diffusion Constant :

 D=(1−x)2LL−1∑k>0k2CN+kLCNLCN−kLCNL(1+xk1−xk).

(We note that this formula was previously derived (30) using the Matrix Product Representation that we discuss in the next section). In the limit of a large system size, , with asymmetry parameter and with a fixed value of , the diffusion constant assumes a simple expression

 D∼4ϕLρ(1−ρ)∫∞0duu2tanhϕue−u2.

Third cumulant: the Skewness measures the non-Gaussian character of the fluctuations. An exact combinatorial expression of the third moment, valid for any values of , and , was calculated by S. Prolhac in (76). It is given by

 E36L2 = 1−xL−1∑i>0∑j>0CN+iLCN−iLCN+jLCN−jL(CNL)4(i2+j2)1+xi1−xi1+xj1−xj − 1−xL−1∑i>0∑j>0CN+iLCN+jLCN−i−jL(CNL)3i2+ij+j221+xi1−xi1+xj1−xj − 1−xL−1∑i>0∑j>0CN−iLCN−jLCN+i+jL(CNL)3i2+ij+j221+xi1−xi1+xj1−xj − 1−xL−1∑i>0CN+iLCN−iL(CNL)2i22(1+xi1−xi)2 + (1−x)N(L−N)4(L−1)(2L−1)C2N2L(CNL)2 − (1−x)N(L−N)6(L−1)(3L−1)C3N3L(CNL)3.

For , and keeping fixed, this formula becomes

 E3ϕ(ρ(1−ρ))3/2L5/2≃−4π3√3+ 12∫∞0dudv(u2+v2)e−u2−v2−(u2+uv+v2)e−u2−uv−v2tanhϕutanhϕv.

This shows that the fluctuations display a non-Gaussian behaviour. We remark that for the TASEP limit is recovered:

 E3≃(32−83√3)π(ρ(1−ρ))2L3.

A systematic expansion procedure that completely solves the problem to all orders and yields exact expressions for all the cumulants of the current, for an arbitrary value of the asymmetry parameter , was carried out by S. Prolhac in (78). Using the functional Bethe Ansatz, S. Prolhac derived a parametric representation of the cumulant generating function similar to the one given for the TASEP in equations (31) and (32), but where the binomial coefficients are replaced by combinatorial expressions that are related to some tree structures. A closed expansion of w.r.t. was derived and the coefficients that appear at each order were given a combinatorial interpretation. This expansion, valid for any finite values of , and , was then used to study the large system size limit with various scalings of the asymmetry. Various regimes were found and the corresponding expressions for the cumulants were fully worked out:

• For , the model falls into the Edward-Wilkinson universality class.

• The range , where is a finite number, defines the weakly asymmetric regime (to be discussed below).

• The intermediate regime, corresponding to , exhibits a specific scaling behaviour that, to our knowledge, can not be represented by a continuous stochastic equation.

• For the system is in the strongly asymmetric regime.

• Finally, corresponds to the KPZ universality class, which contains the TASEP. This limit was also studied by Lee and Kim (59).

We conclude this section by some remarks specific to the weakly asymmetric case, for which the asymmetry parameter scales as in the limit of large system sizes . In this case, we also need to rescale the fugacity parameter as and the following asymptotic formula for the cumulant generating function can be derived

 ~E(γ,ν)≡E(γL,1−νL) ≃ ρ(1−ρ)(γ2+γν)L−ρ(1−ρ)γ2ν2L2+1L2ϕ[ρ(1−ρ)(γ2+γν)], (44) with ϕ(z) = ∞∑k=1B2k−2k!(k−1)!zk, (45)

and where the ’s are Bernoulli Numbers. We observe that the leading order (in ) is quadratic in and describes Gaussian fluctuations. It is only in the subleading correction (in ) that the non-Gaussian character arises. This formula was also obtained for the symmetric exclusion case in (6). We observe that the series that defines the function has a finite radius of convergence and that has a singularity for . Thus, non-analyticities appear in as soon as

 ν≥νc=2π√ρ(1−ρ).

By Legendre transform, non-analyticities also occur in the large deviation function defined in (30). At half-filling, the singularity appears at as can be seen in Figure 6. For the leading behaviour of is quadratic (corresponding to Gaussian fluctuations) and is given by

 G(j)=(j−νρ(1−ρ))24Lρ(1−ρ). (46)

For , the series expansions (44) and  (45) break down and the large deviation function becomes non-quadratic even at leading order. This phase transition was predicted by T. Bodineau and B. Derrida using macroscopic fluctuation theory (17); (18). These authors studied the optimal density profile that corresponds to the observation of the total current over a large time . They found that this optimal profile is flat for but it becomes linearly unstable for . In fact, when the optimal profile is time-dependent. The critical value of the total current for which this phase-transition occurs is and therefore one must have for this transition to occur. One can observe in Figure 6 that for , the large deviation function becomes non-quadratic and develops a kink at a special value of the total current .

## Iv Multispecies exclusion processes and Matrix Ansatz

From the mathematical point of view, the ASEP is one of the simplest but non-trivial model for which the hydrodynamic limit can be rigorously proved. At large scales, the distribution of the particles of the ASEP emerges as a density field that evolves according to Burgers equation with a vanishingly small viscosity (91). The Burgers equation is the textbook prototype for shock formation: a smooth initial distribution can develop a singularity (a discontinuity) in finite time. A natural question that arises is whether this shock is an artifact of the hydrodynamic limit or whether, under some specific conditions, the original ASEP does display some singularity at the microscopic scale (5); (37); (19); (36); (51); (69). The question was answered positively: an abrupt transition does exist at the level of the particle system and its width is of the order of the lattice size. However, defining precisely the position of the shock at the microscopic level requires some thought, and this was achieved by introducing a second-class particle (see e.g., P. A. Ferrari, 1992 (37)). This second class particle behaves as a first-class particle with respect to holes and as a hole with respect to first-class particles. The dynamical rules thus take the following form: where all transitions occur with rate 1. This dynamics corresponds to coupling two TASEP models. In order to locate the shock it is enough to introduce a single second class particle in the system.

A straightforward generalization of the two species case of first and second class particles is the multispecies process where there is a hierarchy amongst the different species: first class particles have highest priority and can overtake all other classes of particles; second class particles can overtake all other classes except the first class etc… hence, during an infinitesimal time step , the following processes take place on each bond with probability :

 I 0→0 I forI≠0 I J→J I for1≤I

Note that particles can always overtake holes (= 0-th class particles). This model will be called the -TASEP and it can be obtained by coupling ordinary TASEP models (64); (65).

Suppose that there are particles of class on a ring of size . Then, the total number of configurations is given by

 Ω=L!P0!P1!P2!…PN!

The total number of particles is given . We warn the reader that in this section denotes the number of species and not the number of particles.

The object of this section is to provide an algebraic description of the stationary measure of this system.

### iv.1 Matrix Ansatz for Two Species

The idea of representing the stationary weights of a configuration as traces over a quadratic algebra goes back to the seminal article of Derrida, Evans, Hakim and Pasquier in 1993 (27). These authors were studying the 1-species TASEP with open boundary conditions but they realized that the same idea could also be applied to the Two Species ASEP on a ring (Derrida, Janowski, Lebowitz and Speer, 1993 (28)). This technique, known as the Matrix Product Representation (or Matrix Ansatz), has been exceptionally fruitful in the field of one-dimensional stochastic models and has led to a very large number of exact solutions. This method seems to be complementary to the Bethe Ansatz in many problems and the exact relations between the two techniques is not yet fully understood (3); (47). A solver’s guide to the Matrix Ansatz method has recently been written by R.A. Blythe and M. R. Evans (12) and contains many applications to various models.

Here, we explain how the Matrix Ansatz works for the two species TASEP on a ring, with dynamical rules given by (47). A Configuration of the model can be represented by a string made from the ’letters’ 1,2 and 0, e.g., . According to the Matrix Ansatz, the stationary weight of is given by

 P(C)=P(01220211)=1Z