Parameter permutation symmetry in particle systems and random polymers

Many integrable stochastic particle systems in one space dimension (such as TASEP - Totally Asymmetric Simple Exclusion Process - and its various deformations, with a notable exception of ASEP) remain integrable when we equip each particle $x_i$ with its own jump rate parameter $\nu_i$. It is a consequence of integrability that the distribution of each particle $x_n(t)$ in a system started from the step initial configuration depends on the parameters $\nu_j$, $j\le n$, in a symmetric way. A transposition $\nu_n \leftrightarrow \nu_{n+1}$ of the parameters thus affects only the distribution of $x_n(t)$. For $q$-Hahn TASEP and its degenerations (namely, $q$-TASEP and beta polymer) we realize the transposition $\nu_n \leftrightarrow \nu_{n+1}$ as an explicit Markov swap operator acting on the single particle $x_n(t)$. For beta polymer, the swap operator can be interpreted as a simple modification of the lattice on which the polymer is considered. Our main tools are Markov duality and contour integral formulas for joint moments. In particular, our constructions lead to a continuous time Markov process $\mathsf{Q}^{(\mathsf{t})}$ preserving the time $\mathsf{t}$ distribution of the $q$-TASEP (with step initial configuration, where $\mathsf{t}\in \mathbb{R}_{>0}$ is fixed). The dual system is a certain transient modification of the stochastic $q$-Boson system. We identify asymptotic survival probabilities of this transient process with $q$-moments of the $q$-TASEP, and use this to show convergence of the process $\mathsf{Q}^{(\mathsf{t})}$ with arbitrary initial data to its stationary distribution. Setting $q=0$, we recover the results about the usual TASEP established recently in arXiv:1907.09155 [math.PR] by a different approach based on Gibbs ensembles of interlacing particles in two dimensions.

Ever since the original works on TASEP it was clear [ITW01], [GTW02] that integrability of some particle systems like TASEP is preserved in the presence of countably many extra parameters, e.g., when each particle has its own jump rate. We will refer to such more general systems as multiparameter ones. This notion should be contrasted with, e.g., the q-deformation (by means of just one parameter) taking TASEP to q-TASEP. The latter is much more subtle and relies on passing to a different class of symmetric functions (q-Whittaker vs Schur symmetric functions).
It should be noted that TASEP in inhomogeneous space (when the jump rate of a particle depends on its location) does not seem to be integrable [JL92], [Sep01], [CLST13] (cf. recent asymptotic fluctuation results [BSS14], [BSS17] which required very delicate asymptotic analysis). Moreover, it is not known whether ASEP has any integrable multiparameter deformations. The stochastic six vertex model [GS92], [BCG16] which scales to ASEP admits such a multiparameter deformation [BP18a], but this deformation is destroyed by the scaling. Recently other families of spatially inhomogeneous integrable stochastic particle systems in one and two space dimensions were studied in [BP18b], [KPS19], [Pet19], [Ass19].
All known multiparameter integrable stochastic particle systems display a common feature. Namely, certain joint distributions in these systems are symmetric under (suitably restricted classes of) permutations of the parameters. This symmetry is far from being evident from the beginning, and is often observed only as a consequence of explicit formulas. The main goal of the present paper is to explore probabilistic consequences of parameter symmetries in integrable particle systems.
1.2. Distributional symmetry of the q-Hahn TASEP. The most general system we consider is the q-Hahn TASEP introduced in [Pov13] and studied in [Cor14], [Vet15]. Its multiparameter deformation appears in [BP18a]. Under this deformation, each particle x n carries its own parameter ν n ∈ (0, 1) which determines the jump distribution of the particle. The distribution of x n (t) at any time moment t ∈ Z ≥0 in the q-Hahn TASEP started from the step initial configuration 2 depends on the parameters ν 1 , . . . , ν n in a symmetric way. (We recall the definition of the q-Hahn TASEP and check this symmetry in Section 3.2.) The main structural result of the present paper can be vaguely formulated as follows: 1 We say that an event in continuous time happens at rate α if P(waiting time till event occurs > t) = e −αt for all t ∈ R ≥0 .
That is, every site of Z<0 is occupied by a particle, and every site of Z ≥0 is empty.
Theorem 1.1 (Theorem 3.8 in the text). The elementary transposition ν n ↔ ν n+1 , ν n+1 < ν n , of two neighboring parameters in the q-Hahn TASEP started from step is equivalent in distribution to the action of an explicit Markov swap operator p qH n on the particle x n . The equivalence in distribution holds at any fixed time t ∈ Z ≥0 in the q-Hahn TASEP, while the swap operator p qH n does not depend on t. See Figure 1 for an illustration.
We prove this result in Section 3.3 using known q-moment contour integral formulas and duality results for the q-Hahn TASEP.
Remark 1.2. There are certain other classes of initial data (for example, half-stationary) for which the q-Hahn TASEP displays parameter symmetry. Moreover, via the spectral theory [BCPS15] one sees that for fairly general initial data the swap operators should be simultaneously applied to the particle system and the initial distribution. For simplicity, in this paper we focus only on the step initial configuration.
1.3. Applications. We explore a number of interesting consequences of the distributional symmetry of the q-Hahn TASEP realized by the swap operators. Let us briefly describe them. We take a continuous time limit of the q-Hahn TASEP and the swap operators. Denote by M qH q,ν;t the distribution of the parameter homogeneous (i.e., ν n ≡ ν), continuous time q-Hahn TASEP at time t ∈ R ≥0 started from step (see Section 4.2 for a detailed definition). The q-Hahn TASEP evolution acts on M qH q,ν;t by increasing the time parameter t. We find that a suitable continuous time limit of the swap operators produces a (time-inhomogeneous) process B qH which acts on M qH q,ν;t by simultaneous rescalings (ν, t) → (νe −τ , te −τ ) (here τ ∈ R ≥0 is the time parameter of the process B qH ). See Theorem 4.7 for a detailed formulation and Figure 6 for an illustration of the two actions. When ν = 0, the backward process becomes time-homogeneous, and we discuss this case in more detail in the next Section 1.4.
When q = ν = 0, Theorem 4.7 recovers one of the main results of the recent work [PS19] on the existence of a time-homogeneous, continuous time process mapping the distributions of the usual TASEP back in time. We remark that the proof of this result following from the present paper is completely different from the argument given in [PS19]. The latter went through the well-known connection of the TASEP distribution and a Schur process [OR03] on interlacing arrays (about this connection see, e.g., [BF14]). For Schur processes, the two-dimensional version of the swap operator is accessible by elementary means.
In a scaling limit q, ν n → 1, the q-Hahn TASEP turns into the beta polymer model introduced in [BC16a]. In Section 6 we construct swap operators for the multiparameter beta polymer model (this multiparameter beta polymer can be read off [BP18a]). The argument is formally independent of the rest of the paper, but proceeds through the same steps. The swap operator can be realized as a simple modification of the lattice on which the beta polymer is defined. See Theorem 6.4 for a detailed formulation of the result, and Figures 10 and 11 for illustrations of lattice modifications.
1.4. Stationary dynamics on the q-TASEP distribution. The last application concerns q-TASEP, which is a ν = 0 degeneration of the q-Hahn TASEP. Let us focus on this case in more detail. Denote by M qT q;t the time t distribution of the continuous time q-TASEP [BC14], [BCS14] started from the step initial configuration step. Here we assume that the q-TASEP is parameterhomogeneous (that is, a i ≡ 1 in the notation of [BCS14]). See Section 5.1 for the definition of the q-TASEP.
When ν = 0, Theorem 4.7 produces a new time-homogeneous, continuous time Markov process which we denote by Q (t) , with the following properties: • The process Q (t) is a combination of two independent dynamics: the q-TASEP evolution, and the (slowed down by the factor of t) backward q-TASEP evolution. The latter is a suitable degeneration of the backward q-Hahn process B qH . Under this degeneration, the backward process becomes time-homogeneous. See Section 5.2 for the full definition of the process Q (t) . • (Proposition 5.3) The process Q (t) preserves the distribution M qT q;t . (The time parameter t ∈ R ≥0 of the q-TASEP distribution is fixed and is incorporated into the definition of Q (t) .) • (Theorem 5.7) Start the process Q (t) from an arbitrary particle configuration on Z which is empty far to the right and densely packed far to the left. Then in the long time limit the distribution of this process converges to Q (t) .
We establish Theorem 5.7 by making use of duality for the stationary process Q (t) which extends the duality between the q-TASEP and the stochastic q-Boson process from [BCS14]. (The stochastic q-Boson process dates back to [SW98].) In fact, we are able to use the same duality functional (corresponding to joint q-moments) for Q (t) . As a result we find that the process dual to Q (t) is a new transient modification of the stochastic q-Boson process. The long time limit of this transient process is readily accessible, and Theorem 5.7 follows by matching the long time behavior of all q-moments of Q (t) (with arbitrary initial configuration) with those of the q-TASEP (with step initial configuration).
Let us illustrate the transient modification in the simplest case of the first q-moment. Consider the continuous time random walk n(t) on Z ≥0 which jumps from k to k − 1, k ≥ 1, at rate 1 − q. When the walk reaches zero, it stops. From the q-TASEP duality [BCS14] we have E qT step q xn(t)+n = P(n(t) > 0 | n(0) = n), n = 1, 2, . . . . (1.1) Here the left-hand side is the expectation over the q-TASEP distribution M qT q;t , and the right-hand side corresponds to the random walk n(t). Similarly to (1.1), joint q-moments of the q-TASEP are governed by a multiparticle version of the process n(t) -the stochastic q-Boson system.
Let us now fix the q-TASEP time parameter t ∈ R >0 , and consider a different random walk n (t) (τ ) on Z ≥0 (here τ is the new continuous time variable) with the following jump rates: This process has a single absorbing state 0, and otherwise is transient. In other words, after a large time τ , the particle n (t) (τ ) is either at 0, or runs off to infinity. Note however that this process does not make infinitely many jumps in finite time. The duality for the stationary process Q (t) which we prove in this paper states (in the simplest case) that Here the left-hand side is the expectation over the stationary process started from step, and the right-hand side may be called the survival probability (up to time τ ) of the transient random walk n (t) . See Corollary 5.6 for the general statement which connects joint q-moments of the stationary process Q (t) with a multiparticle version of n (t) (τ ). We call this multiparticle process the transient stochastic q-Boson system. Taking the long time limit of (1.2), we see that The right-hand side is the asymptotic survival probability that the transient random walk eventually runs off to infinity and is not absorbed at zero. This probability, viewed as a function of the initial location n, is a harmonic function 3 for the transient random walk n (t) (τ ), which, moreover, takes value 1 at n = +∞. This in fact identifies the harmonic function uniquely. From the stationarity of M qT q;t under the process Q (t) one can check that (1.1) satisfies the same harmonicity condition, which implies that (1.3) equals (1.1).
A general multiparticle argument involving identifying the "correct" harmonic function (asymptotic survival probability) of the transient stochastic q-Boson system with the joint q-moments of the q-TASEP leads to the proof that Q (t) converges to its stationary distribution M qT q;t (Theorem 5.7).
1.5. Outline. In Section 2 we give general definitions related to parameter-symmetric particle systems and swap operators. In Section 3 for q-Hahn TASEP we present an explicit realization of the parameter transposition in terms of a Markov swap operator corresponding to a random jump of a single particle. In Section 4 we pass to the continuous time in the q-Hahn TASEP, and define the q-Hahn backward process. This also implies the results about the TASEP from [PS19]. In Section 5 we define and study the dynamics preserving the q-TASEP distribution, and show its convergence to stationarity. In Section 6 we obtain swap operators for the beta polymer.
1.6. Acknowledgments. I am grateful Vadim Gorin for helpful discussions. I am grateful to the organizers of the workshop "Dimers, Ising Model, and their Interactions" and the support of the Banff International Research Station where a part of this work was done. The work was partially supported by the NSF grant DMS-1664617.

