Combinatorial mappings of exclusion processes

We review various combinatorial interpretations and mappings of stationary-state probabilities of the totally asymmetric, partially asymmetric and symmetric simple exclusion processes (TASEP, PASEP, SSEP respectively). In these steady states, the statistical weight of a configuration is determined from a matrix product, which can be written explicitly in terms of generalised ladder operators. This lends a natural association to the enumeration of random walks with certain properties. Specifically, there is a one-to-many mapping of steady-state configurations to a larger state space of discrete paths, which themselves map to an even larger state space of number permutations. It is often the case that the configuration weights in the extended space are of a relatively simple form (e.g., a Boltzmann-like distribution). Meanwhile, various physical properties of the nonequilibrium steady state - such as the entropy - can be interpreted in terms of how this larger state space has been partitioned. These mappings sometimes allow physical results to be derived very simply, and conversely the physical approach allows some new combinatorial problems to be solved. This work brings together results and observations scattered in the combinatorics and statistical physics literature, and also presents new results. The review is pitched at statistical physicists who, though not professional combinatorialists, are competent and enthusiastic amateurs.

Appendix F Proof of determinantal form of the partition function 54 1. Motivation: nonequilibrium stationary states and combinatorics A physical system is said to be in a nonequilibrium steady state (NESS) when the probability distribution over microstates is independent of time, but there are nevertheless nonzero currents of some quantity (such as energy, mass, or, more generally, probability). In the presence of these currents, the microstate probability distribution is not calculable by the conventional equilibrium statistical mechanics approach, where a partition function of Boltzmann weights would give access to all macroscopic observables. An open problem is to find an equivalent unified approach to solve a NESS. This review focuses on a family of models with steady states that do permit exact solution: the asymmetric simple exclusion process (ASEP). These models, which involve hard-core particles hopping on a one-dimensional lattice, have endured for several reasons. First, they are motivated by physical phenomena, such as biophysical or vehicular transport [1,2], which are of interest in their own right. Second, their analytical tractability has provided a number of fundamental insights into the behaviour of systems that are driven out of equilibrium, and how they contrast with equilibrium systems. For example, a major finding was that one-dimensional driven systems with shortrange interactions can undergo phase transitions while their equilibrium counterparts can not [3,4]. Finally, the structure of the distribution over microstates that arises in the nonequilibrium setting is of interest from a more mathematical perspective.
It is this latter perspective that is the focus of this short review, with the particular aim of bringing together mathematical results that are scattered across the literature to the attention of statistical physicists. Our starting point is the matrix product solution of the ASEP, which was obtained by statistical physicists in the early 1990s [5][6][7][8]. The first solution [5] was obtained for special values of model parameters and was based on recursion relations between the statistical weight of configurations on lattices of different lengths. These recursion relations provide a clue that a deeper mathematical structure underlies this NESS. A major development was the introduction of matrix product expressions for the stationary weights, wherein the recursion relations are replaced by algebraic rules that the matrices and vectors involved must obey [7]. This approach allowed for physical quantities such as particle currents, density profiles and correlation functions, to be calculated more easily. Moreover, it also allowed solvable generalisations of the ASEP to be identified. The scope and application of the matrix product method is reviewed in detail in [9].
Furthermore, the exact solution of the ASEP stationary state reveals, with increasing system size, number sequences that are ubiquitous within enumerative combinatorics. Specifically, in Section 2 below, we note the appearance of factorials (which count permutations), Catalan numbers (which count a variety of objects, including the number of legal combinations of nested brackets) and ballot numbers (which count subsets of nested bracket combinations). In other contexts, discussed in Section 3, one also finds Narayana numbers (which count different subsets of nested brackets to ballot numbers) and Eulerian numbers (which count subsets of permutations).
The main contribution of this review is to bring together and unify various results and observations in the literature which may shed light on the question: why should these combinatorial sequences arise in a nonequilibrium physics problem?. At the broadest level, these sequences emerge from one-to-many mappings of configurations of particles in the ASEP to a larger set of objects (often, but not exclusively, paths on a lattice). For this we also introduce the idea of dominated paths. This is a new way of illustrating the larger set of objects, which offers an intuitive interpretation of the statistical weight of ASEP configurations, based on the paths they map to. We show how this is equivalent to mappings such as Motzkin paths, that are already known in the literature.
This mapping fully describes the configuration space for when the hopping of the particles is entirely asymmetric. When the particles can hop either left or right, the set of objects vastly grows in size, to permutations of numbers. We outline a mapping from the combinatorial literature where the statistical weight of an ASEP configuration is found by enumerating permutations that follow certain rules. We then make an observation that these permutations can in fact be viewed as a mapping of the dominated paths discussed earlier. We show that this effectively interpolates between the totally asymmetric and partially asymmetric variants of the ASEP.
We also review some cases where dynamical rules can be defined on these larger state spaces that, when projected onto the ASEP configuration space, recover the appropriate ASEP dynamics. In particular, the dynamical transitions in the larger space may allow a bijection between the set of transition rates out of and into any configuration, implying a uniform distribution on the larger state space.
Finally we review generalisations to multispecies processes, where in certain cases the stationary state may be expressed as a matrix product. The combinatorial implications for the multispecies case are currently being explored and our aim is to point out emerging connections with queueing theory, integrable systems and associated algebraic combinatorics.
On a more practical level, we provide various examples throughout the review where physical results for the ASEP (such as correlation functions and density profiles) can be derived very efficiently using known results for these enumeration problems. Conversely Figure 1. The exclusion process that constitutes the bulk of this review. Particles attempt to enter the system at rate α, attempt to move to the right at rate 1, to the left at rate q, and exit from the rightmost site at rate β. Particles may only move into free sites.
the matrix product approach may help solve otherwise challenging combinatorial problems [10].
We begin the review by discussing the the totally asymmetric, single particle species exclusion process, which turns out to be the easiest model to interpret in terms of combinatorics. We will introduce the mappings of this process, upon which the symmetric and partially asymmetric variants later generalise. We also discuss measures of entropy of these systems, and multispecies extensions. To begin, however, in the next section we give a summary of the matrix product formalism that solves these exclusion processes. The reader with some familiarity with these concepts may proceed to Section 2.3 where we make some simple observations of the matrix product formalism that allude to this combinatorial structure, or to Section 2.4 where we set out the structure of this review in more detail.

The asymmetric simple exclusion process
The ASEP ( Figure 1) is a stochastic open system of hopping particles. On a lattice of N sites, particles can make unit steps left and right, at respective rates q ≥ 0 and 1. They enter onto site 1 from the left at rate α, and leave from site N at rate β. The system is referred to as the symmetric simple exclusion process (SSEP) if q = 1, totally asymmetric (TASEP) if q = 0, and partially asymmetric (PASEP) for all other q. The system is sometimes referred to as weakly asymmetric (WASEP) if q → 1 in a system-size-dependent way e.g. q = 1 − O(N −1/2 ). In all variants the particles can only move onto vacant sites. In the long time limit, this system approaches a NESS, whereby the current of particles running left-to-right through the system stabilises.
This system is exactly solvable by a mathematical formalism known as the matrix product approach, introduced in [7]. The weight W of any configuration C of particles and vacant sites (or holes) can be expressed as an ordered product of matrices, one per site, that represent the sequence of occupied and vacant sites along the lattice. As we now review, this approach permits the derivation of quantities such as the nonequilibrium partition function, average density profile and the steady-state current, as a function of q, α, β. We note that the ASEP and matrix product solution can be generalised to a five-parameter model with two additional parameters γ and δ that correspond to exit of particles at the left and entry of particles at the right boundary [7,11,12]. However in this review we restrict ourselves to the three parameter case.

Overview of matrix product formalism of the ASEP
Denote an ASEP configuration of size N by C = (τ 1 , τ 2 , . . . τ N ), where τ i = 0 if site i is empty, and 1 if it is occupied. The weight of C is an ordered product of matrices, reduced to a scalar by two boundary vectors W |, |V : Now, if we define matrices X i (1) = D, X i (0) = E to represent occupied and vacant sites respectively, then the matrix product form (1) for the weight of a configuration is a solution to the underlying steady-state master equation, given that these matrices and vectors have the following properties: Note that for the TASEP where q = 0, relation (2) simplifies to The first proof that the configurational weights (1) form a stationary solution of the master equation for the ASEP was given in [7]. The quadratic algebra (2)(3)(4) can be used to reduce any string of matrices {D, E} N into strings that can be directly evaluated. As an example, the configuration C = (1, 0, 1) (a 3-site ASEP) has weight where we have taken W |V = 1. Here we see that repeated application of (2) generates a sum of matrix products. Once each term is in a normal-ordered form with all Es to the left and Ds to the right, Eqs. (3) and (4) can be used to convert each term to a scalar. In the following we will associate each term with the weight of a configuration in an expanded space of configurations. If we think of the rates α, β and q as each being the exponential of an energy-like quantity, then configurations in the expanded space can be thought of following a Boltzmann distribution with an appropriately-defined energy function. Then, each weight in the nonequilibrium ensemble is given by the sum over weights in an equilibrium ensemble.
The partition function Z N of the nonequilibrium ensemble is the total weight of the 2 N distinct configurations of particles and holes. This can be written as since the expansion of (D + E) N yields all possible strings of D and E of length N . Explicit expressions for Z N are given below via Eqs. (24) and (26). The case for general q is more complicated [13]; for completeness we present this in Appendix A.
2.1.1. Phase diagram. For 0 ≤ q < 1, the ASEP exhibits three phases that depend on α, β, q. The transitions between phases are typified by sharp changes (nonanalyticities in the limit N → ∞) in the steady-state particle current J and density profile τ i . See Figure 2 for the phase diagram and plots of the density profile in the three phases.
In detail, the phases are: The slow rate of exit from the lattice, β, creates an accumulation of particles, queueing to leave the system at the right hand boundary, which propagates back into the bulk of the system. Consequently the system is characterised by a high occupation density overall.
The particles enter the system at a low rate α and are able to freely move through the bulk and exit without jamming. Consequently the system is characterised by a low occupation density overall.
• Maximal current (MC), α, β > (1 − q)/2. Here, the current J ∼ (1 − q)/4 is maximised. In this phase, particles enter and leave the system at high rates, and the current is restricted by how efficiently the particles can move through the bulk.
We refer the reader to [9] for further discussion of the phase diagram. Notice that in the limit q → 1 the phase boundaries approach the axes. Consequently, the SSEP exhibits no phase transitions. We also mention the reverse bias case q > 1. The reduction relations (2-4) still hold in this phase (and as such the q-general mappings to be presented will still hold), however like the SSEP there are no phase transitions [13].
Finally, we have a particle-hole symmetry: the dynamics of particles moving to the right are identical to the dynamics of holes moving to the left [9]. Thus the system of particles moving through the lattice from left to right with entry rate α and exit rate β is identical to a system of holes moving from right to left with entry rate β and exit rate α. Then the high density phase is the counterpart, for holes, of the low density phase and the low density phase is the counterpart, for holes, of the high density phase.

