Stationary distributions of the multi-type ASEP

We give a recursive construction of the stationary distribution of multi-type asymmetric simple exclusion processes on a finite ring or on the infinite line $Z$. The construction can be interpreted in terms of"multi-line diagrams"or systems of queues in tandem. Let $q$ be the asymmetry parameter of the system. The queueing construction generalises the one previously known for the totally asymmetric ($q=0$) case, by introducing queues in which each potential service is unused with probability $q^k$ when the queue-length is $k$. The analysis is based on the matrix product representation of Prolhac, Evans and Mallick. Consequences of the construction include: a simple method for sampling exactly from the stationary distribution for the system on a ring; results on common denominators of the stationary probabilities, expressed as rational functions of $q$ with non-negative integer coefficients; and probabilistic descriptions of"convoy formation"phenomena in large systems.


Introduction and main results
The asymmetric simple exclusion process (ASEP) is one of the simplest examples of an interacting particle system. In physcial terms it can be seen as a basic model of non-reversible flow. Despite its simplicity, it displays an extremely rich behaviour, and attracts intense study in physics, probability, combinatorics, and beyond.
In this article we consider ASEPs with multiple types of particle. We start by defining the model on the ring Z L = {0, 1, 2, . . . , L − 1} (with cyclic boundary conditions). The system is a continuous-time Markov chain with state-space {1, 2, . . . , N, ∞} Z L , for some N ≥ 1. For a configuration η = (η i , i ∈ Z L ), we say that η i is the type of the particle at site i. (We may sometimes refer to particles of type ∞ as holes -we could equally give them label N + 1 rather than ∞, but the different notation is sometimes helpful.) Fix some q ≥ 0. The dynamics of the process are as follows. If η(i) > η(i + 1) (where the addition i + 1 is done mod L), then the values η(i) and η(i + 1) are exchanged at rate 1. If instead η(i) < η(i + 1), then the values η(i) and η(i + 1) are exchanged at rate q. That is, neighbouring particles rearrange themselves into increasing order at rate 1, and into decreasing order at rate q. We will assume q ∈ [0, 1) throughout. Such a process was perhaps first described by Alcaraz and Rittenberg [1] as an example of a more general class of multi-type particle systems.
The dynamics preserve the number of particles of each type. Restricted to a state space with a given number of particles k n of each type, for n = 1, . . . , N and k 1 + k 2 + · · · + k N ≤ L, the process is irreducible. Our particular focus is on the stationary distribution.
Consider first the case N = 1. This is the standard (one-type) ASEP (for an extensive introduction, see Part III of the book of Liggett [34]). For the one-type ASEP on the ring, each site of Z L contains either a particle (type 1) or a hole (type ∞). A particle exchanges places with a hole to its left at rate 1, and with a hole to its right at rate q. Fix the number of particles to be k. Then it's well known (and straightforward to show) that the stationary distribution of the process is uniform on all L k configurations.
The ASEP with types {1, . . . , N, ∞} can be viewed as a coupling of N one-type ASEPs. Specifically, for any n = 1, . . . , N , we can consider a projection under which types r ≤ n are considered "particles" and types r > n are considered "holes"; for each n, the image of the process under this projection is a one-type ASEP. In particular, projecting the stationary distribution in this way gives a uniform distribution on the corresponding state-space.
However, despite the fact that all these projections are uniform, it is far from the case that the full stationary distribution is uniform. See for example Figure  1.1 for samples from the stationary distribution of a system with 1000 particles all of different types, on the ring with 1000 sites, for different values of q. A striking feature of the configurations observed is the appearance of long strings of particles with similar labels ("convoys"). This "clustering" is more pronounced for small q, and more pronounced around the middle of the range of particle types.
Before stating the main results of the paper, we mention certain previous results -see Section 1.4 for further background. The stationary distribution of the 2-type ASEP (where the non-uniformity mentioned in the previous paragraph already appears) was studied by Derrida, Janowksy, Lebowitz and Speer [19] (see also [18]) using the matrix product ansatz. Further combinatorial descriptions in the case of the TASEP (totally asymmetric, i.e. q = 0) were given by Ferrari, Fontes and Kohayakawa [23], Angel [3], and Duchi and Schaeffer [20]. In [26], Ferrari and Martin introduced a new recursive method which gives a construction of the stationary distribution of the N -type TASEP for any N . Samples of the stationary distribution of an ASEP on a ring of 1000 sites, with 1000 different types of particle. The horizontal axis is labelled by site, and the vertical axis by particle type. Top: q = 0; middle: q = 0.1; bottom: q = 0.8. The "convoys" are more pronounced when q is small (and the effect is stronger around the middle of the range of particle classes). See the results in Section 6.
The construction by Ferrari and Martin is written in terms of a system of queues in tandem (or so-called "multiline diagrams"), and the proof uses time-reversal arguments. A different approach, but exploiting a related recursive structure, was then found by Evans, Ferrari and Mallick [21], to give a representation of the stationary distribution of the N -type TASEP using the matrix product ansatz (this was elaborated further for example by Arita, Ayyer, Mallick and Prolhac in [6]).
This matrix-product representation was then extended to q > 0 by Prolhac, Evans and Mallick in [37]. They raised the question of whether the more probabilistic/combinatorial construction of [26] could also be generalised to the ASEP. In this paper we show that indeed it can, giving a multi-line queue construction of the stationary distribution of the ASEP.
This construction has a number of nice consequences: • algorithmic: it gives an efficient method for sampling exactly from systems even with a large number of types (see Algorithm 2 below, which was used to generate the samples in Figure 1.1 with 1000 types of particle); • algebraic: it leads directly to an expression for a common denominator of the stationary probabilities of an ASEP system, expressed as rational functions of q (see Theorem 1.3); • probabilistic: it can be exploited to give interesting qualitative and quantitative information about the properties of the stationary distribution, including the clustering properties described in Theorem 1.5.
A central object in the construction for the TASEP in [26] is a simple model of a priority queue in discrete time, with a Markovian service process. When a service occurs in the queue, the customer in the queue with highest priority departs (unlesss the queue is empty, in which case the service is unused). To generalise to q > 0, we introduce a queue with "rejected services". When a service occurs, it is offered to each customer in turn, in order of priority. Each customer accepts the service with probability 1 − q; the first customer to accept it departs (unless they all reject it, in which case the service is unused). The multi-line process can be seen as a collection of several such queues in series; such a model is closely related to the so-called q-TASEP (especially the discrete-time version considered for example in [11]).
It would be possible to give a proof of our main result by an elaboration of the methods involving dynamic reversibility which were used in [26] [27]. Instead, in this paper we rely directly on the matrix-product representation of Prolhac, Evans and Mallick, which we explain in Section 2. We make use of particular instances of matrices which satisfy the quadratic algebra of [37]; these can be related to Markov transition matrices which govern the evolution of a queue.

Algorithms and multi-line diagrams
We now begin to state the main results of the paper. We start with an algorithm for sampling from the stationary distribution of the N -type ASEP on Z L , with given particle counts. Later in the paper we will be able to restate it in terms of a function which assigns weights to "multi-line diagrams", such that the probability distribution on the diagrams proportional to these weights projects to the ASEP stationary distribution.
We describe the algorithm using queueing language, although we postpone until later the introduction of notation which describes queue-length processes more formally (see Section 3). The sites of Z L are treated as "times" in the queueing process. Hence time is cyclic (a setting in which time is more straightforwardly indexed by Z is introduced in Section 3 where we consider ASEPs on Z rather than Z L ).
First consider the case N = 2. Fix L, k 1 , k 2 with k 1 + k 2 ≤ L. The following algorithm can be used to generate a sample from the stationary distribution of the 2-type ASEP on Z L with k 1 type-1 particles and k 2 type-2 particles. The algorithm does the following thing: (i) it randomly chooses a set of k 1 arrival times, and a set of k 1 + k 2 service times; (ii) it assigns to each arrival time a different departure time chosen from among the service times. This gives an output of the queue, in which each time i ∈ Z L is either a departure time, an unused service time, or a time which had no service. This is used to define a configuration in {1, 2, ∞} Z L . Algorithm 1.
• Choose a set of "arrival times" uniformly among all subsets of Z L of size k 1 , and independently a set of "service times" uniformly among all subsets of Z L of size k 1 + k 2 .
• Now process the arrivals one by one, in an arbitrary order (for example left to right), and assign a departure time to each one from among the service times, in the following way.
Suppose we have already processed r of the k 1 arrivals, where 0 ≤ r < k 1 . This means we have already assigned r of the k 1 + k 2 service times. Now look at an arrival we have not processed yet -say it occurs at time i. We wish to assign a service time to be the departure time of the arrival at time i, which has not yet been assigned to any other arrival. If there is a service at time i which has not yet been assigned, then assign it to this arrival, and we are done. Otherwise, let the remaining k 1 + k 2 − r available potential service times be i 1 , i 2 , . . . , i k 1 +k 2 −r ; we list these sites in cyclic order around the ring starting from i, so that Now assign the arrival at i to the service at i j with probability q j−1 /(1 + q + q 2 + · · · + q k 1 +k 2 −r−1 ). (1.1) • Having done this for all k 1 arrivals, we have chosen k 1 departure times, and the remaining k 2 − k 1 service times will be unused. Now define an output configuration D ∈ {1, 2, ∞} Z L as follows: for each i ∈ Z L , Theorem 1.1. Algorithm 1 generates a configuration D i , i ∈ Z L whose distribution is the stationary distribution of the ASEP on Z L with k 1 type-1 particles and k 2 type-2 particles.
Now we generalise this algorithm to all N ≥ 2, in a recursive way. Fix L and k 1 , . . . , k N . The following algorithm generates an algorithm from the stationary distribution of the N -type ASEP on Z L with particle counts k 1 , . . . , k N (that is, with k n particles of type n for 1 ≤ n ≤ N ). The input into the algorithm is a collection of k 1 + · · · + k N service times, along with an arrival process A = (A i , i ∈ Z L ) ∈ {1, 2, . . . , N, ∞} Z L , which is itself drawn from the stationary distribution of the (N − 1)-type ASEP on Z L with particle counts k 1 , . . . , k N −1 . If A i = n where 1 ≤ n ≤ N , there is an arrival of type n at time i. If A i = ∞, there is no arrival at time i. The algorithm gives an output in which each time in Z L is a departure type of some customer with a type n ∈ {1, 2, . . . , N − 1}, or an unused service time, or a time which had no service. This is used to define a configuration in {1, 2, . . . , N } Z L .