From symmetry to swap operators
This section contains an abstract discussion of stochastic particle systems on Z which depend symmetrically on their parameters. The main notions which we use in other parts of the paper are parameter-symmetric stochastic particle system and swap operators.
2.1. Parameter-symmetric particle systems. Let Conf f in (Z) be the space of particle configurations x = (. . . < x 3 < x 2 < x 1 ), x i ∈ Z, which can be obtained from the step configuration step := (. . . , −3, −2, −1) by finitely many operations of moving a particle to the right by one into the nearby empty spot. The space Conf f in (Z) is countable. By a multiparameter interacting particle system x(t) we mean a Markov process on Conf f in (Z) evolving in continuous or discrete time, such that x(0) = step. Assume that this Markov process depends on countably many parameters ν = {ν i } i∈Z ≥1 , ν i ∈ R d . 4 One should think that ν i is attached to the particle x i , but the distribution of each x j (t) may depend on all of the ν i 's. We denote the process depending on ν by x ν (t). In this section we assume that all the parameters ν i are pairwise distinct.
By imposing a specific distributional symmetry of x ν (t) under the action of S(∞) on the parameters ν, we arrive at the following definition: Definition 2.1. A multiparameter particle system x ν (t) is called parameter-symmetric, if for all n and t we have the following equality of joint distributions: x snν n+1 (t), x snν n−1 (t), . . . , x snν 1 (t)). (2.1) That is, the transposition s n preserves the joint distribution of all particles except x n .
Here is a straightforward corollary of this definition: Corollary 2.2. In a parameter-symmetric particle system, for any t and any σ ∈ S(n) × S n (∞), the random variables x ν n (t) and x σν n (t) have the same distribution. Remark 2.3. In Section 6 below we consider the beta polymer model, which may also be viewed as a particle system, but the particles live in (0, 1]. For concreteness, in the general discussion in this section we stick to particle systems in Z.

2.2.
Coupling. Let m 1 , m 2 be two probability measures on the same measurable space (E, F). A coupling between m 1 and m 2 is, by definition, a measure M = M (dz, dz ) on (E × E, F ⊗ F) whose marginals are m 1 (dz) and m 2 (dz ), respectively: A coupling is not defined uniquely, but always exists (the product measure M = m 1 ⊗ m 2 is an example).
In the notation of the previous subsection, start from a parameter-symmetric particle system x ν (t). Fix time t and index n ∈ Z ≥1 , and consider two distributions x ν (t) and x snν (t) on the same countable space Conf f in (Z). We would like to find a coupling M = M n between the distributions of x ν (t) and x snν (t) which satisfies an additional constraint corresponding to (2.1): Such a coupling also might not be defined uniquely. An example of a coupling satisfying (2.2) can be obtained by adapting the basic product measure example. Denote x ν n (t) := {x ν k (t) : k = n}, and similarly for x snν n (t). Define 3) Here δ(·) is the Dirac delta, P (·) is a shorthand notation for the probability distribution of x ν n (t), and the two quantities P (· | ·) are the conditional distributions of x ν n (t) (resp. x snν n (t)) given the locations of all other particles. Note that both conditional distributions P (· | ·) in (2.3) are supported on the same interval (if n = 1, then, by agreement, x 0 ≡ +∞ and the interval is infinite; for n ≥ 2 the interval is finite). The next statement follows from the above definitions: Lemma 2.4. The distribution M ind n defined by (2.3) is a coupling between the distributions of x ν (t) and x snν (t), and satisfies (2.2).