Explicit matrix representation
In the previous section, we saw how to use use the reduction relations (2-5) in a formal way to calculate configurational weights-that is, without reference to any matrix explicit representation. One can go on to calculate physical observables, such as the current and density profile, in this way [7]. However, for the purposes of identifying mappings to combinatorial enumeration problems, it is often helpful to write out D, E, W |, |V explicitly. Generally, no finite-size matrices obey (2-5) (except along special parameter curves [14,15]) and we instead resort to semi-infinite representations, of which several are known. For our purposes, the most useful representation is [9,13,16]: with the boundary vectors κ is a normalisation factor defined such that W |V = 1. Here we have employed the shorthand Although the objects in (12)(13)(14)(15) diverge as q → 1, any product of the form lim q→1 W |DEDDE . . . |V is well defined, and the limit gives the correct SSEP configuration weight. Other representations are possible, see [9]. The matrices (12), (13) are reminiscent of ladder operators for the quantum harmonic oscillator. That is, these matrices act on a state ket to transform it into a superposition of |n and |n ± 1 . In fact, these matrices can be related to the ladder operators for a q-deformed quantum harmonic oscillator, a fact that was exploited in [13,16] to calculate physical properties for the PASEP. We now make three further observations of the ASEP, all of which motivate a combinatorial interpretation of the matrix product representation.