Algorithm 2.
• Choose the arrival process according to the stationary distribution of the (N − 1)-type ASEP on Z L with particle counts k 1 , k 2 , . . . , k N −1 , and independently choose the set of service times uniformly among all subsets of size K = k 1 + k 2 + · · · + k N .
• Now consider the arrival times one by one, as in Algorithm 1, but importantly this is now done in order of type. We start by considering each of the type-1 arrivals, then each of the type-2 arrivals, and so on. The order in which arrivals of the same type are considered is arbitrary (for example, left to right).
We assign a different service time to each considered arrival in turn, as follows. Suppose we have already processed r of the arrivals, where 0 ≤ r < k 1 + k 2 + · · · + k N −1 . This means we have already assigned r of the K service times. Now look at an arrival we have not processed yet -say it occurs at time i. We wish to assign it a departure time from among those service times which are not yet assigned. If there is a service at time i which has not yet been assigned, then assign it to this arrival, and we are done. Otherwise, let the remaining K −r available service times be i 1 , i 2 , . . . , i K−r ; we list these sites in cyclic order around the ring starting from i, so that Now we assign the arrival at i to the service at i j with probability • Having done this for all arrivals, we have chosen k 1 + k 2 + · · · + k N −1 departure times, and the remaining k N service times are unused. Now define the configuration D ∈ {1, 2, . . . , N } Z L as follows: Theorem 1.2. Algorithm 2 generates a configuration (D i , i ∈ Z L ) whose distribution is the stationary distribution of the N -type ASEP on Z L with particle counts k 1 , . . . , k N .
To generate a sample from an N -type TASEP, we need to apply Algorithm 2 N − 1 times in all. The nth iteration takes as arrival process a configuration whose distribution is stationary for an n-type ASEP, along with an independent service process, and outputs a configuration whose distribution is stationary for an (n + 1)type ASEP. (The first iteration, with n = 1, is equivalent to Algorithm 1.) We can combine the N − 1 iterations into a "multi-line diagram" with N lines, as was done for the q = 0 case in [26]. The nth line of the diagram is a n-type ASEP configuration, i.e. a configuration in {1, . . . , n, ∞} Z L , with particle counts k 1 , . . . , k n . It is the arrival process for the nth iteration of the algorithm (if 1 ≤ n ≤ N − 1) and the output configuration of the (n − 1)st iteration of the algorithm (if 2 ≤ n ≤ N ). Specifically, the last line of the diagram gives a sample from the N -type stationary distribution.
On the nth line, the number of occupied sites (i.e. the sites where the value of the configuration is finite rather than infinity) is k 1 + · · · + k n . Ignoring the types, the set . The configuration of occupied sites is the same in both cases. The upper diagram then shows the deterministic assignment of types to particles in the case q = 0. The lower diagram shows one of the (many) other possible assignments of types to particles when q > 0, and "rejected services" are allowed.
of occupied sites is uniform among all subsets of Z L of size k 1 + · · · + k n ; furthermore, the sets of occupied sets on different lines is independent. Note that for q = 0 (the case of the TASEP), the only randomness in Algorithm 2 is in the first step when the arrival and service processes are chosen. In that case, whenever the algorithm looks to assign a departure time to an arrival at time i, it chooses the first so-far unassigned service time in i, i + 1, i + 2, . . . . In this case, given the sets of occupied sites on the different lines of the multi-line diagram, the assignment of types to the occupied sites is deterministic. This gives the original algorithm of [26]. Figure 1.2 shows two multi-line diagrams for the case N = 4, L = 10, and (k 1 , k 2 , k 3 , k 4 ) = (1, 3, 3, 1). The sets of occupied sites are the same for the two diagrams. The upper diagram is the one which results (deterministically) from the configuration of occupied sites in the case q = 0. The lower diagram shows another possible configuration (one of many) when q > 0.

Common denominators
By considering the various sources of randomness that go into the construction of a multi-line diagram using repeated applications of Algorithm 2, we can obtain the following result giving a common denominator for the ASEP probabilities. We write [n] q = 1 + q + · · · + q n−1 , and [n] q ! = [1] Theorem 1.3. Consider the N -type ASEP on Z L , with particle counts (k 1 , . . . , k N ). The stationary probability of every state is a polynomial in q with non-negative integer coefficients divided by the common denominator A special case is that of a process on Z L with L different classes of particle, as depicted in Figure 1.1. All other processes on Z L can be obtained as projections of this process. We can take N = L−1 and k 1 = k 2 = · · · = k L−1 = 1 (it is equivalent to regard the final particle as a particle of type L or as a hole). This gives the following result: Corollary 1.4. The stationary probabilities of the ASEP on Z L with one particle of each of L different classes are given by polynomials in q with non-negative integer coordinates divided by the common denominator The expression in (1.6) can also be written as We note that a related formula appears in work of Cantini, de Gier and Wheeler [14] in a more general setting of Macdonald polynomials; specialising their formula in their Section 5 to the ASEP, one obtains an expression like (1.7) but with an extra factor of (1 − q) L(L−1)/2 , and without the restriction to non-negative integer coefficients in the numerator. The denominators given by (1.5) and (1.7) may well not be the best possible. In fact, these are common denominators for the probabilities of all multi-line diagrams; the projection from the multi-line diagrams to the bottom line giving a single ASEP configuration may introduce further common factors. 1 + 7q + 7q 2 + 9q 3 96(1 + q)(1 + q + q 2 ) Figure 1.3: Stationary probabilities for a system with 4 sites (numbered clockwise around the ring). The particles are labelled with their type; the single hole may equivalently be seen as a particle of type 4. There are 6 possible configurations, up to rotation. Note the symmetry between q and 1/q; the probabilities remain unchanged under replacing q by 1/q and reversing the order of the particles (since replacing q by 1/q is equivalent to exchanging left and right and multiplying time by a factor 1/q; the time-change has no effect on the stationary distribution). Note also that all expressions are equal when q = 1, since this gives the symmetric exclusion process whose equilibrium is uniform over all configurations.
As an example, we can consider the case of the 4-type system on a ring of size 4. Corollary 1.4 gives a common denominator of 96(1 + q) 2 (1 + q + q 2 ) = 96(1 + 3q + 4q 2 + 3q 3 + q 4 ). However, in fact we can dispense with one of the factors of 1 + q; the stationary probabilities as a function of q are shown in Figure 1 The extra symmetries of the multiline diagram may be easier to explore in the context of the alternative queueing construction mentioned in Section 7, where we do not need to treat particles differently according to whether a service is available to them at the time they enter the queue.