Swap operators.
With a coupling one can typically associate two conditional distributions. In our situation, a coupling M n satisfying (2.2) leads to the two distributions on I n (2.4) which we denote by p n (x snν n (t) | x ν (t)) and p n (x ν n (t) | x snν (t)). Indeed, under, say, p n it suffices to specify only the conditional distribution of x snν n (t), as all the other locations in x snν (t) stay the same. Thus, a coupling M n satisfying (2.2) is determined by either p n or p n .
In the particular example M ind n (2.3), the distribution p n simply corresponds to forgetting the previous location of x ν n (t), and selecting independently the new particle x snν n (t) ∈ I n (according to the distribution with the parameters s n ν) given the remaining configuration x snν n (t) = x ν n (t). This distribution p n corresponding to M ind n can be quite complicated as it may depend on the whole remaining configuration x ν n (t). This dependence may also nontrivially incorporate the time parameter t.
In this paper we describe specific parameter-symmetric particle systems for which there exist much simpler conditional probabilities p n or p n . Let us give a definition clarifying what we mean here by "simpler": Definition 2.5. The conditional probability p n is said to be local if p n (x snν n (t) | x ν (t)) depends only on n, ν, and three particle locations x ν n+1 (t), x ν n (t), x ν n−1 (t). The definition for p n is analogous.
We will interpret the local conditional probability p n as a Markov operator. When applied, p n leads to a random move x ν n (t) → x snν n (t) given x ν n+1 (t), x ν n−1 (t). In distribution the application of p n is equivalent to the swapping of the parameters ν n ↔ ν n+1 . Due to this interpretation, we will call p n the (Markov ) swap operator.
In the examples we consider, swap operators will also be independent of t.
Remark 2.6. Typically, only one of the probabilities p n and p n can be local (and thus correspond to a swap operator). Indeed, assuming that p n is local, we can write where x ν n (t) = x snν n (t), and we also assume that the probability in the denominator is nonzero. If one wants p n to be local, too, it is necessary that the ratio of the probabilities P (x ν (t)) P (x snν (t)) (in which the configurations x ν (t) and x snν (t) differ only by the location of the n-th particle) depends only on the four particle locations . This (quite strong) condition on the ratio of the probabilities does not hold for the particle systems considered in the present paper. (In particular, using the explicit Rakos-Schutz formula [RS06] expressing transition probabilities in TASEP with particle-dependent speeds as determinants one can check that the condition fails for the usual TASEP.)