Three observations
2.3.1. Combinatorial factors in α = β = 1 partition functions. Beginning with the SSEP, we first directly calculate the partition function Z N . Several methods are known [17,18], and here we present perhaps the most straightforward of these, following [19]. The partition function is written where C = D + E. Now, the commutation property (2) for q = 1, implies that DC = C(D + 1). Repeating this (N − 1) times gives which we insert into (18) and apply the reduction relations Given that Z 0 = 1, this recursion is easily solved to give In the case α = β = 1, (24) reduces to The factorial is the most familiar combinatorial number, counting the number of permutations of (N + 1) integers in this case. This suggests that the SSEP with α = β = 1 can be related to a uniform distribution over the space of permutations. In Section 4 we will see that this is the case.
We now turn to the TASEP (the case q = 0). Although the reduction relation (2) simplifies when setting q = 0, evaluation of the partition function proves more challenging in the case of general α and β. Nonetheless, it can be calculated directly using the reduction relations [7] or by formal expansion of its generating function [9,[20][21][22] (see Section 6.2 for an example). One finds This time, setting α = β = 1 reduces (26) to the summation where we have used the Chu-Vandermonde identity [23] ∞ Finally we obtain where C N +1 is the (N + 1) th Catalan number. These numbers are very well-known in combinatorics, solving at least 60 counting problems [24]. For example, C n is the number of ways to match n pairs of brackets, accounting for all the different ways they may be nested. In Appendix B we see that one can obtain (29) rather directly from the form of matrices and vectors. To summarise, the results (25), (29) give two integer sequences (N + 1)! = 1, 2, 6, 24, 120, 720, . . . and C N +1 = 1, 2, 5, 14, 42, 132, . . . , ubiquitous in enumerative combinatorics and now arising in this nonequilibrium physics problem.
2.3.2. TASEP partition function in terms of bicoloured Motzkin and Dyck paths. One way to show a connection between the TASEP and enumeration problems is to examine the explicit matrix and vector representations (12)(13)(14)(15). Here, we use this approach to show a link to counting bicoloured Motzkin paths and Dyck paths [21,25,26].
First, we denote a ladder operator as g, and state ket vectors |n such that g|n = |n − 1 , g † |n = |n + 1 , with boundary condition | − 1 = 0. We also define the scalar product with a bra vector n|k = δ nk . From this, and on setting α = β = 1, the representations (12-15) simplify to and the partition function is an enumeration of all walks in the nonnegative plane that start and return to the origin, with N steps from the set { , , ×, ·}, where '×' and '·' are distinct nonmovement steps. These are bicoloured Motzkin paths. Each bicoloured Motzkin path of length N is equivalent to a Dyck path of length 2(N + 1) and vice versa. Dyck paths comprise only up and down unit steps, starting and ending at the origin without going below zero. To transform a Motzkin path to a Dyck path we associate to: • each × an up step followed by a down step ( , ), • each · a down step followed by an up step ( , ), • each two down steps ( , ), • each two up steps ( , ), and finally bookend each walk an with up and down step. See Figure 3 for an example. The number of Dyck paths of length 2n is known to be the Catalan number C n [24,27]. To see this note that the total number of paths that start and terminate at zero is 2n n . The number of invalid paths-paths that cross below zero-is counted by reflecting these paths about the axis at the point they first hit −1. These reflected paths all terminate at −2 ( Figure 4), and the total number of such paths is 2n n−1 . The number of valid paths is then 2n n − 2n n−1 = 1 n+1 2n n = C n . For the case α = β = 1, each of these paths are equally weighted. Therefore the normalisation, Z N just counts the total number of paths and is equal to C N +1 , consistent with Eq. (29). See Appendix B or [28] for more details.
which is sometimes referred to as a ballot number [25,29,30], and is the solution to the following enumerative problem: the number of Dyck paths that can be drawn of length 2N that return to the origin p times (including the final return). An example of this is shown in Figure 5. Now, for each of these walks with p returns, we create a set of (p+1) walks whereby the walk is inverted about zero at the q th return, taking q = (0, 1, . . . p), again see Figure  5. Finally we associate to each of these new inverted walks a factor of (1/α) q (1/β) p−q . By this construction, these walks can return to the origin multiple times, but cross it at most once. Such walks have been considered in the context of the TASEP in [25] and is called a one-transit walk. A weight 1/α is applied to each return from above, and 1/β to each return from below. Summing the weights over all such walks then gives the TASEP partition function in (26) [21,25].
In this picture, we see very clearly the connection to an equilibrium partition function over an extended configuration space. Recall that in the TASEP, there are 2 N configurations of particles and holes. The corresponding set of one-transit walks contains C N +1 configurations, which exceeds 2 N : asymptotically, C n ∼ 4 n √ πn 3/2 . Each walk has a weight that can be interpreted as a Boltzmann factor; rewritingα = ln α,β = ln β, the weight for a walk with given p and q can be written as e −qα−(p−q)β . Summing over multiple such Boltzmann-like weights gives the TASEP partition function (26).
As we further discuss in Section 3, the mapping from TASEP configurations to one-transit walks and other combinatorial objects is one-to-many. That is, while each walk can be uniquely identified with a TASEP configuration, the converse is not true. Another way to look at this is as the TASEP defining a partitioning of an extended configuration space. The partition function is invariant under this partitioning; however other measures, such as entropies, are sensitive to it (see Section 6).

Structure of this review
The three observations in the previous subsection point to a deeper underlying combinatorial structure of the matrix product solution to the ASEP stationary state. In what follows we formalise and develop equivalent combinatorial interpretations of the matrix product weights. See Figure 6 for a schematic illustration of the mappings between these interpretations.
As implied in Section 2.3.1, we see the combinatorial structure is at its clearest in the α = β = 1 TASEP where Catalan numbers arise. The case α = β = 1, q = 0 is the focus of Section 3, where we show a mapping between nonequilibrium configurations and path enumeration problems. In Section 4 we discuss a combinatorial problem of permutations that the SSEP (q = 1) maps to. In Section 5 we show how these mappings generalise to the full α, β, q parameter space. It turns out one can associate a q-dependent weight to the permutations of Section 4, thus generalising to the PASEP. It then remains to encode the other two parameters α, β into these mappings, which we discuss in Section 5.4.1.
In Section 6 we consider a more general calculation in the TASEP, of a quantity H λ known as the Rényi entropy. This is a measure of the partitioning of the state space, and e H λ is in turn a measure of an effective number of participating configurations in the NESS. For the case λ = 2 we show how the equivalence to path enumeration problems maps the calculation onto another combinatorial problem, involving enumerating random walks in the upper quadrant. We may then use results and techniques from the random walk literature to compute the Rényi entropy in this case.
Finally, in Section 7 we summarise generalisations of the matrix product solution to the case of multispecies ASEPs such as those involving second-class particles. Our aim is to point out further interesting combinatorial connections such as priority queues and connections with integrable systems and algebraic combinatorics that are currently being explored. We conclude in Section 8. Figure 6. Schematic of the combinatorial mappings to be outlined in this review. An ASEP configuration (an example is illustrated in the left column) has a one-tomany mapping to certain dominated paths (illustrated in the middle column), which we propose in Section 5.5 to in turn map one-to-many to permutations of integers that follow certain rules (illustrated in the right column).

α = β = 1 TASEP
The α = β = 1 TASEP proves to be the most analytically tractable system as the weights of configurations are integers. We outline a one-to-one mapping between configurations of the TASEP and a class of length-N paths, and introduce a measure of dominance [31][32][33] to find the weight of the configuration. This mapping proves equivalent to several others, including that to Motzkin paths, which arises naturally from the explicit matrix representation (12)(13)(14)(15) in the discussion of Section 2.3.2. We frame the state space in terms of the path dominance mapping, as the translation from TASEP configuration to path is simple in this case and offers an intuitive link between the weight of a configuration and the area its path encloses.  (35), (36) and (37). Right: two paths T , T . Here, T dominates T .

Mapping to a path dominance problem
Consider the set of discrete paths T ∈ {↑, →} N that begin at (0, 0), and end at (Q, P ), with P + Q = N steps in total. The total number of paths is A path T can be defined by its set of steps. For example, the path shown in Figure 7, left, can be specified as Alternatively, we can specify, for each value of the x-coordinate 0, 1, . . . Q the maximal y-coordinate of the path. For the path T above, we would have Equally, we could specify the maximal x-coordinate for each value of the y-coordinate 0, 1, . . . P : With this formalism established, we can now define what is meant by dominance [31][32][33]. Take two paths T , T which both terminate at (Q, P ). T dominates T (denoted T T ) if T lies completely on or below T (see Figure 7, right). In terms of the maximal x and y coordinates, for all i. By this definition, T dominates itself, and it is also possible that for two paths, neither dominates the other; if the paths cross then neither path lies completely under the perimeter of the other. We emphasize that this formalism only applies to paths of the same length that have the same start and end points.
This leads to the following combinatorial problem: how many paths W(T ) in total does T dominate? This quantity can be written out iteratively, accumulating all possible dominated paths as T grows step by step. Formally, this is which is a set of Q nested sums. Take, for example, T = (↑, ↑, →, →), y T = (2, 2, 2). This has a weight of 6, found by manually drawing all dominated paths ( Figure 8), or from the summation This problem is of interest to us as each path of length N maps uniquely to a length-N ASEP configuration. Specifically, an ASEP configuration with occupied sites (j 1 , j 2 , . . . j P ) maps to a path T where steps (j 1 , j 2 , . . . j P ) are ↑, and the remaining steps are →. In other words, for a given dominant path T , we can read off the TASEP configuration by going along the path and translating each upward step to a particle, and each rightward step to a hole.
The weight of a TASEP configuration then turns out to be given by W(T ), the number of paths that T dominates. For brevity, we will also refer to W(T ) as the weight of the path. For example, the path in Figure 8 maps to C = (1, 1, 0, 0), which indeed has a weight of 6: This mapping is proven by showing that weight of a path (38) satisfies a set of reduction relations equivalent to (3)(4)(5)(6). More formally, we require where the notation W(a, b, . . .) denotes concatenation of the path segments a, b, . . . . Eqs (44) and (45) are equivalent to W |E = W | and D|V = |V respectively, and  Figure 9. Graphical representation of (44), (45). Adding a → to the start or a ↑ to the end of a path does not change its weight (i.e., the number of paths it can dominate). are trivial by inspection ( Figure 9). Relation (46) is the equivalent of DE = D + E (see Figure 10) and requires more work, but can be derived by brute force from the summation formula (38), see Appendix C. We now highlight three results that first originated in the path dominance literature and that we can exploit to give insights into the TASEP without any additional work. A fourth result concerns the λ = 2 Rényi entropy, which we reserve for Section 6.3.1.
3.1.1. Most probable configuration. The first result is a simple observation, and is that for a length-N path containing P ↑ steps, the most dominant path is as this rectangular path encloses all others. The equivalent TASEP configuration C = (1, . . . 1, 0, . . . 0) is P particles stacked to the left, and is the most probable configuration with P particles. Furthermore, the most probable configuration overall will be N/2 particles followed by N/2 holes (if N is odd, the N/2 and N/2particle configurations are equally most probable). In the matrix formalism, this weight corresponds to the decomposition of the string W |D P E N −P |V using the algebraic rules.
At the other extreme, any configuration with P particles stacked to the right has the minimum weight of 1. This is because the only path that T * = (→, . . . →, ↑, . . . ↑) dominates is itself, though in the matrix product formulation this is already trivial given W |E . . . ED . . . D|V = 1.
3.1.2. Weight with fixed particle number and Narayana numbers. Given this mapping, the total weight of configurations C P with P particles is the total weight of all paths that terminate at (N − P, P ). In the path dominance literature this is known [33]: where T (n, k) is a Narayana number [34] ( Table 1, sequence A001263 in the OEIS [35]). One combinatorial interpretation of these numbers is the number of ways to match n pairs of brackets with k innermost pairs. That is, the sequence (()()) would be counted by T (3,2). Recall that the total number of ways to match n pairs of brackets (without the restriction on innermost pairs) is the Catalan number C n . Consequently n k=1 T (n, k) = C n , and we find from (48) that as previously. In Appendix B we provide a derivation of these results within the matrix product formalism.

Determinant form of configuration weight.
Finally, and most significantly, Narayana [34] (and later Kreweras [32]) has shown in this path dominance problem that the weight of a path can be written as a determinant: or equivalently ('turning the path on its side') With the mapping from paths, this in turn provides an analytic formula for the weight of any TASEP configuration. For example, recall the example path from Eq. (35), Figure  7. Using the first determinant formula, this has weight from its y-coordinates (36) and equivalently using its x-coordinates (37) This path T maps to the TASEP configuration which then implies that the determinants in (53) and (54) are equivalent to the matrix product This reveals a deeper link between the matrix product approach and the reduction relations (3-6) with an elegant determinant structure. Probing these determinants further, notice from this example that reading down each column reveals the (y m−1 +1) th row in Pascal's triangle. It is also 'nearly' a lower-diagonal matrix, and a simple example of a Hessenberg matrix. Taking (53), this gives a simplified recursive determinant formula (adapted from Theorem 2.1 of [36]): where M (r−1) is the (r − 1) th minor of M .
In the context of the TASEP, this determinant formula has since been improved upon to encode α and β, see Section 5.2 and [37].

Other representations
We refer to an ordered pair of two paths where one dominates the other as a dominated path. As previously noted in Section 3.1.2, the total number of dominated paths is given by the Catalan number C N +1 . This set of dominated paths is the extended configuration space. This space can be equivalently expressed in terms of bicoloured Motzkin paths or "complete configurations", which we now discuss.
The full partition function Z N is the number of unique dominated paths. Consider one such configuration with two paths T T . Denote the i th steps of T , T as T (i), T (i) respectively.
Comparing the two paths, on each step we have four possible outcomes, which we track with a height difference h ≥ 0, that must start and end at zero: The paths run parallel vertically, ∆h = 0, • T (i) = T (i) = →. The paths run parallel horizontally, ∆h = 0.
Over each step, h can therefore change by ±1, or zero in two distinct ways (denoted with '·' and '×'). The partition function is then equivalently the number of paths moving left-to-right of length N , from the step set { , , · , ×}, that start and end at zero, without going below zero (as T T ). This is then an enumeration of bicoloured Motzkin paths.
Extending this idea, the weight of a length-N configuration C with sites (j 1 , j 2 , . . . j P ) occupied is the number of length-N bicoloured Motzkin paths that can be drawn from { , ·} at steps (j 1 , j 2 , . . . j P ), and { , ×} in the remaining steps. Thus each Motzkin path is mapped one-to-one to a dominated path. See Figure 11 for an example, where we draw all N = 3 ASEP configurations in terms of Motzkin paths.
This Motzkin path interpretation aligns neatly with the explicit representation we quote in (12)(13)(14)(15). For other explicit representations, other path interpretations naturally arise. Brak et al. present a comprehensive set of these alternative walks in [38], as well as encoding weights to generalise for α, β, q.
3.2.2. Markov chain of "complete configurations". Duchi and Schaeffer [39] express this same space of C N +1 configurations as a set of closed, two-row systems, which they term complete configurations. Furthermore, they define a Markov process in this space that reproduces ASEP dynamics on the top row of the system.
Each of these complete configurations comprise N particles and N holes (which they refer to as "black" and "white" particles), arranged across two rows of length N . The particles may be arranged in any way across both rows, given the constraint that there are always at least as many particles as holes in the first i columns, i = 1, 2, . . . N (the "positivity condition"). The top-row configuration is the ASEP state that the complete configurations map to.
The Markov process that the authors construct is a clockwise flow of these N particles around both rows, with ASEP-like hopping on the top row, and a set of bottom-row dynamics (involving long-range "sweeps" of clusters of particles or holes) so as to preserve the positivity condition. The top row of these closed configurations replicate open TASEP dynamics. In particular, we note that a feature of a complete configuration is that if the top-left site is empty, the bottom-left row is occupiedotherwise the positivity condition would be violated. This means that a particle can always enter the top row at a rate α, just as in the TASEP. Similarly, if the top-right site is occupied, the bottom-right site must be empty, allowing particles to exit the top row at rate β, again as in the TASEP.
In the general α, β case, the bottom row dynamics that yields the desired weighting Figure 12. Simplified two-row dynamics for the case α = β = 1 which generates a uniform distribution over the space of complete configurations. Particles hop clockwise into empty spaces around the lattice at unit rate. Vertical dashed lines indicate zone boundaries, wherein each zone contains an equal number of particles and holes. If, on the upper row, a particle crosses a zone boundary when it hops (shown by the curved blue arrow), the particle below the receiving site is also forced to hop (shown by the straight blue arrow).
of configurations in the extended space is rather complex. However, there is a simplified dynamics, involving only local moves, that generates a uniform distribution over the extended set of complete configurations in the special case where α = β = 1. These dynamics, which we derive in Appendix D, are illustrated in Figure 12. Each particle hops at unit rate, as long as the site in front of them is empty. As noted above, the positivity condition already guarantees that particles enter and leave the top row at unit rate.
The one nontrivial aspect of the dynamics involves the notion of zones. Going from left to right across the lattice, we mark a zone boundary at each point where the total number of particles to the left of the zone boundary (on both rows, and all the way back to the left boundary) is equal to the total number of holes. These are shown with vertical dashed lines in Figure 12. Now, if a particle on the top row can hop across a zone boundary, then the positivity condition implies that the site below must be empty, and that the site to the right on the bottom row must contain a particle. If the particle on the top row were to hop, then the positivity condition would be violated. Thus, to preserve positivity, the bottom row particle to the right must hop to the left when the particle on the top row hops. This forced bottom row move is shown with a straight arrow in Figure 12. Note that if we cover up the bottom row, we obtain the usual dynamics of the TASEP. However, if we cover up the top row, we get a different dynamics as a bottom row particle may move at double the usual rate, depending on the configuration of the top row.
By construction, the dynamics in the full system remain within the space of complete configurations. It can also be shown that every complete configuration can be reached from any other (see Appendix D). This means that all complete configurations have a nonzero probability in the steady state. Finally, the number of ways out of each configuration is equal to the number of ways in (see again Appendix D). Together, these three properties imply that the distribution over complete configurations is uniform. One can also show that these dynamics are dynamically reversible [40], that is, they Figure 13. Example of a complete configuration in [39] (left) and its equivalent Motzkin path (centre) and dominated path (right). The top row of the complete configuration shows that these correspond to C = (1, 0, 1, 0, 0, 0, 1). satisfy a generalisation of detailed balance (see Appendix D). When we sum over the bottom row configurations for any given top row, we obtain the weight of a physical configuration of the TASEP.
The basic principles that apply to this simplified dynamics can also be shown for the general α, β case [39]. In particular, one can uniquely associate to each way into a configuration a way to leave it. Not all of the moves take place at the same rate when α or β is not equal to 1; consequently complete configurations have different weights in this case (see Section 5).
Here we expand on how each complete configuration in this two-row system maps to a Motzkin path or dominated path. With reference to Figure 13, if we assign to each column with (τ top ; τ bottom ) entries a for (1; 1), a for (0; 0), a '·' for (1; 0) and '×' for (0; 1), then configurations are once again a set of bicoloured Motzkin paths once the positivity condition is imposed. Corteel and Williams [41] have since introduced a Markov chain that reproduces PASEP dynamics (where the additional parameter q is introduced), using a larger set of (N +1)! configurations. The broader idea of connecting stochastic transport processes to simpler "dual" systems has also been applied to analyse models away from the ASEP, see [42].
Our discussion so far has been limited to the TASEP (q = 0). We now move from the totally asymmetric case to the totally symmetric case where particles can hop either direction in the bulk at equal rates, by setting q = 1.
We previously showed that the SSEP partition function in the case α = β = 1 is Z N = (N + 1)!, see Eq. (25). This combined with the analysis of the TASEP in Section 3 suggests that the 2 N configurations of the SSEP may map to an even larger set of (N + 1)! > C N +1 > 2 N configurations. This indeed turns out to be the case; consider the integers (0, 1, . . . N ), of which there are (N + 1)! permutations. The 2 N configurations of the SSEP define a partitioning of these (N + 1)! permutations.

Mapping to a permutation problem
This mapping was first identified and formally proven by Corteel and Williams [41] in the context of a Markov chain of permutations. Here we focus only on the mapping from SSEP to permutations using a slightly different but equivalent formalism to [41]. This is complemented with a more detailed analysis in Appendix E.
Consider a permutation of the (N + 1) integers in (0, 1, . . . N ), denoted (i 1 , i 2 , . . . i N +1 ). Reading this string of integers left-to-right, we say that i n has been raised by i n+1 if i n+1 > i n . This time, we are interested in the following problem: how many permutations are there where only a particular set of integers (j 1 , j 2 , . . . j P ) are raised?
This proves to be equivalent to the weight of a length-N SSEP configuration with particles at sites (j 1 , j 2 , . . . j P ) + 1. We illustrate this with an example. The SSEP configuration C = (0, 1, 0, 0), has N = 4 sites, and P = 1 particle at position j 1 + 1 = 2. This has a weight of 7, calculated directly with the reduction relations (2)(3)(4)(5): As anticipated, there are also 7 permutations of (0, 1, 2, 3, 4) where only j 1 = 1 is raised (the underline highlights where an integer has been raised): If the SSEP indeed maps to these permutations, we should expect to find an equivalent set of reduction relations like that of the SSEP (46), which we show in Appendix E, in fact for the more general DE = qED + D + E, where weights as powers of q are associated to each permutation (see Section 5).
Having established this mapping, we can quickly derive the steady-state density profile and arbitrary-order correlations between sites. We also use a result in the literature on the combinatorics of permutations, which allows us to find the probability of finding P particles in the system. 4.1.1. Steady-state density profile. We can now identify the average steady-state occupation of site i as being the fraction of permutations of (0, 1, . . . N ) where integer (i − 1) is raised. Note that we do not care whether any other integers are raised. One slight complication is that (i − 1) can only be raised if it is not at the final position within the permutation. From this interpretation we can very quickly calculate the full density profile. If (i − 1) is not in the final position, (i − 1) can be raised by any of (i, i + 1, . . . N ) from the N integers greater than (i − 1), giving fraction N −(i−1) N . We then multiply by the fraction of permutations where i is not in the final position which is N N +1 . We thus obtain recovering the known linear profile [43].
4.1.2. Arbitrary-order correlation functions. We can extend this approach to calculate higher-order correlations between different sites without having to perform any explicit matrix calculation. First, consider the correlation where i 2 > i 1 . This is equivalently the fraction of permutations of (0, 1, . . . N ) where both (i 1 − 1) and (i 2 − 1) are raised. First, (i 2 − 1) can be raised by any of the (N + 1 − i 2 ) integers from (i 2 , . . . N ), and the fraction of suitable permutations is then (N + 1 − i 2 )/(N + 1). In this subset, (i 1 − 1) can be raised by any of the (N + 1 − i 1 ) integers from (i 1 , . . . N ), excluding the integer that raised (i 2 − 1). The fraction of valid permutations here is then (N − i 1 )/N . Combined, we then recover the result from [17,43] By the same interpretation this can be extended to an arbitrary-order correlation between K different sites i K , i K−1 , . . . i 1 , where i K > i K−1 > . . . > i 1 [18,44]: 4.1.3. Weight with fixed particle number and Eulerian numbers. The sum of all weights of configurations C P with P particles is the number of permutations of (0, 1, . . . N ) with a total of P integers raised (again we do not care which integers in particular). We state the result from the combinatorial literature [45,46]: where is known as an Eulerian number (Table 2, sequence A008292 in the OEIS [47]), and has several neat properties reminiscent of binomial coefficients, such as the recursion [48] n + 1 k = (n + 1 − k) n k − 1 + (k + 1) n k . The generating function of (70) is succinct [48]: where the coefficient {t N z P }G(t, z) is the probability of finding P particles in a length-N SSEP. Finally, the summation over Eulerian numbers for fixed N is equivalent to the summation of all N -site SSEP weights, and gives the factorial [49] N P =0 This is trivial in the context of Eulerian numbers, as it is simply the summation of all permutations of (N + 1) integers.

Generalised parameter mappings
Up to now, we have focused on the parameter restriction α = β = 1, q = 0, 1. To generalise for α, β, q, we do not need to expand beyond the state spaces of dominated paths and permutations already introduced, however we now associate weights, as products of α, β, q. The highlight of the following is the pleasing result that a closedform formula is available for the weight of a general TASEP configuration in the form of a matrix determinant. Finally in Section 5.5 we propose a mapping for general q that would effectively interpolate between bicoloured Motzkin paths and weighted permutations, and in turn the weights of TASEP, PASEP and SSEP configurations.

α = β = 1 PASEP and weighted permutations
Let us recall the exact expression for the PASEP partition function [13] which is included in Appendix A. This expression interpolates between q = 1 and q = 0 suggesting that there may be some combinatorial entities which interpolate between the mappings we have identified in Sections 3 and 4. We first expand on the results of Section 4.1 by showing how an arbitrary q is encoded into the mapping of SSEP configurations to permutations, first shown in [50]. We remain with our slightly different formalism introduced earlier, again referring the reader to Appendix E for further details of the mapping.
Take a permutation of the (N + 1) integers (0, 1, . . . N ) in which a set of integers (j 1 , j 2 , . . . j P ) in the permutation are raised. Let k be the integer that raises j . For example, if the permutation contains the integer 1 immediately followed by 4, then we will have j = 1 and k = 4 for some = 1, 2, . . . P . Now, for a given raise, we look for the integers with values between j and k and that sit to the right of the pair j k in the permutation. Let the number of such integers be r . Then, to each raise we associate a weight q r . The total weight of the permutation is the product of these weights.

Determinant form of TASEP weight with general α, β
Away from permutations and returning to the TASEP, Mandelshtam has generalised the determinant form (52) of a TASEP configuration weight for arbitrary α, β (Corollary 5.2 in [37], modified to be consistent with notation used here): where the entries of M with n, m = 1, . . . P , and the x n , x m are the coordinates associated to an ASEP configuration in Section 3.1.
Using this formula, we are able to write down an expression for the TASEP partition function Z N in the form of a determinant. As far as we are aware the expression we now derive has not previously appeared in the literature. Given that the partition function is the weight of a single 'staircase' path of length 2N (see Figure  14). For this path, x m = m, and Mandelshtam's formula (75), (76) eventually reduces to where we see rows of Pascal's triangle in the coefficients of 1, [1/α + 1/β], 1/αβ when reading down columns of M . We verify in Appendix F that (78) and the partition function are equivalent, through generating functions.

α, β generalisation of path dominance problem
Following on from this determinant formula, there is a straightforward generalisation to α, β in the dominated path interpretation of TASEP weights. In the context of the original reference [37] these are referred to as "weighted Catalan paths", which translate into our formalism as follows: each dominated path has an associated weight (1/α) p (1/β) q , where p is the number of horizontal steps where both paths run together, and q is the number of 'up' steps the dominated path takes at the end of the walk. See

General α, β, q
We have arrived at the most general case of general α, β, q. We shall discuss a natural interpretation in terms of bicoloured Motzkin paths that arises from an explicit matrix representation. Otherwise, the most notable progress here has been by Corteel and Williams [50], who derive a generalised version of the path representation of configurations in Section 3.1, termed permutation tableaux. At this level of generality, there are few new physical insights that have been made other than establishment of the mapping.
The weight of the path is then the product of these weights. See Figure 16 for an example. Note that and always appear in pairs, eliminating the square root in factors of √ c n .

Mapping between Motzkin paths and permutations for α = β = 1 and general q.
We conclude this section by proposing a one-to-one mapping between a set of decorated bicoloured Motzkin paths and a permutation of (N + 1) integers that applies for general q when α = β = 1. These permutations have the integers (j 1 , j 2 , . . . j P ) raised by the integers (k 1 , k 2 , . . . k P ), and have the a weight q r as described in Section 5.1. In addition to the usual up steps ( ), down steps ( ) and horizontal steps of two colours ( · and ×), the decorated bicoloured Motzkin path features at each horizontal position i = 0, 1, . . . N a bauble hanging at a height m i = 0, 1, . . . n i , where n i is the height of the path at position i. An example of such a path is shown in Figure 17. The weight of the path is q r where r = N i=0 m i , i.e., the sum of the bauble heights. Summing over all bauble positions for a given Motzkin path yields the weight of Section 5.4.1 after putting α = β = 1 (and therewith a = b = −q). To see this we note that each . . . pair between heights n and (n + 1) contributes a weight (1 + q + . . . + q n )(1 + q + . . . + q n+1 ) while a × and a · at height n each contribute a weight (1 + q + . . . + q n ). Each possible decorated bicoloured Motzkin path can be translated into a permutation of the integers (0, 1, . . . N ) through an iterative procedure that we now describe. The basic idea is that the horizontal position at the start of each up step or horizontal step of type · determines the integers j that are raised, and the horizontal position at the end of each down step or horizontal step of type · determines the integers k that do the raising. The baubles determine the order that integers appear in the permutation. The algorithm ensures that the constraint that only the integers j are raised is never violated.
The iterative procedure begins with an empty string, and at step i involves placing the integer i into the string. At intermediate steps, the string may contain placeholder elements which receive the integers k that do the raising. For each successive integer i = 0, 1, . . . N , the iteration comprises two sub-steps: (a) If the previous segment (that connecting i − 1 to i) is a down-step or a horizontal step of type · , replace the (m i + 1) th placeholder from the left with the integer i. Otherwise, if m i > 0, insert integer i after the m th i placeholder. Otherwise, place integer i at the start of the string. (Note: when i = 0 the only option is for the string to become 0.) (b) If the next segment (that connecting i to i + 1) is an up-step or a horizontal step of type · , insert a placeholder element after the integer i. (Note: when i = N , there is no next segment, and the algorithm terminates.) In Table 3 we illustrate this procedure with the example path and decoration shown in Figure 17, thereby determining that it corresponds to the permutation 6 1 4 3 0 2 5. Since each integer is added to the string exactly once, we must end up with some permutation of the (N + 1) integers (0, 1, . . . N ). Sub-step (b) of the algorithm ensures that a placeholder is inserted immediately after any integer that must be raised. The first clause of sub-step (a) ensures that those integers that do the raising are inserted into the placeholders: these integers must necessarily be larger than those to the left of the placeholder, and so these are indeed raised as required. The remaining clauses of sub-step (a) insert the non-raising integers: these are always entered either at the start of the string (where no raising is possible) or after a placeholder element that will receive a higher integer at a later stage of the algorithm (and therefore also do not raise). Thus the resulting string has integers (j 1 , j 2 . . . j P ) raised, and integers (k 1 , k 2 , . . . k P )  doing the raising, as required. The height of the path keeps track of the number of placeholders and therefore the number of places where each successive integer can be placed while respecting the constraints. The position of the bauble determines how many placeholders lie to the left of the inserted integer, and therefore the power of q that contributes to the weight. To show that the decorated paths map one-to-one to permutations, one needs to show that each permutation obtained from the above algorithm is distinct. For two paths of the same shape, but different bauble positions, this seems plausible, because the bauble position determines the position of each new integer i relative to those already present in the string. The relative order of the first i integers can not be changed by later insertions, so one expects each bauble configuration to map to a distinct permutation. Meanwhile, different permutations correspond to different sets of raised or raising integers, and therewith a partitioning of the permutations into distinct subsets. Finally, it is already established in (25) that lim q→1 Z(α, β, q) = (N + 1)! for α = β = 1. Since this normalisation counts the total number of decorated bicoloured Motzkin paths, it would then follow that every possible permutation is represented by one of the paths.
A formal proof of the proposition that there is a one-to-one mapping between decorated Motzkin paths and permutations would be welcome. In particular, it would demonstrate the one-to-many mapping from ASEP configurations to bicoloured Motzkin paths which themselves map one-to-many to permutations as was illustrated in Figure 6. For each bicoloured Motzkin path, exactly one mapped permutation has weight q 0 = 1, that is, the decorated path with baubles at m = (0, 0, . . . 0)), while all others have a positive power of q. This would connect the path dominance mapping of the TASEP (q → 0), and the permutation mapping of the SSEP (q → 1). [50,54,55] have mapped the most general case of α, β, q to a problem in an area known as tableaux combinatorics. We refer the reader to [50] for the original work, and [54,55] for a more generalised case of staircase tableaux that encodes two extra parameters γ, δ (so particles may also enter from the right, and leave from the left).