Clustering
Our final result explains the "clustering" phenomenon visible in Figure 1.1, showing samples from the stationary distribution of the ASEP with 1000 distinct labels 1, 2, . . . , 1000 on the ring with 1000 sites, for different values of q. A striking feature of the configurations is the appearance of long strings of nearby particles with similar labels ("convoys"). In this section we explain some aspects of that picture, by considering an appropriate local limit of stationary distributions as L → ∞. The analysis will rely in an integral way on the queueing construction of the stationary distribution of the 2-type system on Z which we develop in Section 3 (see Theorem 3.2).
Consider the ASEP on Z L , with sites now written for convenience as − L 2 , . . . , −1, 0, 1, . . . , L 2 − 1 in a cyclic way, and with L particles with distinct types 1, 2, . . . , L. (We could equally write the largest-numbered type as ∞, i.e. a hole.) The process is irreducible, and its stationary distribution can be constructed via an "(L − 1)-line process" involving repeated applications of Algorithm 2, as described in Section 1.1.
Let Y (L) i be the type of the particle at site i, in a configuration distributed according to this stationary distribution. We will rescale linearly, so that all labels lie in [0, 1], and also pad the configuration with infinite strings of zeros on each side, to define a configuration W (L) = (W The following results describe some aspects of the clustering phenomenon observed in the configurations in Figure 1.1: i , i ∈ Z) converges in distribution (with respect to the product topology on [0, 1] Z ) to a limit which we denote by W = (W i , i ∈ Z). The distribution of W is translation-invariant.
(b) (W 0 , W 1 ) has joint distribution on [0, 1] 2 with density given by (1.10) [By this we mean that if A is a measurable subset of [0, 1] 2 , then and if B is any measurable subset of [0, 1], then It follows that (c) With probability 1, there are infinitely many k such that W 0 = W k . is uniformly distributed on {1, 2, . . . , L}). However, it is definitely not the case that the W i are independent; indeed, two neighbouring values have positive probability to be equal! Note that the density f of two neighbouring values (W 0 , W 1 ) on x > y decays to zero as the diagonal is approached if q = 0, but not if q > 0.
Going beyond part (c), the methods of Section 6 in fact explain how to describe the distribution of the "convoy" {k ≥ 1 : W k = W 0 } conditional on {W 0 = x} for x ∈ [0, 1] somewhat explicitly in terms of a random walk on Z whose transition probabilities depend on x.
In the TASEP case (q = 0), the process W of Theorem 1.5 was considered in a related context by Amir, Angel and Valkó [2]. Consider a TASEP on Z started from an initial state in which site i contains a particle of type −i, for each i ∈ Z; the dynamics of the process are that each pair of neighbouring particles sort themselves into increasing order at rate 1. Let X i (t) be the position at time t of the particle with label i. The limit U i = lim t→∞ with probability 1 for all i (as follows from a result of [36]), and has Uniform[−1, 1] distribution. The TASEP speed process (U i , i ∈ Z) was introduced and studied by [2]. Rescaling by W i = (1 − U i )/2, we get the process (W i , i ∈ Z) with the distribution described in Theorem 1.5. The properties (1.9)-(1.11) for the case q = 0 are given in Theorem 1.7 of [2], and further results concern for example the joint distribution of more than two entries.
Still in the case q = 0, particularly sharp results for stationary distributions on finite rings were then given by Ayyer and Linusson in [9]. Their results include closed-form expressions for "three-point" probabilities of the form P(Y = n), as well as more general "two-point" probabilities of the form P(Y (By taking appropriate limits, one can regain some of the formulas from [2] for the speed process.) The proofs involve an intricate combinatorial analysis of the multi-line queue construction for the TASEP. It would of course be interesting to explore whether an analogous application of the multi-line constructions presented here could lead to similar results for multi-point probabilities in the case q > 0.

Related work
Let us mention some further related work, in addition to that discussed above. The recursive approach to the construction of multi-type particle systems has been extended to a variety of different particle systems, including discrete-time TASEPs [35], inhomogeneous (or "multi-rate") versions of the multi-type TASEP [7,8,13], and a variety of zero-range processes [32,33,31], and also to the description of joint distributions of Busemann functions for last-passage percolation [22].
Connections between the multi-type ASEP on the ring and families of symmetric polynomials such as Schubert polynomials and Macdonald polynomials have been studied by various authors, for example by Cantini in [13] and in recent work by Corteel, Mandelshtam and Williams in [16]; the latter paper involves a description of Macdonald polynomials using objects related to the multi-line diagrams described in this paper.
Stationary distributions for the multi-type ASEP on a finite interval with open boundary conditions have also been widely studied. Not all sets of boundary rates are expected to lead to exactly solvable models; however, a variety of classes of integrable boundary conditions have recently been established -see for example work of Crampe, Finn, Ragoucy and Vanicat [17,28], and also work of Cantini, Garbali, de Gier and Wheeler [14,15] where further connections to families of orthogonal polynomials are made.
Systems with closed boundaries on finite or half-infinite intervals have recently been considered by Belitsky and Schütz [10] (see also [5]); as well as describing sta-tionary distributions, they obtain duality properties and use them to study hydrodynamic limits (including the behaviour of shocks). The stationary distributions of these systems are closely related to certain stationary distributions for the system on Z which (unlike those considered in this paper) are not translation-invariant; instead their projections onto one-type distributions are blocking measures in the sense of [12]. The connections between distributions of this type and the Mallows measure on permutations were previously studied by Gnedin and Olshanskii [29,30], and recent work by Angel, Holroyd, Hutchcroft and Levy [4] describes a link between such processes and a model of stable matchings.
Integral formulas for multi-type ASEPs on Z with finitely many particles are given by Tracy and Widom in [38].

Plan of the paper
The rest of the paper is organised as follows.
In Section 2 we introduce the recursive matrix product construction of Prolhac, Evans and Mallick [37]. We give specific realisations of the relevant matrices, which are closely related to Markov transition matrices for the queueing systems that we use to construct the multi-type ASEP stationary distribution.
These queueing models are introduced in Section 3. We also give results concerning stationary distributions of multi-type ASEPs on Z, which may also be considered as main results of the paper in their own right.
The proofs of Theorems 1.1, 1.2 and 1.3 are then developed in Section 4. We cover the particular case N = 2 in some detail, with the aim of making the argument in the more general case as easily comprehensible as possible.
The results concerning stationary distributions on Z are proved in Section 5; these are deduced from the corresponding results for processes on the ring using a rather intricate coupling argument, which may be of independent interest.
Those results form a central part of the proof of Theorem 1.5, given in Section 6. The limit process W of Theorem 1.5 can be identified as a stationary distribution for an ASEP on Z whose particle-types are continuous, and distributed uniformly on [0, 1]. This process is a generalisation of the "TASEP speed process" studied by Amir, Angel and Valko in [2] (see also [24] for closely related results). The process W can be studied via the projection of the continuous-type ASEP onto an N -type process; much useful information can be extracted already from the case N = 2. The hardest part of the proof is the argument to establish Theorem 6(c), where a subtle argument involving stochastic domination between random walks is required.
Finally in Section 7 we briefly discuss an alternative construction involving a modified queueing discipline; this would result in a more complicated matrix-product structure, but a rather simpler and more natural formulation in terms of multi-line diagrams and their weights.

Matrix product framework
In this section we describe the matrix product representation for the stationary distribution of the process on the ring given by Prolhac, Evans and Mallick [37], which is the starting-point of the proof of the multi-line construction for the ASEP stationary distribution. Our notation will be slightly different from in that paper; since we consider the ASEP with jumps left rather than right, the matrices which appear in the matrix product solution are the transposes of those used in [37]. This change is convenient since now time in the associated queueing systems flows from left to right, and the matrices which appear are closely related to Markov transition matrices.
Suppose the matrices δ, and α satisfy the following relations: At (2.5) below, we'll give specific examples of appropriate matrices (which will be infinite-dimensional) α, δ, . Now we define matrices X ∞ be scalars and equal to 1.
where the matrices a (N ) m,j are given by The dynamics of the ASEP on the ring Z L preserve the number of particles of each type. Consider the process with k n particles of type n, for n = 1, 2, . . . , N , where k n > 0 for all n and also k 1 + k 2 + · · · + k N < L.
Theorem 2.1 (Prolhac, Evans and Mallick [37]). The stationary distribution of the N -type TASEP on Z L , with k n particles of type n for n = 1, 2, . . . , N , is given by where Z (N,L) k 1 ,...,k N is a normalizing constant chosen such that the sum of the right-hand side over all configurations with the particle counts k 1 , . . . , k N is 1.
Matrices satisfying the relations (2.1) can be realised in many ways. We give a particular version with direct links to the queueing interpretations that we will develop: Note immediately that is a stochastic matrix (that is, a matrix whose entries are non-negative and whose row-sums are all equal to 1), and also δ + α is also a stochastic matrix. Of course, so is the identity matrix I. The rows and columns of all these matrices are considered to be indexed by Z + = {0, 1, 2, . . . }. The indices can be seen as queue lengths, and the matrices are interpreted as transition matrices in a queueing process evolving over time.
Each matrix a (N ) m,j is a tensor product of N − 1 of the fundamental matrices α, δ, , I, and can be seen as indexed by vectors in Z N −1 + , representing queue-lengths of customers of types 1, 2, . . . , N − 1 in an (N − 1)-type queueing process. The quantity a (N ) m,n will represent the weight of a transition associated with the arrival of a customer of type m and departure of a customer of type n (with ∞ representing no customer).
In fact, the non-zero contributions to X (N ) n are all terms of the form where m r ∈ {1, . . . , r, ∞} for r = 1, . . . , N − 1. This quantity represents the weight of a transition in a system of N queues in series, associated to the arrival of a customer of type m 1 in the first queue, a transfer of a customer of type m r from queue r − 1 to queue r for 2 ≤ r ≤ N − 1, and a departure of a customer of type n from queue N − 1 (and the value ∞ indicates the absence of a customer).

The multi-type ASEP on Z
The multi-type ASEP on the whole line Z is defined analogously to the process on the ring Z L . The N -type system is a continuous-time Markov process with state-space {1, 2, . . . , N, ∞} Z . For a configuration η = (η i , i ∈ Z), say that η i is the type of the particle at site i. The dynamics are as follows: if η(i) > η(i + 1), then the values η(i) and η(i + 1) are exchanged at rate 1. If instead η(i) < η(i + 1), then the values η(i) and η(i + 1) are exchanged at rate q.
In the case N = 1, the ergodic translation-invariant stationary distributions are all Bernoulli product measures, in which each site contains a particle (type 1) with probability λ, and otherwise a hole (where λ ∈ [0, 1]). (There are also non-translationinvariant stationary distributions, the blocking measures considered for example by [12]; these are concentrated on configurations with only finitely many holes to the left of the origin and only finitely many particles to the right of the origin.) As shown in [25], for given λ 1 , . . . , λ N with λ 1 + · · · + λ N < 1, there is a unique ergodic translation-invariant stationary distribution for the N -type ASEP on Z in which the intensity of type-n particles is λ n for 1 ≤ n ≤ N . We denote it by ν (N ) λ 1 ,...,λ N . For each n = 1, . . . , N , we can consider a projection under which types r ≤ n are considered "particles" and types r > n are considered "holes". Under this projection, the N -type ASEP becomes a one-type ASEP, and so in particular the image of ν (N ) λ 1 ,...,λ N is Bernoulli product measure with intensity λ 1 + · · · + λ r . However, although all these projections are product measures, ν (N ) λ 1 ,...,λ N is not itself a product measure! In the case q = 0, this stationary measure on Z was constructed in [26].

Single-type queue
We now define the basic model of a discrete-time queue including rejected services, which will be used to describe the stationary distribution of the two-type ASEP on Z. In later sections we consider systems consisting of several such queues in series, with multiple types of customer -these will be used to describe multi-type stationary distributions. We then describe analogous systems where the "time" index is cyclic, in order to describe the stationary distributions of systems on the ring Z L .
The queue is Markovian. Write A i = 1 if there is an arrival at time-step i, and A i = ∞ otherwise. Write S i = 1 if there is a service available at time-step i, and S i = ∞ otherwise. The processes of arrivals and of services are independent Bernoulli processes, with rates λ and µ respectively. That is, at each time-step, an arrival occurs with probability λ, and then independently a service is available with probability µ (with independence between different time-steps).
Suppose the queue-length at the start of the time-step is k. There are four possibilities: • No arrival occurs, no service available, with probability (1 − λ)(1 − µ). The queue-length remains k.
• Arrival occurs and service available, with probability λµ. A departure occurs, and the queue-length remains k.
• No arrival occurs, service is available, with probability (1 − λ)µ. With probability 1 − q k , a departure occurs, and the queue-length goes down to k − 1. With probability q k , an unused service occurs, and the queue-length remains k.
Note that an unused service is allowed to occur only if no arrival has occurred at that time-slot. We can imagine the service mechanics as follows. When a service is available, it is offered to each customer in turn. Each one in turn accepts it with probability 1−q and rejects it with probability q, except that a customer who has just arrived must always accept. As soon as a customer accepts the service, that customer departs and we stop. If all k of the customers reject the service (with probability q k if no arrival has occurred) then it remains unused.
The transition matrix of the queue-length process is given by where α, δ, are given in (2.5). The four terms in (3.1) correspond to the four possibilities for the evolution of a time-step listed above. The term in δ corresponds to transitions in which a departure but no arrival occurs. α corresponds to transitions in which no arrival occurs and a service is offered but unused. corresponds to transitions in which arrival occurs but no service is offered. Finally, the terms with I correspond to transitions where either there is arrival and service, or there is no arrival and no service is offered.
For stability we assume λ < µ. In that case, the queue-length process is positive recurrent, and there is a unique equilibrium version.
Note that at each time i ∈ Z, one of three possibilities occurs: a departure, an unused service, or no available service. As at (1.2), we define a departure process Let Q i be the number of customers in the queue at the beginning of time-step i. Then we have the following recursive formula: λ,µ−λ (the unique ergodic stationary distribution of the two-type ASEP on Z with parameter q and with density λ of first-class particles and µ − λ of second-class particles).
Our proof of this result will be based on the matrix product representation in Theorem 2.1.

Multi-type queues
Next, we consider the extension of the queue described above to a multi-class queue with priorities. The queue will contain N − 1 classes (or types) of customer, labelled 1, 2, . . . , N − 1. The lower the number of the class, the higher the priority.
The state of the queue at a given time i is now a vector (Q denotes the number of customers of class n present in the queue at the beginning of time-step i. At each time-step i, at most one customer arrives (with some given class). Write A i = n if a customer of type n arrives, and A i = ∞ if there is no arrival. Initially we don't specify the arrival process, but concentrate on describing the action of the queueing server.
As before, the process of available services is a Bernoulli process of rate µ; write S i = 1 if a service occurs (which happens with probability µ) and S i = ∞ otherwise. As above, the service will be offered in turn to each customer in the queue. This is now done in order of priority; the service is offered to each of the first-class customers, then to each of the second-class customers, and so on until some accepts it or all customers have rejected it. Each customer accepts the service with probability 1 − q and rejects it with probability q, with the exception that if a customer has just arrived at the queue, then that customer will always accept.
We give two brief examples. Suppose that at the beginning of a time-slot, the queue contains 3 first-class customers, 1 second-class customer and 4 third-class customers. Suppose that there is no arrival, and a service is available. Then a departure of a first-class customer occurs with probability 1 − q 3 , a second-class departure occurs with probability q 3 − q 4 , a third-class departure occurs with probability q 4 − q 8 and an unused service occurs with probability q 8 .
Suppose instead that an arrival of a second-class customer occurs, increasing the number of second-class customers to 2. Now a first-class customer departs with probability 1 − q 3 , and a second-class customer departs with probability q 3 . It is impossible for a third-class departure or an unused service to occur.
Note the particular case q = 0. Here it is always a customer of the highest priority-type present who departs, and a service is unused if and only if the queue is empty. As

Multi-type stationary distribution
We now explain how to construct stationary distributions for the multitype ASEP on Z recursively. The departure process of a queue with N − 1 types, whose arrival process corresponds to the stationary distribution for the (N − 1)-type ASEP, is used to give a stationary distribution for the N -type ASEP.
Hence the stationary distribution with N types can be seen as the output of a series of N − 1 queues in tandem. The rth of the N − 1 queues contains r types of customer in its arrival process. Its departure process also contains r types, to which we add an (r + 1)st type corresponding to unused services.

"Queues" on the ring: two-type
To construct the stationary distribution of the ASEP on Z L = {0, 1, . . . , L − 1} (the ring with L sites) we consider "queues" in which the time is cyclic; i.e. we replace the time index Z by Z L . All addition and subtraction is to be understood modulo L In this section we cover the two-type ASEP. The numbers of first-class and of second-class particles are conserved by the dynamics. For each k 1 and k 2 with k 1 + k 2 < L, the set of configurations with k 1 first-class and k 2 second-class particles (and hence L − k 1 − k 2 holes) forms a single communicating class; there is a unique stationary distribution concentrated on such a set of configurations.
The construction of the two-type stationary distribution uses a queue with one type of customer. We will have k 1 arrivals and k 1 + k 2 services.
The set of times at which arrivals occur is chosen uniformly from all subsets of Z L of size k 1 , and the set of times at which services are available is chosen uniformly from all subsets of Z L of size k 1 + k 2 , independently of arrivals.
Analogously to the queue described in Section 3.2, we want the following rules. If no service is available, then no departure occurs. If an arrival and an available service both occur, then a departure must occur. Finally, if a service is available but an arrival does not occur, there may be either a departure or an unused service.
Formally, suppose we are given arrival and service processes A and S in {1, ∞} Z L , Then we say that a queue-length process If such a D exists it is unique. (The first line determines when D i = 1 and the next two lines together determine when D i = ∞.) For (A, S, Q) ∈ R (2) , we write D(A, S, Q) for the unique D satsifying the properties in (3.5). Now, given A and S, define a weight function on valid queue-length processes as follows. For valid Q, define the weight w i (Q|A, S) associated to site i ∈ Z L by where D = D(A, S, Q). (It's straightforward to check that if the properties in (3.5) are satisfied, then exactly one of the cases in (3.6) occurs.) Now define the weight of the whole process Q by w(Q|A, S) = i∈Z L w i (Q|A, S). Given A and S, we now take the probability of the queue-length process Q to be proportional to w(Q|A, S), by It's easily seen that the denominator is finite as long as q < 1. In fact, we will show later (Lemma 4.2) that this normalizing constant depends only on k 1 and k 2 , and not on particular A and S. The weight has the following interpretation. At each time slot where a service occurs, this service is offered to the customers in turn. Each customer accepts an offer with probability 1 − q and rejects it with probability q, with the exception that any customer who has just arrived at the queue must accept it. The values of w i in the second, third, and fourth lines of (3.6) above are the corresponding probabilities for observing a departure or an unused service, taking into account the current composition of the queue and the information about whether an arrival has occurred. Now we use the weight w to give a distribution on departure configurations D; namely, choose A and S uniformly, as described above; choose Q in proportion to P (Q|A, S) (equivalently, in proportion to w(Q|A, S)); finally, take D = D(A, S, Q). Theorem 3.3. The distribution of the configuration D = D i , i ∈ Z L induced by the weight (3.6) on queue-length processes is the stationary distribution of the two-type ASEP on Z L with k 1 first-class and k 2 second-class particles.

"Queues" on the ring: multi-type
In this section we construct the stationary distribution for the ASEP on Z L with several types of particle, generalising the result of the previous section for 2-type systems.
The construction is done recursively. To describe the N -type equilibrium, we use the (N − 1)-type equilibrium and a priority "queue" (with cyclic time as in the previous section).
The arrival process contains N − 1 types of particle, and holes. Let k n be the number of particles of type n. Write A i = n if there is an arrival of type n at time i, and A i = ∞ if there is no arrival. Write S i = 1 if there is a service available at time i and otherwise S i = ∞.
We choose the process A according to the stationary distribution of an (N − 1)type system with k n particles of type n for 1 ≤ n ≤ N − 1. Independently of the arrivals, the times of potential services are chosen uniformly from all subsets of Z L of size k 1 + k 2 + · · · + k N .
We now consider queue-length processes Q We want the queue to obey the following rules. If no service occurs, then no departure can occur. If an arrival and a service both occur, then a departure must also occur, and the departing customer must have type no larger than that of the arrival. If a service occurs but an arrival does not occur, there may be either a departure or an unused service.
Formally, suppose we are given arrival and service processes A ∈ {1, 2, . . . , N − 1, ∞} Z L and S in {1, ∞} Z L , Then we say that a queue-length process Q = (Q  Then for valid Q, define the weight w i (Q|A, S) associated to site i ∈ Z L by Given A and S, the probability P (Q|A, S) of the queue-length process Q is now taken to be proportional to w(Q|A, S), just as at (3.7). (We will show in Lemma 4.4 that the denominator w(Q|A, S) is again finite, and indeed depends on A and S only through the particle counts k 1 , . . . , k N .) The weight has the following interpretation. At each time slot where a service occurs, this service is offered to each of the customers, in order of priority (starting with those of type 1, then those of type 2, and so on). Each customer accepts an offer with probability 1 − q and rejects it with probability q, with the exception that any customer who has just arrived at the queue must accept it. The values of w i above are the corresponding probabilities for observing a particular type of departure or unused service, taking into account the current composition of the queue and the type of arrival (if any).
Then the weight w yields a distribution on departure configurations D; namely, choose A from the (N − 1)-type stationary distribution, and independently choose S uniformly, as described above; choose Q in proportion to P (Q|A, S) (equivalently, in proportion to w(Q|A, S); finally, take D = D(A, S, Q).
Theorem 3.4. The distribution of the configuration D = D i , i ∈ Z L induced by the weight (3.9) on queue-length processes is the stationary distribution of the N -type ASEP on Z L with particle counts k 1 , . . . , k N .

Proofs of main results for systems on Z L
In this section, we prove Theorem 3.3 and Theorem 3.4 concerning the construction of the ASEP on Z L (in the 2-type and general N -type cases respectively), and then Theorem 1.1 and Theorem 1.2 justifying Algorithms 1 and 2.

Proofs of Theorem 3.3 and Theorem 3.4
We start with the 2-type case of Theorem 3.3. We explain this case somewhat thoroughly, so as to indicate as clearly as possible the extension of the argument to N types, where the notation is more complicated and fewer details will be included. The form of the matrix product solution in the case N = 2 is particularly simple. From (2.2) and (2.3), we obtain We work specifically with the forms of α, and δ given at (2.5).
> 0, then there exists a unique pair A, S such that (A, S, Q) ∈ R (2) and D = D(A, S, Q).
Proof. Note that the non-zero entries of the matrices X r , r = 1, 2, ∞, are precisely the diagonal and super-diagonal entries of X and the four lines of (3.6). For part (i), we can verify that when any one of these four cases hold, the relevant matrix entry, namely the (Q i , Q i+1 ) entry in X D i , equals the weight defined in (3.6).
For part (ii), if we are given D, Q such that for each i one of the lines in (4.2) is satisfied, then defining A i and S i for each i according to the corresponding line of (3.6) satisfies (3.5), so that (A, S, Q) ∈ R (2) and D = D(A, S, Q) as required.
Finally, from (3.5) it's also easy to see that for a given D, Q there could not be more than one pair A, S such that (A, S, D, Q) ∈ R (2) (since the first line of (3.5) determines A, and the second and third lines determine S).
We will need one further property, that the denominator in (3.7) does not depend on A and S. This property is a central (and perhaps non-obvious) part of the argument. Its proof will do much of the work needed for the justification of Algorithm 1 in Theorem 1.1, and is given below. Meanwhile we complete the proof of Theorem 3.3: Proof of Theorem 3.3. Write P (2,L) k 1 ,k 2 for the distribution described by (3.6) and (3.7). Take any configuration D with k 1 first-class and k 2 second-class particles. We restrict to A and S with i 1(A i = 1) = k 1 and i 1(S i = 1) = k 1 + k 2 . Then Then there exists a unique S such that (A, S, Q) ∈ R (N ) and D = D(A, S, Q).
Proof. Similarly to the proof of Lemma 4.1, we need to show that the lines of (3.9) correspond to the cases where the entries in the tensor products defined in (2.3) are non-zero. This correspondence will be line-by-line between the first six lines of (2.3) and the six lines of (3.9). We will not go through every case, but will give an example.
First consider the four matrices I, α, δ and involved in the tensor products in (2.3), using as before the explicit forms at (2.5). The non-negative entries of I and α are precisely the diagonal entries; those of δ are the entries on the subdiagonal, and those of are those on the superdiagonal.
Then consider for example the fourth line of (2.3), which is the case of a (N ) m,n for n < m ≤ N − 1. Inspecting the tensor product, we see that the subdiagonal matrix δ appears in the nth component, and the superdiagonal matrix appears in the mth component; all the other matrices involved are diagonal. So the only non-zero terms a are those where Q i+1 = Q i − e n + e m . This corresponds to the fourth line of (3.9). We then need to check that the weight defined in the fourth line of (3.9) is the same as the entry in the tensor product. The components r = 1, 2, . . . , n − 1 of the tensor product contribute values α Q (r) The nth component contributes the value δ Q (n) i . The remaining components all contribute the value 1 from the relevant entries in the matrix or I. Multiplying together, we indeed obtain the weight defined in the fourth line of (3.9) as required.
In a similar way each of the first six lines of (2.3) accords with the corresponding line of (3.9). So for part (i), we have that if (3.8) is satisfied, then one of the lines of (3.9) holds, and then the weight defined by that line is the same as the tensor entry in the corresponding line of (2.3).
For part (ii), suppose we have A, D, Q such that the tensor entry a is positive for each i. That is, for each i, one of the lines of (2.3) is positive; then one can verify as above that if we take S i = 1 if D i < ∞ and S i = ∞ if D i = ∞, the corresponding line of (3.9) also holds, and that this is the only such choice of S i . This choice of S then satisfies (3.8) (and is the only such choice).
The corresponding result to Lemma 4.2 is: Again this result is a central part of the subsequent argument, and we prove it below.
Proof of Theorem 3.4. Write P N,L k 1 ,...,k N for the distribution resulting from the weight defined at 3.9. Let D be any configuration with particle counts k 1 , . . . , k N . Note that if D = D(A, S, Q) then A has particle counts k 1 , . . . , k N −1 and i 1(S i = 1) = k 1 + · · · + k N . We have  We start with the 2-type case. Throughout, we fix some arrival process A with k 1 arrivals, and some service process S with k 1 + k 2 potential services. Let S be the set of times where service is offered; that is, S = {i ∈ Z L : S i = 1}.
We first define a system of random marks attached to the service times, which will be used to determine which services are used and which remain unused. Let B i , i ∈ S be i.i.d. geometric with parameter q; that is, B i takes value b i ∈ {0, 1, 2, . . . } with probability q b i (1 − q). We write P B for the joint law of the B i , so for a given vector Recall that under our queueing model, if there is no arrival at i and the current queue-length Q i is k, the service at i is unused with probability q k and used with probability 1 − q k . We will arrange that the service at i is used if B i < Q i , and unused if B i ≥ Q i . Hence B i may be interpreted as the number of times that the service at i is refused by a customer. Suppose Q is a valid queue-length process (i.e. (A, S, Q) ∈ R (2) ). We will say that Q is compatible with a mark b i at i, where S i = 1, if one of the following three conditions holds: If Q is compatible with b i for each i ∈ S, we simply say that Q is compatible with the collection b = {b i , i ∈ S}, and we write Q ∼ b. Now we may rewrite the definition of w around (3.6) as follows. Given A, S and a valid queue-length process Q, w(Q|A, S) = P B (Q ∼ B). Proof. Let D = D(A, S, Q) and D = D(A, S, Q ) be the departure processes associated to Q and Q respectively. Without loss of generality, suppose that for some i, Q i ≥ Q i . We wish to show that also Q i+1 ≥ Q i+1 . Consider two cases. Suppose Q i > Q i . Since both processes have the same arrivals, and there is at most one departure at any time, certainly also Q i+1 ≥ Q i+1 . If instead Q i = Q i then from the fact that Q and Q are compatible with the same marks b, we have D i = D i and hence also Q i+1 = Q i+1 .
Hence we obtain Q i ≥ Q i for all i. Again since Q and Q are compatible with the same marks, we then obtain that D has a departure whenever D has a departure. But the total number of departures is k 1 in each process, so the set of departure times is the same for both. Since the set of arrival times is also the same, the queue-length processes differ by some constant.
Before stating the next result, we formulate two further equivalent versions of Algorithm 1. Recall that Algorithm 1 considers the k 1 arrivals in turn; when considering the (r + 1)st arrival, it chooses between the k 1 + k 2 − r services still available. If a service coincides with the arrival, that one is chosen. Otherwise, we number the available services i 1 , i 2 , . . . , i k 1 +k 2 −r in cyclic order, and service i j is chosen with probability q j−1 /(1 + q + q 2 + · · · + q k 1 +k 2 −r−1 ).
Equivalently, we can do the following. Now, when an arrival does not find an available service coinciding with it, we extend the list of available services i 1 , . . . , i k 1 +k 2 −r cyclically to an infinite sequence i j , j ≥ 1, by setting i j = i j whenever j ≡ j mod k 1 + k 2 − r. Now we choose service i j with probability q j−1 /(1 + q + q 2 + . . . ) = q j−1 (1 − q). That is, we choose service i J where J is a geometric random variable with parameter q. We have the following interpretation: we "offer" the arrival to each service i j in turn, proceding in cyclic order around the ring. Each service accepts the arrival with probability 1 − q. If a service accepts the arrival, then the arrival is assigned to that service; otherwise we continue to the next service in the list. If we have toured all the way around the ring without any offer being accepted, we continue in cyclic order, starting again at the beginning. It's easy to check that this procedure gives the same distribution as the one in the previous paragraph.
Secondly, we can allow a set of marks b = (b i , i ∈ S) to control the procedure. As the algorithm proceeds, a service at i will "reject" the first b i times it receives an offer, and "accept" the next time (if it is ever offered that many). If the marks b i are randomly chosen, with distribution i.i.d. geometric with parameter q, then this is again equivalent to accepting each offer independently with probability 1 − q. Proof. The second reformulation above can also be used to give a process Q which is compatible with a given set of marks b. We generate this Q by considering separately the contribution of each customer in the system. Suppose that a given customer arrives at time i and departs at time j. This customer contributes 1 to the queue lengths at times Q i+1 , Q i+2 , . . . , Q j−1 , Q j . In addition, suppose that the match between this customer and the service at j was rejected of r times before being accepted. (In that case, the algorithm has done j "complete tours" of the ring in the process of assigning the departure time of this customer). Then this gives a further contribution of r to the queue-length Q i at all times i .
We claim that summing up these contributions from each customer gives the desired queue-length process. Each arrival gives an increase of 1 in the queue-length, and each departure gives a decrease of 1. Further, each time that a service at i is unused by a customer, that customer gives a contribution of 1 to the queue-length at i; as a result, if the service remains unused at the end of the procedure, we know that Q i ≤ b i . In addition, any departing customer still contributes at the moment of departure; if a service at i is used for a departure after being rejected b i times, we know that Q i ≥ b i + 1. From these properties it follows Q is a valid queue-length process for (A, S), and also the Q is compatible with the marks b.
Hence there is at least one compatible queueing process. Also Lemma 4.5 states that any two such processes differ by a constant. Finally, all processes are nonnegative. So indeed a minimal such process exists. Proof. Suppose Q = Q min (B). From Lemma 4.5, if also Q ∼ B then Q and Q must differ by a constant, say Q i = Q i + m for all i, and since Q is minimal, we must have m ≥ 0. So it's enough to show that in that case P(Q ∼ B|Q = Q min (B)) = q mk 2 . From (4.3), we have that if Q = Q min (B), then B i ≥ Q i for all i ∈ S such that A i = ∞ and D i = 2, and further that B ∼ Q also holds iff B i ≥ Q i + m for all such i. For each such i, we have P(B i ≥ Q i + m|B i ≥ Q i ) = q m , since B i is geometric with parameter q. Since the variables B i are independent, and since there are k 2 such i, the overall conditional probability equals q mk 2 as required.
At this point we can deduce the fact that the total weight of all processes Q depends on A and S only through the total number k 1 of arrivals and k 2 of services, which was an important element of the argument in Section 4.1: Proof of Theorem 1.1. We wish to show that the distribution on configurations obtained from Algorithm 1 (or either of the variants described before Lemma 4.6) is the same as that given by (3.6) and ( As observed at (4.4), w(Q) = P B (Q ∼ B). We now decompose this weight by writing w(Q) = Q w(Q, Q ), where we define w(Q, Q ) = P(Q ∼ B, Q = Q min (B)). For all Q and Q such that this weight is positive, Lemma 4.5 gives that Q and Q differ by a constant and that D(Q) = D(Q ). Hence equivalent to (3.6) and (1.2) is the following: choose the pair Q and Q in proportion to the weight w(Q, Q ), and then take D(Q ).
Next we writew(Q ) = Q w(Q, Q ). Then we have another equivalent version: choose Q in proportion to the weightw(Q ), and take D(Q ). (4.7) Since 1/(1−q k 2 ) is a constant, we have another equivalent version: choose Q with probability proportional to the quantity P B (Q = Q min (B)) and take D(Q ). But of course this is the same as the following; generate a sample of the weights B, and take D(Q min (B)).
Finally, using Lemma 4.5 again, we may do the following: generate B, choose any Q such that Q ∼ B, and take the departure process D(Q). But, by the argument preceding Lemma 4.6 above, this is equivalent to what Algorithm 1 does.
Hence indeed Algorithm 1 leads to the same distribution as (3.6) and (1.2), as desired.
This completes the proof of the two-type result in Theorem 1.1. The structure of the proof of the multi-type result in Theorem 1.2 is entirely analogous. We outline the generalisation of the argument.
As before, we consider marks B i at each of the service times i, which are i.i.d. geometric with parameter q. The mark B i represents the number of times that the service at time i may be rejected. A multi-type process Q is compatible with a set of marks b (for which we write Q ∼ B) if each of its embedded one-type queues is compatible with b. That is, for each n and for each i such that S i = 1, one of the following conditions holds, just as at (4.3): • A (≤n) (i) = 1;  The equivalent of Lemma 4.6 also holds, so that for any set of marks b there exist queueing processes compatible with b and among them a minimal process Q min (b) among those queueing processes compatible with the marks b. Here the process Q min is minimal in the sense that for any other compatible Q and any n and i, Q (≤n) min (i) ≤ Q (≤n) (i). (It may not necessarily be the case that Q (n) min (i) ≤ Q (n) (i)). Indeed, the multi-type Q min can be obtained by taking each Q (≤n) min to be the minimal process for the embdedded one-type queue, as given by Lemma 4.6.
The following generalisation of Lemma 4.7 then holds: Lemma 4.9. Let A be an arrival process with particle counts k 1 , . . . , k N −1 . Let S be a service process with k 1 + · · · + k N services. Let Q and Q be two valid queue-length processes. If there exist m n ≥ 0, 1 ≤ n ≤ N − 1 such that for all n and i, Otherwise To complete the proof of Theorem 1.2, start from w(Q|A, S) = P B (Q ∼ B). Define again w(Q, Q ) = P B (Q ∼ B, Q = Q min (B)) and thenW (Q ) = Q w(Q, Q ). The equivalent step to (4.7) is that The denominator in the final expression is constant, and the end of the proof goes through just as before.

Proof of the results on common denominators
As noted in Section 1, we can generate a sample from the N -type ASEP on Z L by applying Algorithm 2 N − 1 times in all. The nth iteration takes as arrival process a configuration whose distribution is stationary for an n-type ASEP, along with an independent service process, and outputs a configuration whose distribution is stationary for an (n + 1)-type ASEP. The N − 1 iterations can be combined into a "multi-line" diagram with N lines. In particular, the bottom line of the diagram gives a sample from an N -type system. We can identify two sources of randomness within this procedure; the choice of the configuration of occupied and unoccupied sites on the lines of the multi-line diagram, and the assignment of types to the occupied sites. The set of occupied sites on line n is chosen uniformly from all subsets of Z L of size k 1 + · · · + k n , independently for different lines; hence overall the configuration is uniformly chosen from possibilities.
To assign types to particles in line r, given the types in line r−1, we use Algorithm 2. We treat the arrivals in order of their type (within a single type, the order in which we treat the arrivals is irrelevant -for example, we can work from left to right). If m of the arrivals (from line r − 1) have been assigned to service-times (on line r), then there remain k 1 + · · · + k r − m services. When we come to assign the (m + 1)st arrival, there are two cases. Either there is an unassigned service immediately below it, in which it is assigned there. Otherwise, the choice between the remaining services is done according to a distribution with weights 1 1 + q + · · · + q k 1 +···+kr−m−1 1, q, . . . , q k 1 +···+kr−m−1 . (4.11) Once the set of occupied sites is fixed, the probability of obtaining any particular multiline diagram is given by a product of terms of the form (4.11), where r, m vary over 1 ≤ m ≤ k 1 +k 2 +· · ·+k r−1 . To obtain the polynomial (1.5) we take the product over such r, m of the denominator [k 1 + · · · + k r − m − 1] q from (4.11), and multiply by (4.10). Every numerator term in (4.11) is a power of q, so indeed the probability of any multiline diagram is a polynomial in q with non-negative integer coefficients divided by the polynomial (1.5). But the bottom line of the multiline diagram gives a sample from the stationary distribution of the N -type ASEP. So the stationary probabilities are sums of probabilities of multiline diagrams, and so themselves have the same form, giving Theorem 1.3.
For Corollary 1.4, we consider an (L − 1)-type system with k r = 1 for r = 1, 2, . . . , L − 1. (The single hole plays the role of the particle of type L.) In this case (1.5) reduces to (1.7).

Systems on the line
Now we want to consider systems in which the set of sites is given by the whole integer line Z rather than the ring Z L .
Given λ 1 , . . . , λ N with λ i < 1, there is a unique ergodic translation-invariant stationary distribution where the density of particles of type i is λ i , as shown by Ferrari, Kipnis and Saada [25].
We want to prove Theorem 3.1 for the 2-type case and Theorem 3.2 for the multitype case. In each case we will prove stationarity for the measure on the whole line by showing that it is the limit of the stationary distribution of systems on finite rings.
For T ∈ N, consider systems on the ring of size L = 2T , with the sites labelled as −T, . . . , T − 1. Let ν 2T be any sequence of stationary distributions for these systems. If C is a cylinder event which depends on a set of sites −m, . . . , m, then we can consider the probability of C under ν 2T for any T > m.
Proposition 5.1. If the distribution ν satisfies ν 2T (C) → ν(C) as T → ∞ for every cylinder event C, then ν is a stationary distribution for the system on Z.
Proof. Let C depend on the sites −m, . . . , m and fix any t > 0. Consider a system on the ring of size 2T started from some state η (T ) 0 ∼ ν 2T at time 0, and a system on Z started from some state η 0 ∼ ν at time 0.
Since ν 2T → ν on cylinder events, for any given M , whenever T is large enough, we can couple η (This can be achieved by a simple local coupling of the dynamics. We omit the details but a straightforward approach is as follows: the process of jumps between sites i and i + 1 can be dominated by independent Poisson processes P x of rate 1 + q. For any given i, the probability that P i has no point in the time interval [0, t] is e −t(1+q) ; hence with high probability as M → ∞, we can find such sites i 0 ∈ {−M, . . . , −m − 1} and i 1 ∈ {m, . . . , M − 1}. Then one can couple such that in both systems, no particle enters or leaves the interval {i 0 + 1, . . . , i 1 } during [0, t], and such that inside that interval, the evolutions of the two systems are identical. In particular, the two evolutions coincide on the window {−m, . . . , m}.) Hence, restricted to −m, . . . , m, the distribution of the time-t state in the infinite system is arbitrarily close to that of the time-t of the system on the ring. But on the ring, the distribution of the time-t state is ν 2T , since ν 2T is stationary. Also by choosing large enough T , ν 2T can be made to agree arbitrarily closely with ν on the window −m, . . . , m. Hence the probability of the event C occurring at time t in the system on Z must be ν(C). This holds for all cylinder events C, so indeed the distribution at time t is ν; hence ν is a stationary distribution as required.
Now we turn to the queueing Markov chain. We will give it a rather explicit construction, involving data like the "marks" that we used in the previous section when justifying Algorithms 1 and 2.
We will start by giving complete details in the case N = 2, where the chain consists of a single queue. Then we will indicate the extension to larger N (where the chain consists of several queues in tandem, and so the details are somewhat more complicated to write down in full, although entirely analogous).

N = 2 (single-type queue)
First we introduce a rather concrete representation of the Markov chain used in Theorem 3.1. We consider random vectors R i = (A i , S i , B i ), i ∈ Z, which control the evolution of the queue. All entries are independent, with P(A i = 1) = 1 − P(A i = ∞) = λ, P(S i = 1) = 1 − P(S i = ∞) = µ, B i geometric with parameter q.
Here A i = 1 means that there is an arrival at i. S i = ∞ means there is no service available at i. If on the other hand S i = 1, then there is a potential service available at i, which is used if there are at least B i customers in the queue, or if an arrival has just occurred; otherwise it is unused.
Writing Q i for the queue-length at the beginning of time-step i, we can write a recursion of the form for some appropriate function f (which would be easy to write down, but the precise form is not important). This formal representation is useful because it allows us to couple versions of the system evolving from different initial conditions but using the "same randomness" (R i , i ∈ Z), and, particularly, to couple versions of the system evolving on the whole real line Z or on a finite box [−T, T ]. Similarly, we can write the departure process defined at (3.2) in the form again for some suitably chosen function g.

Cyclic evolution
Let G T be the event that there are more services than arrivals in the finite box; specifically, that T −1 When G T holds, we saw in the proof of 4.6 how to construct a cyclic evolution compatible with the randomness; that is, an evolution Q cyc −T , . . . , Q cyc T with the property that Q cyc −T = Q cyc T , and such that for all i = −T, . . . , T − 1, Q cyc i+1 = f (Q cyc i , R i ) as in (5.1) holds. In Lemma 4.5 we saw that any two such cyclic evolutions differ by a constant, and that furthermore any two such evolutions give the same output configuration D cyc = (D cyc i , i ∈ [−T, T − 1]) (where D cyc i = g(Q cyc i , R i ) as at (5.2)). Conditional on the number of arrivals and services, the distribution of this configuration is stationary for the ASEP on Z 2T . Let ν 2T λ,µ−λ be the distribution of D cyc , conditioned on the event G T . Then ν 2T λ,µ−λ is a mixture of stationary distributions for the ASEP (by Theorem 3.3), and so is itself stationary.

Coupling of cyclic evolution to an evolution on Z
The queue-length Markov chain described in Section 3.2 has a stationary distribution π. (This distribution is easy to obtain using the detailed balance equations, but we don't need its particular form here.) Suppose Q eq −T is drawn from the stationary distribution π of the queue-length Markov chain, independently of the randomness (B i , A i , S i ), i = −T, . . . , T − 1. Then define Q eq i , i = −T + 1, . . . , T by Q eq i+1 = f (Q eq i , R i ) as at (5.1). We call this the equilibrium evolution. Similarly define D eq i , i = −T + 1, . . . , T as at (5.  So Q also evolves as a copy of the queue-length chain on [−T, T ], using the same randomness as Q eq , but started from the particular initial state Q eq T . We claim that with high probability, Q eq and Q couple together on the interval; that is, they reach the same state at some point, and then since they use the same randomness, they stay together for the rest of the interval also. Then Q T = Q eq T = Q −T , and so in particular Q is a cyclic evolution. In fact, we aim to show that with high probability, the coupling point occurs before −m, so that also Q i = Q eq i for all i ∈ [−m, m]. See If this event holds then Q eq i = Q i for all t ≤ i ≤ T , and in particular: If both G T and the event in (ii) occur, then the equilibrium departure process D eq and the cyclic departure process D cyc agree on the window [−m, m]. Since m is arbitrary, we obtain from Proposition 5.2 together with (5.3): Corollary 5.3. For any cylinder event C, ν 2T λ,µ−λ (C) → ν (2) λ,µ−λ (C) as T → ∞. Hence by Proposition 5.1, ν is a stationary distribution for the system on Z as required.
In the rest of this section we prove Proposition 5.2. We want to consider a "pair chain", which consists of two copies of the chain evolving with the same randomness.
If > 0, then there exists t 0 such that Proof. The underlying Markov chain, with updates given by (5.1), is irreducible, aperiodic, and positive recurrent, with a unique stationary distribution π. Fix some k such that k x=0 π(x) > 3/4. Starting from any initial state, we can apply the Markov chain ergodic theorem to deduce that with probability 1, the long-run proportion of time which the chain spends in [0, k] is more than 3/4. Now consider the pair chain with updates given by (5.4). Starting from (x 0 , x 0 ), with probability one, eventually each coordinate of the chain spends at least 75% of its time in [0, k]. Hence the pair-chain spends at least half its time in the set In particular, it visits that set infinitely often.
Consider the pair chain evolving from any state in H k . With some uniformly positive probability, after k more steps the chain is in the state [0, 0]. For example, it suffices that at each of the next k steps, there is no arrival, and a service which is accepted by the first customer in the queue (if any is present). This happens with probability at least [(1 − λ)µ(1 − q)] k .
Since the chain visits H k infinitely often with probability 1, and from any state in H k there is uniformly positive probability to reach the state (0, 0) within k steps, it follows that with probability 1 the chain will eventually visit state (0, 0). So certainly there is some time t 0 such that with probability 1 − , the chain visits (0, 0) before time t 0 . But once it visits (0, 0), the two coordinates of the chain stay the same for ever.
Proof of Proposition 5.2. In the setting of Proposition 5.2, we start the pair-chain from a particular state (Q eq −T , Q −T ). We wish to show that for fixed m, the probability that the two coordinates of the chain couple before time −m, i.e. within T − m steps, tends to 1 as T → ∞.
Note that Lemma 5.4 tells us that started from any fixed state, there is a t 0 sufficiently large that with probability at least 1 − , the two coordinates of the chain couple within t 0 steps. It remains to deal with the fact that the initial condition (Q eq −T , Q −T ) is random rather than fixed.
Note that (Q eq t , t ∈ Z) is an equilibrium version of the chain. In particular, both Q eq −T and Q −T = Q eq T have distribution π. For any given δ, choose k sufficiently large that ∞ x=k+1 π x ≤ δ/3. Then Now using Lemma 5.4, there is some t 0 such that, started from any given state in [0, k] × [0, k], the probability of failing to couple within t 0 steps is at most δ/3(k + 1) 2 . Then since there are (k + 1) 2 possible initial states in [0, k] × [0, k], the probability that there exists at least one initial state in [0, k] × [0, k] such that failure to couple occurs (using the updates (5.4)) is at most δ/3.
Then for all T large enough that T − m > t 0 , indeed with probability at least 1 − δ, the pair-chain started from state (Q eq −T , Q −T ) at time −T , and updated using (5.4), couples before time −m. Since δ is arbitrary, this gives Proposition 5.2 as required.
Hence we have Corollary 5.3, and via Proposition 5.1 we have justified the construction of the 2-type stationary distribution on Z in Theorem 3.1.

N > 2 (multi-type queue)
Everything in the previous section generalises straightforwardly to the case N > 2, to prove Theorem 3.2 constructing the stationary measure of the N -type system on Z. Let us outline the changes.
The queueing Markov chain is now a system of queues in tandem. There are N −1 queues, and the departure process of queue r is the arrival process for queue r + 1, for 1 ≤ r ≤ N − 2. (Note that a customer leaving queue r at time i arrives at queue r + 1 in the same time-slot i, so may pass through queue r + 1, and further queues, in the same time-slot.) The state of the system at time-step i may be written as where Q (n,r) i is the number of customers of type n in queue r at the beginning of time-step i. The evolution of the queues can again be controlled by a vector of random variables corresponding to each time-step. Now we would have for some appropriate function f , and this again allows us to couple versions of the system starting from different initial conditions.
The equivalent of the event G T defined at (5.3) is now that, on the time interval [−T, T − 1], queue 1 has more services than arrivals than services, and queue r + 1 has more services than queue r does for each r = 1, 2, . . . , N − 2. Again the probability of this event tends to 1 as T → ∞.
The Markov chain is now multi-dimensional, and its stationary distribution is no longer easy to obtain. However, it is still positive recurrent (as can be seen by considering the marginal evolution of each queue individually; for any given r, the process of the total number of customers in the rth queue, namely r n=1 Q (n,r) i , i ∈ Z, behaves like a single-type queue with arrivals at rate r n=1 λ n and potential services at rate r+1 n=1 λ n ). Then exactly the same ideas for coupling the cyclic evolution and the equilibrium evolution apply, as illustrated in Figure 5.1.
Define From this cylinder convergence, using Proposition 5.1, we get the result of Theorem 3.2 as required.

Clustering
In this section we prove Theorem 1.5. Recall that Y (L) denotes a configuration drawn from the stationary distribution of an ASEP with L types on the ring with L sites − L 2 , . . . , −1, 0, 1, . . . , L 2 − 1, and as at (1.8), we put We want to show that W (L) converges in distribution as L → ∞. Some of the arguments we need involving convergence of stationary distributions on the ring to those on the line have already been developed in the previous section during the proof of Theorems 3.1 and 3.2. We begin with a variation of Proposition 5.5.
Proposition 6.1. Fix λ 1 , . . . , λ N ∈ (0, 1) with λ 1 + · · · + λ N < 1. Proof. We indicate the differences between this result and Proposition 5.5 above. A first and rather trivial difference is that we no longer assume L is even. In the previous section we took L = 2T throughout, but this was purely for notational convenience and makes no difference to the method of proof.
The more substantial difference is that now we take the number of particles of each type to be deterministic, rather than given in terms of binomial random variables as in Proposition 5.5.
Similarly, for any x 0 , x 1 , . . . , x r−1 ∈ (0, 1), the probability of any event of the form converges as L → ∞ to a limit that does not depend on i (which may be written in terms of a suitable r-type stationary distribution for the ASEP on Z).
More generally for (6.2), let (u 1 , . . . , u N ) be the increasing reordering of (x 1 , . . . , x N ). Then we can apply Proposition 6.1 with λ n = u n − u n−1 , n = 1, . . . , N , to obtain the limit in (6.2) in terms of the N -type stationary distribution ν λ 1 ,...,λ N . Since that distribution is translation invariant on Z, the limit is the same for all i.
Convergence of all events of the form of (6.2) is enough to show that the sequence W (L) has a single distributional limit point, and this gives the distribution of W required for part (a) of Theorem 1.5. Since W (L) i has uniform distribution on {1/L, 2/L, . . . , 1}, we have that W i ∼ Uniform[0, 1]. In fact, with little more work we could obtain that the distribution of W is the unique ergodic translation-invariant distribution for the ASEP on Z whose marginals have Uniform[0, 1] distribution (but we won't need to use this fact directly).
For parts (b) and (c), we will analyse the distribution of W via its projections onto 2-type systems as in (6.1), using the queueing construction of the stationary distribution ν λ 1 ,λ 2 of the 2-type ASEP on Z. Let us recall that construction. We have an stationary queueing system in discrete time. At each time i ∈ Z, we have an arrival with probability λ 1 , and independently a potential service with probability µ := λ 1 + λ 2 . When a potential service occurs, if an arrival has occurred at the same time-step, then a departure occurs with probability 1; if no arrival occurred, then a departure occurs with probability 1 − q k where k is the number of customers present in the queue, and otherwise (with probability q k ) an unused service occurs. Then set η i = 1 if a departure occurs at i, η i = 2 if an unused service occurs at i, and η i = ∞ if no service was available at i. At a given time i, the marginal probability that η i = 1 is λ, that η i = 2 is µ − λ, and that η = ∞ is 1 − µ.
First, for part (c), we use arguments similar to those used for related calculations concerning the TASEP speed process in the q = 0 case in [2], although the analysis is more complicated now that q > 0, since unused services can occur even when the queue is not empty.
Recall that we write f for the density of (W 0 , W 1 ) on {(x, y) ∈ [0, 1] 2 : x = y}, and f * for the density along the diagonal x = y. Lemma 6.3.
Note in particular that as µ − λ → 0, the density of second-class particles goes to 0, but conditional on seeing a second-class particle at a given site, the probability that the next site also contains a second-class particle stays bounded away from 0.
Part (i) is straightforward. The second or third part are then easily seen to follow from each other, but are more complicated to establish. Note that when q = 0, an unused service can only occur when the queue is empty, and the conditional probability of the output of the next time-slot is then easy to deduce. For q > 0 on the other hand, an unused service occurs with probability q n when the queue-length is n. As a result, conditioning on an unused service constitutes an exponential tilting of the stationary distribution of the system.
Proof of Lemma 6.3. For part (i), we need the probability of observing an unused service at time 0, followed by no service at time 1. The probability of an unused service at time 0 is µ − λ, and a service occurs at time 1 with probability µ independently of everything that has occurred earlier, giving (µ − λ)(1 − µ) overall as required.
We turn to part (ii). Since the probability of an unused service at time 1 is µ − λ, we need to show that conditional on an unused service at time 1, the probability of another unused service at time 2 is First let's consider the process of the length of the queue after each service, which is a Markov chain. It's a birth-and-death chain on Z + , and for all k > 0, Hence its stationary distribution (π k , k ≥ 0) satisfies The probability of seeing a second-class particle at a given site is k π k (1 − λ)µq k , while the probability of seeing two such particles at a given pair of successive sites is So we wish to show that (6.4) is equal to k π k (1 − λ) 2 µ 2 q 2k k π k (1 − λ)µq k . We claim that for any α < 1 and q < 1, the following identity holds: The desired equality then follows from (6.6) with α = λ 1−λ 1−µ µ . To obtain (6.6), we may write the difference between the RHS and the LHS as This establishes part (ii). Now part (iii) follows because the expressions in (i), (ii) and (iii) must sum to the probability of observing a second-class particle at site 1, which is µ − λ.
Finally, for part (iv), the claimed expressions for the density f (x, y) are continuous away from x = y and the expression for f * (x) is continuous in x. Then for f (x, y), it suffices to check that for y < x, dxdy ν y,x−y (2, 1), and similarly for x < y, d 2 dxdy P(W 0 ≤ x, W 1 ≤ y) = d 2 dxdy P(x < W 0 ≤ y, W 1 > y) = d 2 dxdy ν x,y−x (2, ∞).
For part (b) of Theorem 1.5, we obtain (1.9) and (1.10) by substituting the expressions in parts (i)-(iii) of Lemma 6.3 into part (iv). Then (1.11) follows by integrating over x and y in (1.9), and over x in (1.10).
Finally we turn to part (c) of Theorem 1.5, which is the most involved. It follows from the following result: Proposition 6.4. Define a random set U ⊂ {1, 2, . . . } in the following way.
Then the set U is infinite with probability 1, and is stochastically dominated by the set {i ≥ 1 : W i = W 0 }.
Here Q plays the role of a queue-length process, and U plays the role of the set of unused service times, in a queue whose arrival rate and service rate are both equal to x. (Such a queue is null recurrent.) Although here we only claim that U is a stochastic lower bound for the set {i ≥ 1 : W i = W 0 }, it actually holds that the two have the same distribution. In fact, we can be more precise: there is a system of regular conditional probabilities such that conditional on W 0 = x, the distribution of {i ≥ 1 : W i = W 0 } is the distribution of U obtained by the construction of Proposition 6.4 given X = x. The argument to establish these stronger statements may be apparent from the proof below, but we won't fill in the details.
The rest of this section is devoted to the proof of Proposition 6.4. First we compare the construction in Proposition 6.4 to the queue-length construction of the 2-type stationary distribution on Z: Lemma 6.5. Fix any λ, µ with 0 ≤ λ < µ ≤ 1.
Write U x,x for a set whose distribution is that of the set U in Proposition 6.4 conditioned on X = x. For λ ≤ x ≤ µ, the set U λ,µ stochastically dominates the set U x,x .
Since λ < µ, the distribution defined by (6.9) is dominated by the distribution defined by (6.7). Furthermore, the "up-step" probability λ(1 − µ) in (6.10) is smaller than the corresponding probability x(1 − x) in (6.8), and the down-step probability (1 − λ)µ(1 − q r ) is larger than the corresponding probability x(1 − x)(1 − q k ) when r = k. It follows that there is a coupling of the walks (Q i , Q λ,µ i , i ≥ 1) which always stays in the set {(k, r) : k ≥ r}.
Further, the probability of including the next time-step i in the set U x,x (or U λ,µ respectively) is smaller at (6.8) than at (6.10) whenever k ≥ r, since then x(1 − x)q k < (1 − λ)µq r . From this it's straightforward to extend the coupling in the previous paragraph to a coupling of (Q, Q λ,µ , U x,x , U λ,µ ) such that with probability 1, U x,x ⊆ U λ,µ , as required.
Lemma 6.6. The distribution of the set U λ,µ defined in Lemma 6.5 is the same as the distribution of the set {i > 0 : λ ≤ W i < µ} conditional on the event {λ ≤ W 0 ≤ µ}.
Proof. Consider the projection f from (6.3), with u 1 = λ and u 2 = µ. The configuration h(W ) = (h(W i ), i ∈ Z) has the 2-type stationary distribution ν λ,µ−λ In particular, the sites of type-2 particles in f (W ) are precisely those i such that λ ≤ W i < µ.
The walk (6.10) is exactly the queue-length process with arrival rate λ and service rate µ used to generate the stationary distribution ν λ,µ−λ of the 2-type ASEP on Z, and the adding of a point to U λ,µ corresponds to an unused service in the queue.
Finally, the initial distributionπ λ,µ of Q λ,µ 1 is the distribution of the queue-length in equilibrium conditioned on an unused service having occurred at the previous timestep. To verify this, note that such a distribution should be proportional to q k π k where π k , the unconditioned equilibrium queue-length distribution, satisfies (6.5). This gives the recursion (6.9) as required.
Combining the last two lemmas, we have that the distribution of the set {i > 0 : λ ≤ W i < µ} conditional on the event {λ ≤ W 0 ≤ µ} dominates the distribution of the set U conditional on λ ≤ X < µ. But also the sets {i > 0 : |W i − W 0 | < 1/N } form a decreasing family as N increases. If they dominate U for every N , then also their intersection must dominate U. But the intersection is exactly the set {i > 0 : W i = W 0 }.
Finally note that U is infinite with probability 1, since conditional on any value of X = x ∈ (0, 1), the walk given by (6.8) is null recurrent, and each time the walk hits 0 the next site is added to the set U with probability x(1 − x) independently of all past choices.
This completes the proof of Proposition 6.4, and hence of Theorem 1.5.

Alternative queueing construction
The queueing discipline used in the constructions in this paper has a somewhat unnatural feature, in that a customer who has just arrived at the queue is treated differently from one who was already present since the previous time-step. A customer who has just arrived will always accept any offered service, while a customer who was previously present rejects any offer with probability q. One can also apply equivalent constructions with this distinction removed. Now each customer rejects each service with probability q, irrespective of time of arrival. We conjecture that this alternative model gives exactly the same distribution, and in particular realises the multi-type ASEP stationary distribution (both on Z and on Z L ).
This has been computationally verified for values of L up to 6. In the simplest case N = 2, it is not hard to verify more generally using a slightly different matrix product realisation. Rather than taking X One can then verify that these satisfy the following quadratic algebra of [19]: (The relation between these identities and (2.1) is that if δ, , α satisfy (2.1), then D = I + δ, E = I + and A = α satisfy (7.1).) To extend this from N = 2 to higher N along the lines of the proof above, we would need a variant of the result of Theorem 2.1 covering a different recursive system from the one presented at (2.3). The new system would have more non-zero terms (in particular, since we can now have a departure of lower priority than arrival, we would also have non-zero terms a (N ) m,n for all m < n ≤ N as in the first line of (2.3)). An alternative approach comes from arguments based on dynamic reversibility, involving dynamics defined on the set of multi-line diagrams, as done originally for the TASEP in [26]. Such arguments may also make it possible to extend to q > 0 the "generalised" multi-line construction described for the TASEP in [6], which extend the construction of [26], restoring the symmetry between particles and holes and establishing a richer "Hasse diagram" structure connecting systems with different numbers of particle types.