Swap operators for q-Hahn TASEP
In this section we examine the parameter symmetry and swap operators for the q-Hahn TASEP introduced and studied in [Pov13], [Cor14]. A multiparameter version of the process preserving its integrability is due to [BP18a].
3.1. The q-deformed beta-binomial distribution. We first recall the definition and properties of the q-deformed beta-binomial distribution ϕ q,µ,ν from [Pov13], [Cor14]. We use the standard notation for the q-Pochhammer symbol ( Everywhere throughout the paper we assume that the main parameter q is between 0 and 1.
5 These conditions do not exhaust the full range of (q, µ, ν) for which the weights are nonnegative. See, e.g., [BP18a, Section 6.6.1] for additional families of parameters leading to nonnegative weights.
Define the following difference operator: . Let a function f (n 1 , . . . , n m ) from Z m to C satisfy the following two-body boundary conditions for all n ∈ Z m such that for some i ∈ {1, . . . , m}, n i = n i+1 . (In (3.2), only the i-th and the (i + 1)-st components of n are changed.) Then we have Here [∇ µ,ν ] i is the operator (3.1) applied in the i-th variable.

3.2.
Multiparameter q-Hahn TASEP. Here we recall the particle-inhomogeneous version of the q-Hahn TASEP from [BP18a, Section 6.6]. Let be parameters. To make the system nontrivial, the ν i 's should be uniformly bounded away from 1.
An example of a one-step transition in the q-Hahn TASEP, together with the corresponding probabilities for each particle. Here The q-Hahn TASEP starts from step and evolves in Conf f in (Z) in discrete time t ∈ Z ≥0 . At each time moment, each particle x i independently jumps to the right by j with probability where x 0 = +∞, by agreement. See Figure 2 for an illustration. For the step initial configuration the q-moments of the q-Hahn TASEP were obtained in [BP18a, Corollary 10.4] (in the homogeneous case ν i ≡ ν a proof using duality and coordinate Bethe ansatz is due to [Cor14]). The q-moments are given in the next proposition.
Proposition 3.4. For any ∈ Z ≥1 and any n 1 ≥ n 2 ≥ . . . ≥ n ≥ 1 with the assumption that we have for the q-moments of the q-Hahn TASEP started from step: (3.5) Here the integration contours are positively oriented simple closed curves which are q-nested around {ν j } j=1,...,n 1 (that is, each contour encircles the ν j 's, and, moreover, the z A contour encircles each z B contour, B > A) and leave 0 and 1 outside. See Figure 3 for an illustration. Figure 3. Possible integration contours in (3.5) for = 3. The contours for qz 3 , q 2 z 3 , and qz 2 are shown dotted.
Remark 3.5. Together with particle-dependent inhomogeneity governed by the parameters ν i , one can make the parameter γ time-dependent. That is, at each time step t − 1 → t, the jumping distribution (3.3) can be replaced by ϕ q,γtν i ,ν i (j | x i−1 − x i − 1). The moment formula (3.5) continues to hold when modified by replacing the term 1−γz i The main result of this section (Theorem 3.8 below) also holds in this generality, but for simplicity we will continue to assume that γ does not depend on t.
Since 0 < q < 1 and we start from step, the random variables j=1 q xn j (t)+n j are all between 0 and 1. Because the moment problem for bounded random variables admits a unique solution, the q-moments (3.5) uniquely determine the joint distribution of all the q-Hahn TASEP particles {x i (t)} i∈Z ≥1 at each fixed time moment. This implies the following statement: Proposition 3.6. The multiparameter q-Hahn TASEP with step initial configuration is a parametersymmetric particle system in the sense of Definition 2.1. Moreover, the distribution of each x n (t) depends on the parameters ν 1 , . . . , ν n in a symmetric way.
Denote the right-hand of (3.5) by f (n 1 , . . . , n ), where now n i ∈ Z are not necessarily ordered. Notice that if n = 0, the integrand has no poles inside the z (i.e., the smallest) integration contour. Therefore, f (n 1 , . . . , n −1 , 0) = 0. The following lemma will be employed in the next subsection.
Lemma 3.7. The function f (n 1 , . . . , n ) on Z defined before the lemma satisfies the two-body boundary conditions (3.2), where for n i = n i+1 ∈ Z ≥1 one should take ν = ν n i .
Proof. This statement essentially appears in [Cor14], see also [BCS14]. Its proof is rather short so we reproduce it here. When n i = n i+1 (denote them both by n), the part of the integrand in (3.5) depending on z i , z i+1 contains The left-hand side of the boundary conditions (3.2) for our function f and ν = ν n is as an integral over the contours as in Figure 3, where the integrand now contains where the important observation is that the denominator z i − qz i+1 has canceled out. Now the contour for z i can be deformed (without picking any residues) to coincide with the contour for z i+1 . However, thanks to the factor z i −z i+1 , the integrand is antisymmetric in z i , z i+1 . Therefore, the whole integral vanishes, as desired.
3.3. Markov swap operators for q-Hahn TASEP. Here we prove that the q-Hahn TASEP admits a local conditional distribution corresponding to the permutation s n = (n, n + 1) when the parameters satisfy ν n+1 < ν n before their swap. This leads to the Markov swap operator which we define now. Fix n ∈ Z ≥1 , and let (3.6) Observe that this probability does not depend on x n−1 . The condition ν n+1 < ν n ensures that the swap operator p qH n (3.6) has nonnegative probability weights. Theorem 3.8 (Theorem 1.1 in Introduction). Let x ν (t) be the q-Hahn TASEP with parameters ν = {ν i } i∈Z ≥1 , started from step. Fix n ∈ Z ≥1 and assume that ν n+1 < ν n . Replace x n (t) by a random x n (t) coming from the Markov swap operator p qH n (3.6). Then the new configuration is distributed as the q-Hahn TASEP x snν (t) with swapped parameters.
Proof. We will prove this theorem by applying p qH n in the q-moment formula. Since the q-moments uniquely determine the distribution, this computation will imply the claim. where m 1 ≥ . . . ≥ m > n + 1, n > m 1 ≥ . . . ≥ m ≥ 1, and k = + a + b + . Assume that the parameters ν i satisfy the contour existence condition (3.4). (In the end of the proof we will drop this assumption.) It suffices to show that is taken with the parameters before the swap), is equal to the expectation with the swapped parameters. Indeed, the sum over x n in (3.8) corresponds to the action of the swap operator p qH n on k j=1 q xn j (t)+n j (t) viewed as a function of {x i (t)}. We thus need to show that the expectation of the result with respect to the original parameters leads to the formula with the swapped parameters.
We now start from (3.8), and in the rest of the proof omit the dependence on t for shorter notation. First, we use the symmetry property (Lemma 3.2) to write for the part of the sum in (3.8) involving x n , x n+1 : (3.10) We thus need to compute where the vector n(r) = (n 1 (r), . . . , n k (r)) is as in (3.7), but with (a, b) replaced by (a + b − r, r). The expectation in (3.11) is given by the contour integral as in the right-hand side of (3.5).
Recall the notation f ( n) for this integral, where now n ∈ Z k , and the components of n are not necessarily ordered. By Lemma 3.7, this function f satisfies the two-body boundary conditions. Thus, the right-hand side of (3.11) can be rewritten by Lemma 3.3 as (recall notation (3.1) for the difference operator ∇ µ,ν ): Observe that now each of the difference operators [∇ ν n+1 νn ,ν n+1 ] +a+j can be applied independently inside the integral. We thus have for every variable w = z +a+j , j = 1, . . . , b: We see that the resulting integral coming from (3.8) contains, for each variable z +a+j corresponding to n +a+j = n in (3.7), the product over the parameters (ν 1 , . . . , ν n−1 , ν n+1 ) = s n (ν 1 , . . . , ν n ). Therefore, this integral is equal to the expectation (3.9) with the swapped parameters s n ν, as desired.
It remains to show that we can drop the contour existence assumption (3.4). The preceding argument implies that under (3.4), where P ν , P snν denote the q-Hahn probability distributions with the corresponding parameters at some fixed time t ∈ Z ≥0 . If n ≥ 2, the sum in the left-hand side of (3.12) is finite, and each probability P ν , P snν is a rational function of ν 2 , ν 3 , . . . (note that since . . . , x n+1 , x n , x n−1 , . . . , are fixed, only finitely many of the ν i 's enter (3.12)). The dependence on ν 1 is also rational after canceling out the common factor (γν 1 ;q) t ∞ (ν 1 ;q) t ∞ from both sides. Therefore, identity (3.12) between rational functions in ν i can be analytically continued, and the assumption (3.4) can be dropped.
For n = 1, the sum in the left-hand side of (3.12) becomes infinite. Remove the common factor (γν 1 ;q) t ∞ (ν 1 ;q) t ∞ from both sides again, then the coefficients by each power γ m , m ∈ Z ≥0 , become rational functions in ν i , i = 1, 2, . . .. Therefore, we can again analytically continue identity (3.12) and drop the assumption (3.4). This completes the proof.
Remark 3.9. When ν n = ν n+1 , we have from (3.6) that p qH n (x n | x n+1 , x n , x n−1 ) = 1 x n =xn (where 1 ··· stands for the indicator). Therefore, the swap operator reduces to the identity map, which is appropriate since for ν n = ν n+1 there is nothing to swap.
If ν n < ν n+1 , formula (3.6) for p qH n also makes sense, but some of these probability weights become negative. One can check that all algebraic manipulations in the proof of Theorem 3.8 are still valid for ν n < ν n+1 , but now they do not correspond to actual stochastic objects. This is the reason for the restriction ν n > ν n+1 in Theorem 3.8.
3.4. Duality for the q-Hahn swap operator. Here let us recall the Markov duality relation for the q-Hahn TASEP from [Cor14]. We will heavily use duality in Section 5 below.
Define the duality functional on the product space Conf f in (Z) × W as follows: (3.14) Let T qH(ν) (x, y), x, y ∈ Conf f in (Z), denote the one-step Markov transition operator of the q-Hahn TASEP with parameters ν = {ν i } and γ. We do not include the latter in the notation and assume that it is fixed throughout this subsection.
LetT qH(ν) ( n, m), n, m ∈ W , be the one-step transition operator of a discrete time Markov chain on W which at each time step evolves as follows. Independently at every site k ∈ Z ≥1 containing, say, y k particles, randomly select j particles with probability ϕ q,γν k ,ν k (j | y k ), and move them to the site k − 1. This Markov chain is called the (stochastic) q-Hahn Boson process. See Figure 4 for an illustration.
Here the operators T qH(ν) ,T qH(ν) act in the x and the n variables, respectively. Equivalently in terms of expectations, we have for all x 0 ∈ Conf f in (Z), n 0 ∈ W , and all times t ∈ Z ≥0 : Here in the left-hand side the expectation is taken with respect to the q-Hahn TASEP's evolution starting from x 0 , and in the right-hand side the expectation is with respect to the q-Hahn Boson process started from n 0 .
Consider now the Markov swap operator p qH k (3.6) on Conf f in (Z), where k ∈ Z ≥1 is fixed, and ν k > ν k+1 . It turns out that this operator admits a dual Markov operator on the space W , by means of the same duality functional H (3.14). Namely, definep qH k as the Markov operator which acts only on the k-th location in the q-Boson configuration. If the k-th location contains y k particles, thenp qH k randomly sends y k − j particles from location k to location k + 1, with probability ϕ q, ν k+1 ν k ,ν k+1 (j | y k ). See Figure 5 for an illustration. Proposition 3.11. If ν k > ν k+1 , then we have , where the operator in the left-hand side acts on x, and in the right-hand side -on n.
Proof. This duality relation immediately follows from computation (3.10) performed (with the help of Lemma 3.2) in the proof of Theorem 3.8.

Continuous time limit of repeated swaps
Here we consider two continuous time limits, one of the original q-Hahn TASEP, and another one of the transition probabilities p qH n leading to the new backward q-Hahn process. These two continuous time processes act on the two-parameter family {M qH q,ν;t } t∈R ≥0 , 0≤ν<1 of distributions of the continuous time q-Hahn TASEP (started from step) by continuously changing the parameters. 4.1. Two expansions of the distribution ϕ q,µ,ν . Let us write down two Taylor expansions of the q-deformed beta-binomial distribution from Section 3.1. Their proofs are straightforward.

4.2.
Continuous time q-Hahn TASEP. Fix q ∈ (0, 1), ν ∈ [0, 1) and consider the Poisson-type limit of the q-Hahn TASEP with homogeneous parameters ν i ≡ ν as where t ∈ Z ≥0 is the discrete time before the limit, and t ∈ R ≥0 is the scaled continuous time after the limit. The resulting process is the continuous time q-Hahn TASEP which evolves as follows. Starting from step, in continuous time t ∈ R ≥0 , each particle x n (t), n ∈ Z ≥1 , independently jumps to the right by j ∈ {1, 2, . . . , x n−1 (t) − x n (t) − 1} at rate (In continuous time there is at most one jump at every instance of time.) Here x 0 ≡ +∞, by agreement. This continuous time process was considered in [Tak14], [BC16b].
Remark 4.3 (Multiparameter continuous time q-Hahn TASEP). One can also consider a version of the continuous time q-Hahn TASEP in which each particle x n jumps at rate ν n ψ q,νn (· | ·). This multiparameter deformation preserves integrability. To get from the multiparameter process to the homogeneous one described above one has to set ν n ≡ ν and rescale the continuous time by the factor of ν. For simpler notation, we will mostly consider the homogeneous continuous time q-Hahn TASEP. Its multiparameter generalization is needed only in the proof of Theorem 4.7 below.
The continuous time q-Hahn TASEP possesses the same q-moment formulas as (3.5), with ν i ≡ ν, and the replacement inside the contour integral.

4.3.
Backward q-Hahn process. Here we define a continuous time limit of a certain combination of the swap operators p qH n (3.6). Let us first explain the main idea. Assume that ν 1 > ν 2 > ν 3 > . . .. By Theorem 3.8, the application of the Markov operators p qH 1 , p qH 2 , . . . (in this order) to a random configuration coming from the discrete time q-Hahn TASEP with parameters ν is equivalent in distribution to the permutation ν → . . . s 3 s 2 s 1 ν which exchanges ν 1 with ν 2 , then ν 1 with ν 3 , and so on (so that ν 1 gets pushed all the way to infinity and essentially disappears). Because the particle configuration to which we apply the p qH i 's is densely packed to the left, this application of the infinite product of the swap operators p qH i is well-defined and is a one-step Markov transition operator on Conf f in (Z). Its continuous time limit as ν i → ν for all i will be our new backward q-Hahn TASEP.
Let us now make this idea precise and take the particular parameters ν i = rν i−1 , i ∈ Z ≥0 , where ν ∈ [0, 1) and r ∈ (0, 1). Denote the infinite product of p qH 1 , p qH 2 , . . . as discussed above by B qH(r) q,ν . In detail, under the Markov operator B qH(r) q,ν , each particle x n jumps to the left into the new location x n ∈ {x n+1 + 1, x n+1 + 2, . . . , x n } chosen randomly from the distribution ϕ q,r n ,νr n (x n − x n+1 − 1 | x n − x n+1 − 1). (4.2) The update is sequential for n = 1, 2, . . .. In other words, the new location x n of each x n depends only on the two old locations x n+1 , x n .
Proposition 4.4. For ν i = νr i−1 and any m, k ∈ Z ≥0 we have q,r k−1 ν = δ step (T qH(r k ν) ) m , where T qH(ν) is the one-step Markov transition operator of the discrete time q-Hahn TASEP, δ step is the delta measure at the step configuration, and r k ν = {νr k+i−1 } i∈Z ≥1 is the parameter sequence shifted by k.
Note that T qH(r k ν) , the q-Hahn TASEP transition operator with the shifted parameter sequence r k ν from Proposition 4.4, is the same as T qH(ν) ν →r k ν , i.e., the original q-Hahn TASEP operator in which ν is replaced by r k ν.
We now aim to take the Poisson-type limit of the product B qH(r) where τ ∈ R ≥0 is the new scaled time. The parameters q, ν are assumed fixed. Observe that taking the limit ε 0 in the probabilities (4.2) and using Lemma 4.2 leads to an infinitesimal Markov generator under which each particle x n jumps to the left into x n ∈ {x n+1 +1, x n+1 +2, . . . , x n −1} at rate n · ψ • q,ν (x n − x n+1 − 1 | x n − x n+1 − 1). The factor n appears from the expansion r n = (1 − ε) n = 1 − nε + O(ε 2 ). The action of this generator is well-defined because the configuration from Conf f in (Z) is densely packed to the left, so only finitely many particles can jump in finite time. Denote this Markov generator by B qH q,ν . To be able to take a Poisson-type limit of the result of Proposition 4.4, we need a timeinhomogeneous Markov process. This is because the action of B qH(r) q,r k−1 ν changes the parameter ν as ν → r k ν. Therefore, let us define In words, B qH q,ν (τ, τ ) is the Markov transition operator from time τ to time τ of a process under which at each time s the jumps are governed by the infinitesimal generator B qH q,νe −s . We call the process corresponding to B qH q,ν (τ, τ ) the backward q-Hahn process. Proposition 4.5. In the regime (4.3) we have as Markov operators acting on the space Conf f in (Z).
Proof. Observe that the space of left-packed configurations Conf f in (Z) is countable, and under the Markov operators in both sides of (4.4) the particles jump only to the left. Therefore, the desired limit as ε 0 reduces to the limit of finite-size Markov transition matrices. For the latter the Poisson-type limit is taken in a straightforward way.
Remark 4.6. A time-and space-homogeneous version of the backward q-Hahn process was considered in [BC16b]. Indeed, one can check that the left jumps of the q-Hahn asymmetric exclusion process from [BC16b] coincide with the jumps at rates ψ • q,ν (x n −x n+1 −1 | x n −x n+1 −1).
However, the spatial inhomogeneity of the backward q-Hahn process does not allow to immediately apply the contour integral q-moment formulas from [BC16b] in our situation.

Action on distributions.
For t ∈ R ≥0 and ν ∈ [0, 1) denote by M qH q,ν;t the time t distribution of the continuous q-Hahn TASEP started from step. The combined results of Sections 4.2 and 4.3 imply the following action of the backward q-Hahn process on these distributions: Theorem 4.7. We have for all ν ∈ [0, 1) and t, τ ∈ R ≥0 : M qH q,ν;t B qH q,ν (0, τ ) = M qH q,νe −τ ;te −τ . In words, the time-inhomogeneous backward q-Hahn process maps the distribution M qH q,ν;t onto the distribution from the same family, but with rescaled parameters t and ν. See Figure 6 for an illustration of the action on parameters.
Proof. Take the continuous time multiparameter q-Hahn TASEP from Remark 4.3 and set the parameters to ν n = νr n−1 , where ν ∈ [0, 1) and r ∈ (0, 1). After rescaling the continuous time by ν, under this process each particle x n jumps by j at rate r n−1 ψ q,νr n−1 (j | x n−1 − x n − 1). Denote the distribution of this inhomogeneous process at time t ∈ R ≥0 started from step by M qH(r) q,ν;t . Clearly, lim r→1 M qH(r) q,ν;t = M qH q,ν;t .
A suitable modification of Proposition 4.4 applies to this r-dependent distribution M qH(r) q,ν;t , when we take a sequence of discrete Markov backward steps: Indeed, the application of all the operators B qH(r) q,r i ν , i = 0, 1, . . . , k − 1, turns r n−1 ψ q,νr n−1 , the jump rate of x n , into r k+n−1 ψ q,νr k+n−1 . The latter is the same as the old jump rate but with the parameters (ν, t) multiplied by r k . Taking the limit as r → 1 in (4.5) and using Proposition 4.5 implies the result. 4.5. Corollary. Mapping TASEP back in time. By setting q = ν = 0 in Theorem 4.7, we recover the main result of the recent paper [PS19] on the existence of a Markov process which maps the TASEP distributions back in time. Indeed, we have for the rates (4.1): where 1 ≤ j ≤ m and 0 ≤ j ≤ m − 1. Moreover, for q = ν = 0 the backward continuous time process B qH 0,0 (τ, τ ) = B qH 0,0 (τ − τ ) is time-homogeneous. Under this process, each particle x n independently jumps to the left into one of the holes {x n+1 + 1, . . . , x n − 1} at rate n per each hole. This dynamics is called the backward Hammersley process in [PS19].
We see from (4.6) that M qH 0,0;t is the distribution of the usual TASEP at time t started from the step initial configuration. Under the action of the backward Hammersley process B qH 0,0 for time τ , the distribution M qH 0,0;t maps into M qH 0,0;te −τ . This corollary of Theorem 4.7 is precisely Theorem 1 from [PS19]. In the latter paper the result was obtained in a completely different way using a well-known connection (e.g., see [BF14]) of TASEP and Schur processes, which are probability distributions on two-dimensional arrays of interlacing particles. In contrast, here we proved the more general Theorem 4.7 involving only observables of the particle systems in one space dimension, and did not rely on Schur like processes in two space dimensions.

Stationary dynamics on the q-TASEP distribution
When ν = 0, the q-Hahn TASEP turns into the q-TASEP introduced in [BC14]. We continue working in the continuous time setting as in Section 4. In this section we introduce and study a time-homogeneous, continuous time Markov process which is stationary on the distribution of the q-TASEP.
5.1. q-TASEP and the backward process. The q-TASEP is a continuous time Markov dynamics on Conf f in (Z) depending on a single parameter q. Under it, each particle x n (t) jumps to the right by one at rate 1 − q x n−1 (t)−xn(t)−1 (by agreement, x 0 ≡ +∞, so the first particle performs the Poisson random walk). When the destination of the jump is occupied, the rate is 1 − q 0 = 0, so the jump is blocked automatically. Denote the infinitesimal Markov generator of the q-TASEP by T. For the q-TASEP started from the step initial configuration step, let M qT q;t denote its distribution at time t ∈ R ≥0 . See Figure 7 (jumps to the right) for an illustration. Setting ν = 0 in the backward q-Hahn process from the previous Section 4, we arrive at the backward q-TASEP. This specialization in particular makes the backward process time-homogeneous. Under the backward q-TASEP process (we denote its continuous time by τ ), each particle x n (τ ), n ∈ Z ≥1 , jumps to the left into x n ∈ {x n+1 (τ ) + 1, x n+1 (τ ) + 2, . . . , x n (τ ) − 1} at rate (recall notation (4.1)) Remark 5.1. One can check that but we will not use this fact.
Denote the infinitesimal Markov generator of the backward q-TASEP by B. See Figure 7 (jumps to the left) for an illustration.

5.2.
Definition of the stationary dynamics. Fix the q-TASEP time parameter t ∈ R >0 . Introduce the notation This is an infinitesimal Markov generator of a process under which the particles move to the right according to the q-TASEP, and independently move to the left according to the backward q-TASEP slowed down by the factor of t. By slightly abusing notation, we denote the continuous time Markov process with the generator Q (t) by the same letter. We also adopt a convention of using the letter τ ∈ R ≥0 for the continuous time in the process Q (t) . Thus, t in Q (t) is a fixed parameter which enters the definition of the process.
Proposition 5.2. The continuous time Markov process with the generator Q (t) is well-defined and can start from any configuration x(0) ∈ Conf f in (Z).
Proof. The only problem in the definition of Q (t) is that it may have infinitely many jumps in finite time. First, observe that the q-TASEP is well-defined starting from any configuration x(0) ∈ Conf f in (Z). Next, under Q (t) particles go to the right not faster than under the q-TASEP. Therefore, with high probability, the random configuration of particles under Q (t) is empty to the right and densely packed to the left outside a bounded region of Z (the size of the region may depend on the time τ in Q (t) ). Because of this, the total jump rate of all particles under Q (t) is bounded. Therefore, Q (t) does not generate infinitely many jumps in finite time when started from any configuration x(0) ∈ Conf f in (Z). This completes the proof.
Proposition 5.3. The process Q (t) preserves the q-TASEP distribution M qT q;t .
Proof. By Theorem 4.7, the backward process B, ran for small time δ > 0, maps M qT q;t to M qT q;te −δ . After evolving this distribution under the q-TASEP for time t − te −δ > 0, we get back the original distribution M qT q;t . Differentiating in δ and sending δ 0, the infinitesimal Markov generator of the combined dynamics is readily seen to be t T + B = t Q (t) . As the factor t by Q (t) simply corresponds to the time scale (of the time variable τ ), we get the desired statement that the dynamics Q (t) preserves the measure M qT q;t .

Dual process: transient q-Boson.
Our aim now is to describe the dual process to Q (t) . by means of the same duality functional H(x, n) = 1 n >0 j=1 q xn j +n j (3.14). Recall that n ∈ W (3.13), and we interpret elements of W as -particle configurations in Z ≥0 . First consider the individual components T and B in Q (t) . The dual process to the q-TASEP is known as the stochastic q-Boson particle system [SW98] (see also [BCS14]). Under this process, particles move in continuous time from site k to k − 1, k ∈ Z ≥1 , independently at different sites. More precisely, if a site k ∈ Z ≥1 contains y k particles, then one particle hops from site k to site k − 1 at rate 1 − q y k . See Figure 8 (jumps to the left) for an illustration. Denote the infinitesimal Markov generator of the stochastic q-Boson process byT. The duality between T andT holds in the same sense as in Section 3.4: for any x ∈ Conf f in (Z) and n ∈ W , where T acts on x, andT acts on n. Figure 8. Examples of rates of possible jumps of the stochastic q-BosonT (dashed jumps to the left) andB, the process dual to the backward q-TASEP (gray jumps to the right). Each left jump involves only one particle, while right jumps may involve any number of particles in the stack (in the figure, "1p" means that one particle leaves the given stack, and so on).
Let us now define another continuous time Markov process on W which will be dual to the backward q-TASEP. Under this new process, particles move in continuous time from site k to k + 1, k ∈ Z ≥1 , independently at different sites. More precisely, if a site k ∈ Z ≥0 contains y k particles, then the process sends y k − j particles to site k + 1, where j ∈ {0, 1, . . . , y k−1 }, at rate k · ψ • q,0 (j | y k ). (In particular, particles cannot leave site 0.) Denote the infinitesimal Markov generator of this dynamics byB. See Figure 8 (right jumps) for an illustration.
Proposition 5.5. The Markov generatorB is dual to the backward q-TASEP generator B: for any x ∈ Conf f in (Z) and n ∈ W , where B acts on x, andB acts on n. Consequently, the infinitesimal Markov generatorQ (t) :=T + 1 tB is dual to Q (t) , the generator of the stationary dynamics on the q-TASEP distribution.
Proof. The duality between B andB follows from Proposition 3.11 via the continuous time limit described in Section 4. Indeed, the backward q-Hahn process is a continuous time limit of the combination of the steps p qH k applied at each site of the lattice Z. The dual processB is the limit of the same type of the combination of the dual stepsp qH k , together with the degeneration ν = 0. This implies the duality of B andB. The claim about the duality of Q (t) andQ (t) follows by linearity.
Because the rates of the right jumps underQ (t) grow as the particles of n get farther to the right, the processQ (t) is transient (except the absorption at n = 0), see the proof of Lemma 5.9 below for details. For this reason we callQ (t) the transient stochastic q-Boson particle system on W (transient q-Boson, for short). Let us record the duality between it and Q (t) in terms of expectations: Corollary 5.6. Fix t ∈ R >0 and ∈ Z ≥1 . Take any x 0 ∈ Conf f in (Z), and any n 0 ∈ W . Let {x(τ )} τ ∈R ≥0 be the process on Conf f in (Z) with generator Q (t) started from x 0 , and { n(τ )} τ ∈R ≥0 be the process on W with generatorQ (t) started from n 0 . Then for any τ ∈ R ≥0 we have H(x 0 , n(τ )). (5. 2) The example of this duality statement (and further discussion towards the results of the next Section 5.4) in the simplest case = 1 may be found in Section 1.4 in Introduction.

5.4.
Convergence to the stationary distribution. In this subsection we use duality to prove the following result: Theorem 5.7. Fix t ∈ R >0 . For any initial configuration x(0) ∈ Conf f in (Z) the Markov process x(τ ) with the generator Q (t) converges, as τ → +∞, to the stationary distribution M qT q;t (in the sense of joint distributions of arbitrary finite subcollections of particles).
The proof of Theorem 5.7 occupies the rest of this subsection. Note that if m = 0, we automatically have S τ ( m) = 0 for all τ . because once n reaches 0, it can never become positive again. Therefore, the quantities S τ ( m) decrease in τ due to monotonicity in τ of the events {n (s) > 0 for all s ∈ [0, τ ]}. Because S τ ( m)'s are nonnegative, the limit exists.
Observe that H(step, m) = 1 m >0 . Therefore, for the stationary process x(τ ) started from step we have E stat(t) step j=1 q xm j (τ )+m j = S τ ( m) for all m ∈ W (in particular, if m = 0 then both sides are zero). For other initial conditions for Q (t) this identity does not hold for finite time τ , but it still holds asymptotically: Lemma 5.9. For any m ∈ W and any initial data x 0 ∈ Conf f in (Z) for the stationary process where S( m) is the asymptotic survival probability of the transient q-Boson.
Proof. If m = 0, then both sides of (5.3) are zero, so it suffices to assume that m > 0. By duality (Corollary 5.6), the left-hand side of (5.3) is equal to Here x 0 n are the particle coordinates under the initial data x 0 . First, we show that the Markov process n(τ ) conditioned to stay in the region {n ≥ 1} is transient. Observe that we can couple the -particle process n(τ ) restricted to this region with a single-particle process Y (τ ) on Z ≥1 with jump rates The coupling is such that Y (0) = n (0), and n (τ ) jumps left or right whenever Y (τ ) jumps left or right, respectively. Also n (τ ) might jump to the right more often than Y (τ ). This implies that under this coupling we have Y (τ ) ≤ n (τ ). The process Y (τ ) on Z ≥1 is a standard example of a transient Markov process: it eventually (with probability 1) reaches the part {A, A + 1, . . .} ⊂ Z ≥1 (where A depends only on and t) where the average drift to the right is bounded away from zero. With positive probability Y (τ ) then never comes back from {A, A + 1, . . .} to the neighborhood of zero, and thus is transient.
Using the transience, we can lower bound S( m) by the (positive) probability of the event that: (1) if there were any particles at site 1 at time τ = 0, then all these particles leave 1 by right jumps; (2) after that, n (τ ) never comes back to 1 (and hence all other particles of n(τ ) also never come back to 1). This implies that S( m) > 0. Now denote by S( m) the event that n (τ ) > 0 for all τ conditioned on n(0) = m. Thus, P(S( m)) = S( m) > 0. We have for the expectation in (5.4): The second summand goes to zero as τ → +∞, since inside the event S( m) c , we have n (s) = 0 for all s ∈ (s 0 , +∞) (where s 0 is random but finite). In the first summand, observe that conditioned on the asymptotic survival S( m), we have almost surely due to transience that n j (τ ) → +∞, τ → +∞, for all j = 1, . . . , . Because the initial configuration x 0 ∈ Conf f in (Z) is densely packed to the left, we thus almost surely have q x 0 n j (τ ) +n j (τ ) → 1 for all j as τ → +∞. Therefore, the expectation of the product in the first summand tends to 1, and we see that (5.4) is equal to S( m), as desired.
The asymptotic survival probabilities S( m) satisfy certain normalization at infinity: Proof. We can assume that m > 0, otherwise both expressions S(·) in the claim are zero. The desired statement follows from the transience of the process n(τ ) (with generatorQ (t) ) as in the proof of the previous Lemma 5.9. Namely, we can lower bound the jump rate of n 1 (τ ) to the right from a site k ∈ Z ≥1 by const · k. Let n(τ ) start from m. If m 1 > R is large, the probability that the first particle n 1 (τ ) ever returns to the R/2-neighborhood of zero is close to zero. Thus, the probability S( m) that the process n(τ ) started from m survives and runs off to infinity is close to the asymptotic survival probability of the process on W −1 with one less particle and started from (m 2 , . . . , m ). In the special case = 1, the claim reads S(m 1 ) → 1 as m 1 → +∞, which clearly holds. This implies the claim. q xm j (t)+m j (5.5) for all and m ∈ W . Here the quantity in the left-hand side is the long time limit of the q-moment of Q (t) , and the right-hand side is the q-moment of the q-TASEP. This suffices since in our situation the q-moments uniquely characterize the distribution. We will establish (5.5) by showing that both sides satisfy the same equations (harmonicity with respect to the transient q-Boson) plus normalization at infinity which uniquely determine the function.
Lemma 5.12. As function of m, the survival probabilities S( m) are harmonic for the transient q-Boson, that is, Here the first identity is simply the expression for the action of the generator on a function, and the claim is that this action gives identical 0.
Proof. The argument is rather standard. Consider the evolution of the process n(·) started from m during short time dτ . Conditioned that the process stepped into m (which happens with probabilityQ (t) ( m, m )dτ ), the survival probability is then S( m ). With the complementary probability 1 − m : m = mQ (t) ( m, m )dτ = 1 +Q (t) ( m, m)dτ , the process did not leave m, and the survival probability did not change. Therefore, Taking the coefficient by dτ leads to the desired identity.
Lemma 5.13. Harmonicity condition (5.6) together with the normalization at infinity of Lemma 5.10 and the condition that S( m) = 0 whenever m = 0 uniquely determine the function S( m), m ∈ W .
Proof. Assume that G ( m), where = 1, 2, . . . and m ∈ W , is a family of harmonic functions satisfying normalization at infinity as in Lemma 5.10, that is, for any ε > 0 there exists R such that for all m ∈ W with m 1 > R we have (5.7) Moreover, we assume that G ( m) = 0 whenever m = 0. We will argue by induction on and show that G ( m) is equal to S( m), the asymptotic survival probability of the transient q-Boson on W . For = 1, it is straightforward to see that the space of harmonic functions vanishing at 0 is one-dimensional. In this case the normalization at infinity is the single condition G 1 (m 1 ) → 1 as m 1 → +∞, which determines the harmonic function uniquely.
Next, observe that because the function G is harmonic, it satisfies the following averaging property: where τ is arbitrary. Assume that m > 0 (this does not restrict the generality). For R ∈ Z ≥1 , take the stopping time T R := inf{τ ≥ 0 : n (τ ) = 0 or n 1 (τ ) = R}, where the process n(τ ) starts from m. This stopping time is almost surely bounded and has finite expectation. Then . . , n (T R )) 1 n 1 (T R )=R . (5.8) The first equality above follows from the optional stopping theorem, and the second one is the splitting into two cases, n 1 (T R ) = R or n (T R ) = 0. In the latter case the function G vanishes by our assumptions. Take any ε > 0 and choose R such that (5.7) holds. Thus, by the normalization at infinity, we have where we have replaced G −1 by S using the induction hypothesis. Since S satisfies the same conditions (harmonicity, normalization at infinity, and vanishing when m = 0) as the family G , identity (5.8) is also valid for S. Thus, Lemma 5.14. The q-moments of the q-TASEP in the right-hand side of (5.5) satisfy all the conditions listed in the previous Lemma 5. 13.
In particular, this implies the existence of a harmonic function satisfying the normalization at infinity and vanishing whenever the last component of a vector is zero.
Proof of Lemma 5.14. The harmonicity of the q-moments follows from the fact that the q-TASEP distribution M qT q;t is stationary under Q (t) (Proposition 5.3), together with duality between Q (t) and the transient q-Boson (Corollary 5.6).
The normalization at infinity follows from the fact that the q-TASEP started from step lives on Conf f in (Z), so for any t ∈ R ≥0 we almost surely have x m (t) + m → 0 as m → +∞.
The fact that the q-moments vanish when n = 0 follows from the agreement that x 0 = +∞.
Combining all the lemmas in this subsection, we get the desired Theorem 5.7.

Beta polymer
The q-Hahn TASEP has a remarkable degeneration -the beta polymer model introduced in [BC16a]. This model is also related to a random walk in dynamic beta random environment, but here we will formulate everything only in terms of the polymer model. In this section we present Markov swap operators for the multiparameter beta polymer. The swap operators can be realized as certain additional layers in the strict-weak lattice on which the beta polymer is defined.
6.1. Multiparameter beta polymer and its joint moments. Take parameters γ > 0 and ν n > 0, t, n ∈ Z ≥1 , such that ν n − γ > 0 for all n, Let B t,n ∼ Beta(ν n − γ, γ) be independent beta distributed random variables. Here by the beta distribution Beta(α, β) we mean the one with the density The (inhomogeneous) beta polymer {Z(t, n)} t,n∈Z ≥1 is a collection of random variables satisfying the random recursion with the initial condition Z(0, n) = 1 for all n ∈ Z ≥1 . See Figure 9 for a graphical interpretation of the beta polymer using the strict-weak lattice. The homogeneous version of the beta polymer was studied in [BC16a] as a scaling limit q → 1 of the q-Hahn TASEP. The multiparameter generalization is obtained through the same limit from our model described in Section 3.
The random variables Z(t, n) are between 0 and 1. 6 Because of this, the joint distribution of the beta polymer random variables {Z(t, n)} n∈Z ≥1 (for every fixed t) is determined by the joint moments. These moments have the following form: The weight of a path is defined as the product of its edge weights. The initial condition along the left boundary is Z(0, n) = 1 for all n, so Z(t, n) = 1 above the diagonal (i.e., for n > t).
Remark 6.2. Similarly to Remark 3.5, one can generalize the beta polymer model so that the parameter γ depends on t. The moment formula (6.1) and our main result (Theorem 6.4 below) would continue to hold with straightforward modifications. For simplicity, everywhere below we take γ independent of t.
Corollary 6.3. The multiparameter beta polymer, viewed as a stochastic particle system Z(t, ·) on [0, 1] with time t ∈ Z ≥1 , is parameter-symmetric in the sense of Definition 2.1. Moreover, the distribution of each Z(t, n) depends on the parameters ν 1 , . . . , ν n in a symmetric way.
In other words, when ν n < ν n+1 , the beta polymer admits a Markov swap operator p beta n (in the sense of Definition 2.5) which acts by splitting the segment [Z ν (t, n), Z ν (t, n + 1)] ⊂ [0, 1] as 1 −B :B, and replacing Z ν (t, n) by the separating point.
Finally, each of the operators ∇ beta νn/ν n+1 can be applied separately under the contour integral in g(t; ·) given by (6.1), and we obtain (with the notation w = z +a+j , j = 1, . . . , b): Thus, we see that (6.4) is equal to the expectation E beta(snν) (Z(t, n 1 ) . . . Z(t, n k )) with the swapped parameters ν n ↔ ν n+1 , where (n 1 , . . . , n k ) is given by (6.3). Since joint moments determine the distribution of the beta polymer, we are done.
6.3. Polymer interpretation. Let us give a polymer interpretation of Theorem 6.4 (assuming that ν n < ν n+1 ). First, observe that the quantitỹ Z(t, n) :=B n Z(t, n + 1) + (1 −B n )Z(t, n), (6.6) whereB n ∼ Beta(ν n+1 − ν n , ν n ), is a beta polymer type partition function on a modified lattice. This modified lattice coincides with the one in Figure 9 in the vertical strip {0, 1, . . . , t} × Z ≥1 , has an additional vertex A, and two additional directed edges (t, n) → A and (t, n + 1) → A with weights 1 −B n andB n , respectively. The partition function from the line {0} × Z ≥1 to A is preciselyZ(t, n). See Figure 10 for an illustration. Let us now iterate the swapping in Theorem 6.4 and interchange the parameter ν 1 with ν 2 , then with ν 3 , and so on up to infinity. Assume that ν 1 < ν 2 < . . .. Let B (1) n ∼ Beta(ν n+1 − ν 1 , ν 1 ), n ∈ Z ≥1 , be independent random variables which are also independent of the environment {B t,n } in the beta polymer. Define Proposition 6.5. The joint distribution of {Z (1) (t, n)} n∈Z ≥1 defined above coincides with the joint distribution of the beta polymer {Z (ν 2 ,ν 3 ,...) (t, n)} n∈Z ≥1 at same t, but with the sequence of parameters shifted by one.
Define Z (s) (t, n) to be the polymer partition function from the line {0}×Z ≥1 to the point (s+t, n), s ∈ Z ≥1 , in the modified strict-weak lattice which coincides with the original lattice in Figure 9 in the vertical strip {0, 1, . . . , t} × Z ≥1 . To the right of this strip, the modified lattice is made out of down-right diagonal and horizontal edges with the weights B (s) n on each (t+s−1, n+1) → (t+s, n), and 1 − B (s) n on each (t + s − 1, n) → (t + s, n). See Figure 11 for an illustration.
Proposition 6.6. For any fixed s and t, the joint distribution of the partition functions {Z (s) (t, n)} n∈Z ≥1 on the modified lattice coincides with the joint distribution of the beta polymer partition functions {Z (ν s+1 ,ν s+2 ,...) (t, n)} n∈Z ≥1 with the same t, but with the parameter sequence ν shifted by s.
We will take the scaling limit of the beta polymer model as ν n = εν n , γ = εγ, whereν n >γ > 0 for all n, and 0 <ν 1 <ν 2 < . . .. The edge weights in the lattice in Figure 11 turn into the ones given in Figure 12.  where the directed paths π are as in Figure 11, and the edge weights are given in Figure 12. If s = 0, then we mean the unmodified first-passage time (as studied in [BC16a]). Above the main diagonal (i.e., for n > t) we have F (0) (t, n) = 0 because due to the presence of the Bernoulli components, there always exists a path with zero total weight between (t, n), n > t, and the vertical axis {0} × Z ≥1 . For the first-passage percolation model, an analogue of Proposition 6.6 holds: Proposition 6.8. For fixed s, t, the joint distribution of the first-passage times {F (s) (t, n)} n∈Z ≥1 with parametersν 1 <ν 2 < . . . in the modified lattice is the same as that of the unmodified ones {F (0) (t, n)} n∈Z ≥1 , but with the shifted sequence of parametersν s+1 <ν s+2 < . . ..