Permutation and staircase tableaux. Elsewhere, Corteel and Williams
The details are beyond the scope of this work, but to sketch their approach the authors take ASEP configurations as paths drawn in Section 3.1, and construct a grid across the area the path bounds (a Young diagram). Each entry of this grid can take a value of α, β, q or 1 (or a generalised hop-right rate u), with a set of rules as to which values can go where. The weight of this permutation tableaux is then the product of all of the entries. For a given ASEP configuration a set of these permutation tableaux can be drawn, and the weight of the configuration is the sum of weights of these tableaux. Further combinatorial interpretations of the PASEP partition function have been obtained by Josuat-Vergès [56]: as the generating function for weighted Motzkin paths known as Laguerre histories, and as the generating function for permutations with respect to maxima and minima and other features.

State space measures of the TASEP
We now shift our focus to functions that involve all weights of the TASEP. This is a different type of problem to the mappings discussed up to now, as we are now calculating measures of the state space as a whole, instead of individual configurations. We will still use the random walk interpretation offered by the explicit matrix representations, but now in higher dimensions.
By exactly solving the TASEP and its partition function, we can calculate steadystate quantities such as the particle current J, the steady-state density profile τ i and higher-order correlations between different sites τ i τ j [6]. However, in the absence of equilibrium statistical mechanics, how do we probe the finer details of the probability distribution?

Rényi entropy
To this end, we introduce the Rényi entropy, which is defined [57] as a measure of a full probability distribution. λ is a nonnegative number. To provide an interpretation of H λ , consider A system with L equally likely configurations has e H λ = L. At the other extreme, the same system with a single configuration with probability one has e H λ = 1. Between these extremes, then, e H λ measures an effective number of participating configurations. By increasing λ, the measure places more weight on the higher P values, as the lower ones are exponentially suppressed. At the extreme, e H 0 is the number of configurations with nonzero probability.
The λ → 1 limit recovers the familiar Shannon entropy: For all λ the Rényi entropy (86) is exactly calculable for any equilibrium system with a known partition function, as (aside from the special λ = 1 case) a ratio of partition functions at different temperatures [58]. This includes the one-transit walk in Section 2.3.3, which is an equilibrium system. We have seen from these combinatorial mappings that weights of the TASEP do not take such a convenient exponential form (e.g. Equation (7)) and are summations that can not generally be factorised. Furthermore, even though the ASEP shares a partition function with an equilibrium system, any sum of weights to a power λ depends on how the C N +1 equilibrium configurations map onto 2 N nonequilibrium configurations-in other words, the partitioning of the partition function.
6.1.1. Enumeration of weights raised to a power. Given the form of the Rényi entropy in (86), we are interested in calculating for the TASEP, for different values of λ. As an introductory example, we first calculate the partition function (for general α, β), in terms of weighted bicoloured Motzkin paths. This is the straightforward λ = 1 case of (89). We then show that the problem of (89) is -for integer λ -equivalent to a lattice enumeration problem in λ dimensions (the partition function being the one-dimensional problem). The λ = 2 case was previously solved by the authors in [22].

λ = 1; sum of weights
With q = 0, the explicit representation of the matrices D, E and vectors W |, |V in (12-15) reduce to Then the boundary vectors a, a 2 , a 3 , . . .) , recalling a = (1 − α)/α, b = (1 − β)/β (16). From this, the partition functiion Z N is then after writing the scalar product explicitly. The RHS of (93) takes the form of a generating function in a, b of the quantity i| g + g † + 2 N |k . Given the binomial expansion of g + g † + 2 N and that g|0 = 0, this quantity is the number of bicoloured Motzkin paths of length N between coordinates i and k. The partition function therefore follows from the path enumeration problem i|(g + g † + 2) N |k . One way of solving this is by generating functions. We present a neat example of this by Depken [20]: first, the generating function is written Now using relation (6), we find taking the negative root of η(z) to ensure Z(0) = 1. We should expect this problem to simplify in the case α = β = 1 (a = b = 0). Indeed this is the case, and (93) reduces to which is the number of such walks that start and end at zero. As discussed in Section 2.3.2 and Appendix B this is of course the Catalan number C N +1 . With this, we return to the central problem of (89) and Rényi entropies. Configurations of the TASEP map to C N +1 paths. The partition function Z N is then a summation of this larger state space. The summation of TASEP weights each raised to a power, however, is a much more complex problem as it depends specifically on which paths map to which TASEP configurations. We first take the sum of squared weights.

λ = 2: sum of squared TASEP weights
The sum of squared weights of the ASEP is compactly written in tensor product formalism [22]: where A ⊗ B is the tensor product of A and B. Taking the example of N = 2 We write the explicit representation in (90), (91), now in two dimensions: defining the ladder operators g 1 , g 2 : We can now explicitly write the tensor product (98) as and again the RHS resembles a generating function in a, a, b, b for the quantity i|⊗ j|(g 1 +g 2 +g 1 g 2 +g † 1 +g † 2 +g † 1 g † 2 +2) N |k ⊗|l . This is a two dimensional lattice walk; the six ladder terms can be interpreted as the step set {↑, ↓, →, ←, , }, alongside two non-movement steps '×', '·'. (105) is thus a generating function of the number of walks in the upper quadrant of length N between (i, j) and (k, l) from the step set {↑, ↓, →, ←, , , ×, ·}. See Figure 18 for such a walk. In practice the non-movement steps are easily integrated out.
In an earlier work [22] we calculated the generating function of (105) The enumeration of these upper quadrant walks for a general step set is is a wellresearched topic, and a technique known as the kernel method finds the generating function for most walks [59]. The step set we are dealt with here proves to be one of the more stubborn, and is solved by a more involved obstinate kernel method. While the details of the calculation are beyond the scope of this review, the symmetry of the sixstep walk is exploited in order to solve what turns out to be a two-parameter recursion relation of the generating function [59].
Having found an explicit form of (106), the asymptotic scaling of the sum of squared Figure 19. Graphical representation of the squared weight of a path in the path dominance formalism.
weights emerges, with a different scaling for the three phases. After normalising, We briefly mention that the scaling in this maximal current phase implies that the effective number scales as 2 N / √ N , which is the same asymptotic scaling as the binomial coefficient N N/2 . Given the MC phase has a density of τ i ≈ 1/2 in the bulk, the effective number is in turn proportional to the number of half-filled configurations.
6.3.1. α = β = 1 simplification; direct solution. Setting α = β = 1 (a = b = 0), (105) is the enumeration of walks that start and end at the origin only. Here, the generating function is succinct [22]: The coefficients of this series expansion are (sequence A196148 in the OEIS [60]) and the sum of squared weights for configurations with P particles is (110) with the summation dropped (sequence A111910 [61]). This can be proven by demonstrating that this number sequence and (105) have the same generating function. These numbers take a similar form to Narayana numbers (48). It is interesting to note that a path dominance problem closely related to the sum of square weights was solved directly by Kreweras and Niederhausen in [31] outside of the context of the TASEP. The equivalent problem is the enumeration of triples of dominated paths: paths that can be drawn where one path dominates the other two ( Figure 19).

Higher orders
This path enumeration approach to sums of TASEP weights can be generalised to arbitrary integer power. The sum of weights to the λ th power can be written where Given the explicit representation of D and E, this then reduces to a problem of enumerating the λ-dimension walks in the upper orthant, from the 2 λ+1 steps arising from Even in two dimensions, the step set {↑, ↓, →, ←, , } proved one of the more challenging step sets to solve. The enumeration of λ = 3 octant walks is a current area of research [62,63], but for this particular classification of walk in λ = 3 or higher, no analytical techniques are known. This problem can be equivalently posed as a path dominance problem, extending the λ = 2 case in [31]. The sum of TASEP weights to the λ th power is the number of distinct (λ + 1)-tuples of length-N paths that can be drawn, where one path dominates the λ others.

Multispecies models and further combinatorial connections
So far we have considered in detail the matrix product formulation of the stationary state of the open ASEP and its combinatorial ramifications. As we have seen, this formulation leads to remarkably far-reaching connections between matrix products, lattice paths and weighted permutations.
In this section we briefly give an outlook on how other exclusion processes also have matrix product solutions, and how matrix products may be generalised. The matrix product solution may be applied to a variety of multispecies problems as reviewed in [9]. Here we focus on recent progress on the family of models that involve a hierarchy of particle species: first-class particles, second-class particles and so on.
Our aim here is to highlight further connections between matrix product states and combinatorial constructions and indeed queueing theory.

Second-class particles
One of the first generalisations of the matrix product state of the open ASEP was to a system of first and second-class particles on a ring [64]. Here we first consider the partially asymmetric case. In this system first-class (normal) particles hop to the right with rate 1 and to the left with rate q. They also overtake second-class particles to the right with rate 1 and to the left with rate q. Second-class particles hop to the right with rate 1 and to the left with rate q. Thus first-class particles effectively view second-class particles as holes and from the point of view of holes, second-class particles behave in the same way as first-class particles. The dynamics may be summarised as We generalise our earlier single-species notation, and use the variable τ i = 0, 1, 2 which now implies that site i is empty, contains a first-class particle or contains a second-class particle, respectively. Let us denote by C = (τ 1 , . . . τ N ), a configuration of the system. The matrix product solution of the stationary state was formulated by Derrida, Janowsky, Lebowitz and Speer [64]. In the matrix product formulation on a periodic lattice we now express the weight as a trace of a product of matrices X τ i The periodicity of the lattice is reflected in the periodicity of the trace operation on a product of matrices. The normalisation Z = Z(N, P 1 , P 2 ) now depends on the numbers P 1 , P 2 of first and second-class particles respectively. The matrices X τ i are given by that is, as before, we write a matrix D or a matrix E when the site contains a first-class particle or hole respectively, and if the site contains a second-class particle we write a matrix A. The matrices D, E, A obey the algebraic rules The first equation (122) is as before. The new element is the matrix A.
A convenient representation is given by choosing D and E as before in (12) and (13), and A as the diagonal matrix In the totally asymmetric case q = 0 (125) reduces to Note that this projector obeys A 2 = A and due to the form of A, (120) reduces to where B j is a binary string and ω(B j ) is given by where X τ i is either a D or E according to whether τ i is 1 or 0, is the length of the binary string B j and i labels the entries in that string. Thus the weight of a two species configuration factorises around the second-class particles and reduces to a product of weights (128) of open TASEP configurations which can be computed as before using the reduction rule (122). The matrices D, E and A with A given by (125) may also be used in the open boundary case under certain conditions on the boundary rates for entrance and exit of particles [65]. In particular the case of 'semi-permeable' boundaries where the secondclass particles do not exit or enter has been widely studied [66][67][68][69]. More general open boundary versions require generalisation of the matrices D, E, A to tensor operators and certain conditions on boundary rates [70,71] which relate to Zamolodchikov-Faddeev and Ghoshal-Zamolodchikov relations in integrable systems [72].

Combinatorial connections
In the case of second and first-class paricles with q = 0 on the ring, which we refer to as the 2-TASEP, Ferrari, Fontes and Kohayakawa [73] were able to use the factorisation property described above to give a complete probabilistic description of the stationary measure. Subsequently, Angel [74] was able to make a two-line construction, on line 1 of which there are P 1 particles distributed randomly and similarly on line 2 there are P 1 + P 2 particles distributed randomly. Angel [74] showed that uniformly sampling twoline configurations, then generating the associated 2-TASEP configurations, samples 2-TASEP configurations according to their stationary measure. Thus it follows that the matrix representation of the stationary state furnishes a way of counting the two-line configurations which correspond to a 2-TASEP configuration. Angel also proved that the uniform measure on the two-line configurations is stationary under the dynamics.
Ferrari and Martin [75] were able to define a two-line process and show that the stationary distribution for the two-line process is the uniform distribution. Specifically one can associate to each transition out of any two-line configuration a transition into that configuration such that there is bijection in the transition rates i.e. each member of the collection of transition events out of a configuration is paired with exactly one of the collection of transition events that lead in to that configuration. Then it follows that the stationary state has uniform probability for all allowed configurations of the two-line system. This is essentially the same argument that was used to show a uniform distribution over complete configurations for the TASEP with open boundaries described in Section 3.2.2 (see also Appendix D).
Furthermore, in [75], Ferrari and Martin showed that the two-line construction could be interpreted as a queueing system. The discrete 'time' of the queueing process corresponds to lattice position in the exclusion process (and has no relation to the continuous time of the dynamics of the exclusion process).
Finally, in [76] it was shown how the trajectories of the queueing system are counted within the matrix product formulation. The vector |0 corresponds to an initial queue of length 0. At each discrete time step of the queuing system a matrix D corresponds to a service time of the queue and a matrix E to a non-service time. For example, the action of D on a queue of length n > 0, represented by |n , is There are two possible events at the service time, which correspond to the two terms in (129): the first is the service of a new arrival, so that n remains unchanged; the second is a service and no new arrival, so that n decreases by one. The full details of the queue dynamics are presented in [76]. The trajectory of the length of the queue in the queueing process is precisely a Motzkin path, as defined in Section 2.3.2. Thus this completes the cycle of mappings from matrix product to multiline configurations to queue trajectories to Motzkin paths for the case q = 0.
Very recently [77] a generalised queueing construction, in which each potential service is unused with probability q k when the queue-length is k was shown to give a recursive construction of the stationary distribution of multispecies process with asymmetry parameter q.

Multispecies generalisation
The second-class particle problem is generalised in a straightforward way to the multispecies problem where the class (or species label) of a particle is an integer between 1 to N . Thus there are N classes of particle along with holes, which could be considered as the lowest class of particle. In the partially asymmetric case the dynamics is The dynamics (119) may also be generalised to have different rates for different particle exchanges-sometimes this is referred to as an 'inhomogeneous' multispecies system [78,79].
The totally asymmetric case (q = 0) was first constructed by Ferrari and Martin [75] where their queueing interpretation was generalised to the N -species system. This involved the introduction of a system of tandem priority queues with increasing number of classes of customers in each queue. Similarly to the 2-line process corresponding to the 2-TASEP, Ferrari and Martin [75] defined an N -line process in which each dynamical event corresponds precisely to a dynamical event in the N -TASEP. The steady-state measure of the N -line system is just a uniform distribution of particles. The N -TASEP stationary state is found by counting the multiline configurations which map onto particular N -TASEP configurations. This combinatorial task is then performed with a matrix product formulation.
Meanwhile, the matrix product solution for the totally asymmetric case (q = 0) was constructed in [76] and hierarchical structure for increasing N was elucidated. In these solutions the 'matrices' are in fact tensor products of the fundamental matrices δ = D − 1, = D − 1 and A and obey more complicated relations than (122-123), involving an auxiliary set of 'hat' operators. In [76] it was shown how the matrices generate the system of priority queues defined in [75].
To illustrate how the matrix product generalises from the 2-TASEP case to higher numbers of species and what we mean by tensor operators, let us write down the operators required for the 3-TASEP.
These tensor product operators act on state vectors |l m n ≡ |l ⊗ |m ⊗ |n where |l = 0 for l < 0. It may be shown [76] how a matrix product using X 0 , X 1 , X 2 , X 3 defined in (134-137) precisely enumerate the possible trajectories of the state of the tandem queues giving rise to a given configuration (τ 1 , τ 2 , . . . τ N ).
The matrix product solution was generalised to the partially asymmetric case (q > 0) in [80]. A class of open boundary multispecies processes has been shown to have a matrix product solution [81]. Further properties of the matrix product formulation and its hierarchical structure have been explored in [82,83] for the totally asymmetric and partially symmetric cases. Algebraic properties have been used to generate alternative matrix representations of the Ferrari-Martin construction [84,85].
Finally we mention recent work [86], in which it was shown how certain Macdonald polynomials, which are a set of orthogonal polynomials with some remarkable properties, may be expressed as a matrix product. The matrix product involves q-deformed bosonic operators as does the matrix product for ASEP, which we have discussed in this review. In turn the matrix product formula for Macdonald polynomials can be interpreted in terms of lattice paths leading to a combinatorial interpretation. Also, Macdonald polynomials may be used to express the partition function of the multispecies ASEP. Continuing in this vein a matrix product formula for Koornwinder polynomials has been obtained [87] and the multispecies ASEP has been used to derive new combinatorial formulae for Macdonald polynomials [88].

Conclusion
In this work we have reviewed the connection between the stationary weights of configurations in a paradigmatic nonequilibrium statistical mechanical system (the asymmetric simple exclusion process) and combinatorial enumeration problems, such as counting lattice paths. The earliest solutions of the TASEP (the version of the process in which particles can hop only to the right) appealed to recursion relations [5,8] between configurational weights which can be expressed more powerfully in terms of reduction relations for matrices [7], as described in Section 2. Both the application of recursion relations and the reordering of matrices implicitly define some kind of counting problem. However it is not necessarily obvious from the outset what is being counted. The generalisation to partially asymmetric hopping resulted in more complicated counting problems involving the parameter q.
The most straightforward way to relate the matrix product solution to a lattice path enumeration problem is to exploit a representation of the matrices in terms of the identity, and (in general, q-deformed) raising and lowering operators (Section 2.2). A particular configuration of the ASEP can then be related to a set of Motzkin paths, in which the identity, raising and lowering operators generate steps that are either horizontal, rise upwards, or fall down. Since the matrices are semi-infinite, the paths may not fall below the origin. Thus one set of objects that are being enumerated by the ASEP normalisation is the set of all paths subject to this constraint. This in turn yields a connection to the Catalan numbers, which solve a large number of enumeration problems [24].
Perhaps one of the most appealing representations of a configurational weight in the TASEP is in terms of dominated paths (Section 3.1). Here, a configuration of the TASEP is converted to a path on the square lattice by drawing (in sequence) a vertical step for each particle and a horizontal step for each empty site (hole). The number of paths that fall below this dominant path, and that have the same start and end point, then gives the weight of the TASEP configuration when α = β = 1.
Here we see clearly the general phenomenon whereby a configurational weight in the TASEP is given by a sum over a set of objects with simpler weights that live in a larger space. In the specific case of the dominated paths, the larger space is the set of all lattice paths of a fixed length, and the weights are a power of α multiplied by a power of β. As discussed in Section 2.1, we can think of this as a Boltzmann weight, in which α and β are the exponential of energetic contributions associated with specific steps along the paths.
From a practical point of view, the mapping to enumeration problems can expedite the calculation of physical quantities. For example, we saw in Section 3 that once the mapping is established, results from enumerative combinatorics can be used to establish certain quantities more easily than deriving them from scratch via the matrix product solution. In the case of the Rényi entropy (Section 6), it is only by appealing to the lattice walk picture that this has (yet) been calculated, and then only for the case λ = 2. Perhaps there is scope to exploit the mappings further to extend to general λ, or to obtain new results for systems other than the variants of the asymmetric exclusion process that we have considered here. Finally, it would be of interest to identify cases where generalisations of enumeration problems provide solutions to nonequilibrium stochastic dynamical systems, particularly if these correspond to models for which a matrix product solution has not previously been found.
We have also seen for the open TASEP in Section 3.2.2 and for the multispecies TASEP in Section 7 that it is possible to construct a Markov process on an extended configuration space that converges to the Boltzmann-like distribution described above, with the additional feature that the dynamics in the physical region or physical projection is ASEP-like. In the extended space, the dynamics satisfies a dynamic reversibility condition [40] (equivalently a bijection in the transition rates-see Appendix D.1), which is essentially a form of detailed balance. The idea that we can have a reversible dynamics in the extended space, but an irreversible dynamics in only part of it, can be understood from an information-theoretic perspective. The dynamics set out in Section 3.2.2 induces correlations between the two rows of particles. Projecting onto the subset of physical configurations amounts to erasing these correlations, which in turn causes a loss of information. Since information loss is associated with entropy production, the dynamics in the subsystem can have an irreversible character, even though the dynamics in the full system does not.
It would be of interest to apply this idea to systems where matrix product solutions exist and try to construct the extended configuration space and dynamics. For example there are several open boundary multispecies models where matrix product solutions exist but for which combinatorial interpretations are not yet clear [65-67, 69, 71, 81]. An intriguing question the conditions under which any irreversible stochastic process might be obtainable from a reversible dynamics on a larger space by projecting onto the physical space.
Most generally, it is well-known that the stationary solution of any Markov process on a finite configuration space can be obtained via the matrix tree theorem [89,90]. Specifically, one enumerates spanning in-trees on the directed graph of configurations where the edges are weighted by the rate at which a transition between the configurations takes place. In this way, one will always arrive at a configurational weight that is a sum of products of the elementary transition rates, and where there is a one-to-one correspondence between terms in the sum and spanning in-trees.
In principle, however, the space of trees is very large: since the configurational weight can be expressed as a determinant of an (n − 1)-dimensional matrix, where n is the number of configurations, we might expect the number of trees to be O(e N !) where N is the number of lattice sites. But in the case of the ASEP, the size of the extended configuration space (either lattice paths which grow exponentially with N or permutations which grow factorially with N ) is much smaller than the number of spanning trees implying a massive degeneracy in the weights of spanning trees. These observations show that it is difficult to predict a priori how large the extended space of configurations will be, and furthermore, that its size can be different in different parameter regimes.
which may be alternatively written as In (A.9) a and b are given, as in the main text, as (16) and we have used the q-binomial coefficient where the 'q-shifted factorial' is defined through In the case α = β = 1, a = b = −q and It is simple to check that the recursion and boundary conditions are satisfied by [7] n|C N |m = 2N and that for n = m = 0 we obtain for Z N the Catalan number C N +1 , using definition (32) We now show how Narayana numbers may be obtained in a similar fashion. Define G(N, P ) as the sum of all products of P D-matrices and (N − P ) E-matrices. Using the form of the matrix representation (30)(31), the following recursion relation holds [28] n|G(N, P )|m = n|G(N − 1, P − 1)|m + n|G(N − 1, P − 1)|m (B.5) Setting n = m = 0, we arrive at the Narayana numbers (48) which is the sum of weights of configurations with P particles. Summing over all configurations recovers (B.4): where we have used the Vandermonde identity [91] ∞ p=−∞ Appendix C. Demonstration of DE = D + E in path dominance problem Referring to Section 3.1, an ASEP configuration C is uniquely defined by the set of x or y-coordinates (x 0 , x 1 , . . . x P ), (y 0 , y 1 , . . . y Q ) that its path T traces. We now show that the weight of a path W(T ) has the same reduction relation (analogous to DE = D + E (6)), therefore the weight of the TASEP configuration W(C) = W(T ). If we split the path T into T = (T (1) , ↑, →, T (2) ), we are compelled to show relation (46) (illustrated in Figure 10), which we repeat here: This is analogous to the matrix relation DE = D + E. If T is defined by y coordinates (y 0 , y 1 , y 2 , . . . y i−1 , y i , y i+1 , . . . y Q−1 , y Q ), we therefore require W(y 0 , y 1 , y 2 , . . . y i−1 , y i , y i+1 , . . . y Q−1 , y Q ) (C.2) = + W(y 0 , y 1 , y 2 , . . . y i−1 , y i+1 , . . . y Q−1 , y Q ) (C.3) + W(y 0 , y 1 , y 2 , . . . y i−1 , y i − 1, y i+1 − 1, . . . y Q−1 − 1, y Q − 1) .
We first write W T (1) , →, T (2) explicitly, and rework this expression so to 'complete' each of the n i , n i+1 , . . . n Q−1 summations sequentially:  This nested expression then telescopes down to two sums, which can be identified as: We first note that for a Markov process in which for each transition out of a configuration we may associate a transition into the configuration occurring with the same rate, the stationary state has uniform probability. This is easy to see from the stationarity condition where P (C) is the stationary distribution and W (C → C ) is the transition rate from C → C . For P (C) to be independent of C, we simply require which is satisfied under the bijection of out-transitions onto in-transitions stated above.

D.2. Stationarity of simplified two-row TASEP dynamics
In order to show that the continuous-time Markov process illustrated in Figure 12 converges to a stationary distribution that is uniform in the space of complete configurations, we must demonstrate that: (i) each transition out of a complete configuration is into another complete configuration; (ii) every complete configuration can be accessed (by one or more transitions) by any other; and (iii) each transition out of a complete configuration can be mapped one-to-one onto a transition into that complete configuration.
Under these conditions the stationary distribution is defined over the full set of complete configurations (and only that set), and is uniform (see D.1).
Recall that the definition of a complete configuration is that the total number of particles on both rows of the lattice up to each site i = 1, 2, . . . N must be at least as large as the number of holes. This is, if n i is the number of particles upto site i, we must have n i ≥ 2i − n i , or equivalently, n i ≥ i. Note that the zone boundaries drawn in Figure 12 indicate the points where n i = 0. The two transitions that move particles between the two rows of the lattice leave all n i unchanged: consequently the final configuration is complete if the initial configuration is complete. If a particle on the bottom row hops to the left, one of the n i increases by 1, and the rest remain unchanged. Thus again, an initial configuration that is complete yields to a final configuration that is complete. Conversely, a particle on the top row hopping to the right causes one of the n i to decrease by 1, and the rest remain unchanged. Thus it could be possible to leave the space of complete configurations if n i = 0, i.e., if a particle crosses a zone boundary. Now, if a top-row particle can hop to the right at site i, site (i + 1) must be empty. The constraint n i+1 ≥ 0 implies that we must have a particle in the bottom row at site (i + 1), and in fact we must have n i+1 = 0 too. Meanwhile, the constraint n i−1 ≥ 0 means that we must have a hole in the bottom row at site i. This means that we can arrange for n i and n i+1 to both remain equal to zero if the bottom row particle at (i + 1) moves to the left at the same time as the top row particle moves to the right. Since none of the n i have changed, this configuration is also complete. This demonstrates statement (i) above.
To demonstrate statement (ii), we start by identifying two key configurations. The first has the top row empty, and consequently the bottom row full. Given any starting configuration, we can always reach this state through combinations of top-row hops, and movement of the top-right particle onto the bottom row. We can now reach any desired arrangement of particles on the top row by successive movement of the bottomleft particle onto the top row, and top-row hops. During these moves, we have n i = 0 for all i: this means that every top-row hop from site i to (i + 1) will be accompanied by a forced bottom-row hop from (i + 1) to i. Thus the complete configuration that is arrived at is the one where a hole sits below each top-row particle, and a particle below each top-row hole. There is no complete configuration that has any bottom row particles further to the right than this. Consequently any desired complete configuration can be reached by moving bottom row particles to the left; moreover, the top row remains fixed as this is done. Therefore, we have found a path from any complete configuration to any other, which implies that the Markov process is ergodic in the space of complete configurations.
To show that the number of ways into each configuration equals the number of ways out, we need to identify with each transition a 'driver' particle. In the cases where only one particle hops, that particle is the driver. These all sit at the front of a domain of particles (in the clockwise hopping direction). In the case where both a top and bottom row particle hop, the top-row particle is the driver (because the bottom row particle is forced to move to maintain positivity). Note that the bottom-row particle is the driver for a different transition. Consequently, the number of ways out of each configuration is equal to the number of domains. This situation with transitions into a configuration is more subtle. Each particle that is at the back of a domain may have moved in the last transition. For particles on the top row, the configuration that has that particle one site to the left (or on the bottom row, if it is the leftmost particle) is also complete, and is therefore one that could have been arrived from. On the bottom row, the configuration that has a particle one site to the right (on on the top row, if it is the rightmost particle) is complete, unless at the target site i, n i = 0. Movement of the bottom-row particle alone would imply having come from an incomplete configuration, which is not allowed. However, this configuration can be reached if a top-row particle in site (i + 1) was the driver in the previous move, and forced the bottom-row particle to the left, a move that has not yet been accounted for. Thus to each particle at the back of each domain, we can associate a unique preceding configuration; and consequently the number of ways into each configuration is also equal to the number of domains.
One can also show that these simplified dynamics satisfy a dynamic reversibility condition [40]. This involves an involution between complete configurations C → C (i.e., a self-inverse mapping of the set of complete configurations onto itself). Dynamic reversibility is said to hold for some set of transition rates W (C → C ) if the following conditions hold: where P (C) is the stationary distribution. If the 'hat' operation is the involution that exchanges the top and bottom rows, the first condition immediately follows from the fact that the number of ways out of each configuration is equal to the number of domains, as this is preserved under the involution. One can check the second condition exhaustively by considering each combination of particles and holes on adjacent sites, and noting that the stationary distribution is uniform.

Appendix E. Demonstration of reduction relations in permutation problem
Here we show that the permutation problem detailed in Sections 4.1, 5.1 has an equivalent reduction relation structure to those in the SSEP and PASEP. Following the formalism in these sections, define as shorthand for the total weight of permutations of the integers in (0, 1, . . . N ) where only the j = (j 1 , j 2 , . . . j P ) are raised. The reduction relation W |E = W | has an equivalent form which is trivial to show: for each permutation on the LHS of (E.2), increase every integer by 1, then append the permutation with a 0. This gives each permutation on the RHS, with all weights unchanged. Similarly for D|V = |V , which can be seen as each permutation on the LHS of (E.3), prepended with an (N + 1), corresponds to a permutation on the RHS. Again, all weights are unchanged. Finally, the reduction relation DE = qED + D + E has an equivalent form W N j 1 , k, j 2 (E.4) = qW N j 1 , k + 1, j 2 + W N −1 j 1 , k, j 2 − 1 + W N −1 j 1 , j 2 − 1 matrix, which allows its determinant, which we define detM N ×N ≡ Z N (F.1) to be expressed in a recursive form using (57), from Theorem 2.1 in [36] We now show that Z N and the TASEP partition function Z N (26) have the same generating function, thus making them equivalent. Define this generating function in η as Z