Multi-point distribution of discrete time periodic TASEP

We study the one-dimensional discrete time totally asymmetric simple exclusion process with parallel update rules on a spatially periodic domain. A multi-point space-time joint distribution formula is obtained for general initial conditions. The formula involves contour integrals of Fredholm determinants with kernels acting on certain discrete spaces. For a class of initial conditions satisfying certain technical assumptions, we are able to derive large-time, large-period limit of the joint distribution, under the relaxation time scale t=O(L3/2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=O(L^{3/2})$$\end{document} when the height fluctuations are critically affected by the finite geometry. The assumptions are verified for the step and flat initial conditions. As a corollary we obtain the multi-point distribution of discrete time TASEP on the whole integer lattice Z\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {Z}}$$\end{document} by taking the period L large enough so that the finite-time distribution is not affected by the boundary. The large time limit for the multi-time distribution of discrete time TASEP on Z\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {Z}}$$\end{document} is then obtained for the step initial condition.


Introduction
Models in the one-dimensional Kardar-Parisi-Zhang (KPZ) universality class are expected to have the same limiting height fluctuations under the t : t 2/3 : t 1/3 KPZ scaling for temporal correlations, spatial correlations and height fluctuations. Describing the universal limiting fluctuating field and proving the convergence of concrete models to the universal limiting field have been the main goals in the field and attracted active research over the last two decades. The convergence to the limiting fluctuating field (at one-point or multi-point level) have been obtained for a large class of models (mostly exactly solvable), either for the whole space models [2,[8][9][10][11]13,15,16,20,21], or for half space models [1,6,7].
More recently similar limit theorems were obtained for models on periodic domains. Such models are defined on a finite ring of size L instead of on the infinite or semiinfinite lattices. The crucial interest here is to understand the effect of this finite geometry on the height fluctuations. For periodic models the one-point marginals of limiting height fluctuations were obtained for classical step, flat and stationary initial conditions in [3,14], and independently in [19] under the so-called relaxation time scale t = O(L 3/2 ) when the height fluctuations are critically affected by the finite geometry. Further extensions including formulas for multi-point space-time joint distributions were subsequently obtained in [4] for step initial condition and [5] for more general initial conditions. These limiting formulas are expected to be universal for all related models on periodic domains. So far all limiting distribution formulas were obtained as scaling limits of finite-time joint distributions of a single model in the KPZ universality class, namely the continuous time periodic totally asymmetric simple exclusion process (TASEP on a ring).
The goal of this paper is to study another classical model in the KPZ universality class, the discrete time totally asymmetric exclusion process with parallel updates, on periodic domains. On the infinite lattice Z, this model has been well studied. The one-point marginal distribution was obtained in [11] for the equivalent geometric last passage percolation model and joint distributions of several locations at equal time was obtained in [9]. Recently the joint distributions along the time direction have also been studied, in [12] for two-time case and [13] for general multi-time joint distributions. However on the periodic domain there are fewer results concerning height fluctuations (see [18] for some results on transition probability and stationary distributions), which are the main focuses of this paper. The main results of this paper are summarized as follows: (1) For general initial conditions we obtain a finite-time multi-point (in both space and time) joint distribution formula for discrete time parallel periodic TASEP. The formula consists of an m-fold (m is the number of space-time points being considered in the joint distribution) contour integral with integrand involving a Fredholm determinant, where the Fredholm determinant has kernels acting on certain discrete sets related to the roots of some polynomials. (2) Under the relaxation time scale t = O(L 3/2 ) and the 1 : 2 : 3 KPZ scaling, we obtain large-time, large-period limits for the multi-point joint distributions under certain assumptions on the initial condition, which are verified for the step and flat cases. These limiting formulas agree with those obtained in [4,5], thus providing an evidence that the height fluctuations for periodic models in the KPZ class are in fact universal. (3) We also obtain the same type of multi-point joint distribution formula for discrete time TASEP on the infinite lattice Z, by relating it with the corresponding periodic model with large enough period so that the finite-time distribution is not affected by the boundary. The large time limit for the multi-time distribution formula under KPZ scaling is then derived for step initial condition, which agrees with the one in [15].
Comparing to the previous work [4,5,15] on the multi-point distributions of continuous time TASEP, on either periodic domain or Z, our work consists of formulas with similar structures but involves an extra parameter p describing the hopping probability, which makes the algebraic properties a bit more complicated. In particular the polynomial whose roots are related to the kernels in the periodic formulas now depends on the extra parameter p. We point out that all the formulas for continuous time TASEP can be obtained from our formulas by rescaling the time and taking p → 0.
On the infinite lattice Z, the multi-time distribution of discrete time TASEP (more precisely the equivalent geometric last passage percolation model) has been considered recently in [13], where a formula involving contour integrals of Fredholm determinants was obtained. We expect our formulas to agree with the formulas in [13], both in finite time and for the large time limit. Currently we are not able to prove the equivalence of the two types of formulas though they share a lot of similarities.

Outline of the paper
In Sect. 2 we describe the discrete time parallel TASEP models and state the main results involving several multi-point joint distribution formulas, for both periodic domain and infinite lattice Z, finite time and large time limit. In Sects. 3, 4 and 5 we derive the main finite-time algebraic formulas for multi-point distribution of discrete time parallel periodic TASEP and we regard them as the main technical novelties in this paper. In particular in Sect. 3 we prove a novel transition probability formula for discrete time parallel TASEP on the periodic domain. The formula involves integral of determinants and is proved using coordinate Bethe ansatz. In Sect. 4 we derive the finite-time multi-point distribution formula by performing a multiple sum over the transition probabilities. The key ingredients are certain Cauchy-type summation identities over the eigenfunctions of the generator, which might be of independent interests so we discuss the proof in Sect. 5. In Sects. 6 and 7 we discuss the large time, large period asymptotics for the multi-point distribution under the relaxation time scale t = O(L 3/2 ). Finally in Sect. 8 we derive the multi-time distribution for infinite discrete time parallel TASEP and perform large time asymptotics under KPZ scaling for the step initial condition.

Models and main results
Let N < L be positive integers. We consider discrete time TASEP with parallel updates with N particles on a spatially periodic domain of size L. It is convenient to view the dynamics as particles moving to the right on the integer lattice Z while periodicity forces particle configurations to be identical copies of each other every L sites. More precisely this means that the occupation functions η j (t) which equals 1 if there is a particle at site j ∈ Z at time t and equals 0 otherwise should satisfy η j (t) = η j+k L (t), for all j, k ∈ Z and t ∈ N. (2.1) We fix a single period of size L and denote the locations of totally N particles in this period at time t as Here x i (t) ∈ Z for 1 ≤ i ≤ N and we index the particles from right to left. The locations of all the particles then satisfy x i+k N (t) = x i (t) − k L for 1 ≤ i ≤ N and k ∈ Z so that we have Thus the natural configuration space for the particles should be X (L) The discrete time parallel periodic TASEP with N particles, period L and hopping probability 0 < p < 1, which we denote by dpTASEP(L, N , p), is the following Markovian dynamics on particle configurations x(t) ∈ X (L) N : at each time step, each particle in a single period hops to its right neighbor site independently with probability p = 1 − q provided that the site is empty, otherwise it stays at its current position. As a special case, the discrete time parallel TASEP on Z which we will denote by dTASEP( p) corresponds to particles following the same evolution rules with the period L → ∞. In this case the configuration space for the first N particles (from right to left, we assume the existence of a right-most particle) becomes (2.3) Notation 2.1 Throughout the paper there will be several very similar formulas and quantities corresponding to either discrete time parallel TASEP on a periodic domain with size L or on the infinite lattice Z, as well as their scaling limits. To avoid confusion we will always follow the convention used in (2.2) and (2.3) by putting superscripts (L) to any quantities related to periodic model in finite-time and (∞) to any quantities related to model on the infinite lattice Z in finite-time. On the other hand quantities describing the large-time limits will have superscripts (per) and (kpz) for periodic and infinite models, respectively.
The main results of the paper consist of several formulas for multi-point joint distributions of discrete time parallel TASEP on periodic domain and on Z as well as their large time scaling limits. All of them are expressed as contour integrals of Fredholm determinants. The main difference is that the formulas for periodic models are expressed in terms of discrete sets related to roots of certain algebraic equations, while the Fredholm determinant appearing in formulas for discrete time TASEP on Z has kernels acting on continuous contours.

Multi-point distribution formulas for dpTASEP(L, N, p)
The first theorem is a finite-time multi-point joint distribution formula for discrete time TASEP on a periodic domain. This is the starting point of all other results in this paper. Theorem 2.2 (Finite-time multi-point joint distribution for dpTASEP(L, N , p, y)). Let N < L be integers. Consider discrete time periodic parallel TASEP with hopping probability 0 < p < 1, N particles and L sites per period (denoted by dpTASEP (L, N , p, y)). Here y ∈ X (L) N is the initial condition, i.e., x i (0) = y i for 1 ≤ i ≤ N . Set ρ = N /L. Fix a positive integer m and let (k i , t i ), 1 ≤ i ≤ m be m distinct points in Z × Z ≥0 satisfying 0 ≤ t 1 ≤ · · · ≤ t m . Then for any integers a 1 , . . . , a m , Hre we denote z := (z 1 , . . . , z m ). The contours are over nested circles centered at the origin: 0 < |z m | < · · · < |z 1 | < r c , where r c is a parameter depending on p and ρ that is defined in (2.16). Here and in all the remaining results we suppress the dependence of the integrand on the parameters a i , k i , t i as well as the hopping probability 0 < p < 1. The function C  [4] for the special step initial condition) on the finite-time multi-point distribution of continuous time periodic TASEP. In fact their formulas can be obtained from our formula (2.4) by taking p = ,t = t and letting → 0.
Next we state the theorem on the large-time, large-period scaling limit of (2.4) under the relaxation time scale t = O(L 3/2 ). To emphasize the dependence on the initial condition, we add the subscript "ic" for the terms in the limit which depend on the initial conditions. For convenience we make the following choice of labelling: we assume that x 1 (0) ≤ 0 < x 0 (0). This is equivalent with assuming that the initial condition satisfies y 1 ≤ 0 < y N + L.
Theorem 2.4 (Relaxation time limit). Consider a sequence of Markovian dynamics dpTASEP(L, N , y(L)) depending on the period L, where ρ = ρ L = N /L stays in a compact subset of (0, 1) and y 1 ≤ 0 < y N + L. Suppose that the sequence of initial conditions y(L) satisfies certain assumptions (see Assumption 6.2). Fix a positive integer m and let p j = (γ j , τ j ) be m points in the region Then for every fixed x 1 , . . . , x m ∈ R and parameters a i , k i , t i , 1 ≤ i ≤ m given by Here the constants c i depend explicitly on particle density 0 < ρ < 1 and hopping probability 0 < p < 1 and are given by where we set ν := √ 1 − 4 p · ρ(1 − ρ) for convenience. We recall that in (2.7) P (L) denotes the probability associated to dpTASEP(L, N , y(L)). The function F per ic (·) agrees with the one defined in Section 6.4 of [5] as the relaxation time limit of multipoint distribution of continuous time periodic TASEP. We recall the definition of F per ic (·) in Sect. 6.3 for completeness. The convergence is locally uniform in x j , τ j , and γ j .

Multi-point distribution formulas for dTASEP(p)
The finite-time multi-point joint distribution for discrete time parallel TASEP on Z can be obtained from Theorem 2.2 as a corollary. The rough idea is to take the period L sufficiently large so that the finite-time formula is unaffected by the boundary. The precise procedure is more delicate and will be explained in Sect. 8. For simplicity we will only state result on step initial condition but we point out that in [15] similar multitime formulas for continuous time TASEP was obtained for general initial conditions and we believe their arguments can also be adapted to our model. The main result is the following: Theorem 2.5 (Finite-time multi-point joint distribution for dTASEP( p)). Assume y = (−1, −2, . . . , −N ). Consider discrete time parallel TASEP on Z with initial particle locations x i (0) = y i for 1 ≤ i ≤ N . Let m ≥ 1 be a positive integer and (k 1 , t 1 ), . . . , (k m , t m ) be m distinct points in {1, . . . , N } × Z ≥0 . Assume that 0 ≤ t 1 ≤ · · · ≤ t m . Then, for any integers a 1 , . . . , a m , where the integral contours are circles centered at the origin of radii less than 1.

Remark 2.6
Under the standard coupling between discrete time parallel TASEP and geometric last passage percolation, we have the following equality in distributions: where G(m, n) is the point-to-point last passage time from (0, 0) to (m, n) in the usual geometric last passage percolation model with parameter q = 1 − p. For the right-hand side of (2.10) a formula of similar form as (2.9) was obtained in Theorem 2 of [13]. We expect the two formulas to be equivalent but currently we are not able to prove it for m ≥ 2.
Starting from the finite-time formula, we can then take large-time limit under the 1 : 2 : 3 KPZ scaling.

Theorem 2.7
Consider discrete time parallel TASEP on Z with step initial condition x i (0) = −i. We rescale the parameters such that (2.12) Here D kpz step (θ 1 , . . . , θ m−1 ) is a Fredholm determinant defined in Definition 8. 4. We remark that the right-hand side of (2.12) agrees with the right-hand side of equation (2.19) in [15], with the re-scaled parameters changed from Our parameters are chosen to be consistent with the periodic cases.

Bethe equations and Bethe roots
For our analysis on the discrete time periodic TASEP, the following polynomial and its roots play essential role: Definition 2.8 (Bethe roots). Given z ∈ C and 0 < p < 1. Define the degree L polynomial q z (w) by (2.13) We call this polynomial the Bethe polynomial associated to z and its roots the Bethe roots. We denote the set of all roots of the Bethe polynomial q z (w) by S z : 14) The Bethe root set S z is contained in the level set {w ∈ C : |w| ρ |1 + w| 1−ρ = |z|·|1+ pw| ρ }, which is sometimes called a deformed Cassini oval. For given jumping rate p and density ρ, there exists a critical radius |z| = r c for the level sets. For such critical z, the level set contains a self-intersecting point at w = w c . Here and It is not hard to check that for |z| < r c , the level set consists of two disjoint contours while for |z| > r c the two contours merge to a single contour. For z = r c there is a self-intersecting point for the contour at w = w c . See Fig. 1 for an illustration. We remark that the Bethe polynomials (and their roots) we are considering here are oneparameter generalizations of those considered in [3][4][5] which is related to continuous time periodic TASEP and corresponds to p = 0 degenerations of (2.13).

Definition 2.9 (Left and right Bethe roots).
For |z| < r c , we define the sets where w c and r c are defined in (2.15) and (2.16). Then it is straightforward to check that |L z | = L − N and |R z | = N . Roots in L z and R z are called left and right Bethe roots, respectively. We also define the left and right Bethe polynomials q z,L (w) and q z,R (w) as the monic polynomials with roots in L z and R z : Then by definition we have

A symmetric function related to initial conditions
In our finite-time multi-point distribution formula (2.4), the quantities encoding information in the initial condition are all related to the following symmetric function: Definition 2.10 (Symmetric function). Given p ∈ C and λ = (λ 1 , . . . , λ N ) ∈ Z N with λ 1 ≥ · · · ≥ λ N . We define

Remark 2.11
The symmetric function F λ defined here is a one-parameter generalization of the Grothendieck-like symmetric function defined in equation (3.6) of [5] which corresponds to the p = 0 degeneration in our situation.
The following two quantities related to F λ encode the initial condition: For |z| < r c , we define the global energy E y (z) by When E y (z) = 0, we define the characteristic function χ y (v, u; z) for a left Bethe root u and a right Bethe root v by

Remark 2.13
A straightforward calculation shows that for step initial condition y = (−1, . . . , −N ), we have F λ( y) ≡ 1. Hence the global energy function and characteristic function are both constant 1 for the step initial condition. In general since the roots of the Bethe polynomial q z (w) depend analytically on z, the function E y (z) is analytic for |z| < r c . Furthermore E y (z) can not be constant zero. In fact E y (0) = 1 since when |z| → 0 all the right Bethe roots converge to 0. As a consequence F λ( y) (R z ; p) is nonzero for all but finitely many z in any compact subset of {|z| < r c }, which means χ y (v, u; z) is a well-defined meromorphic function in z on {|z| < r c } for fixed u, v.

Definition of C (L) y ( z) and D (L) y ( z)
is a Fredholm determinant with kernel acting on certain 2 space over discrete sets related to the Bethe roots. More precisely for m distinct complex numbers z i satisfying |z i | < r c , we define the discrete sets and The functions J (w) and Q i ( j) are defined as follows: (2.34) Finally the functions f i (w) encodes the information of the parameters a i , k i , t i 's: for 1 ≤ ≤ m and F 0 (w) := 1.

Definition of D (∞)
The function D (∞) step (θ 1 , . . . , θ m−1 ) has a similar structure as D (L) y ( z) defined in Definition 2.15. It is also a Fredholm determinant of the form det(I − K The kernels also share certain similarities and we will illustrate the relationship between the two in Sect. 8.

Spaces of the operators
Instead of discrete sets, the operators appearing in D (∞) step (θ 1 , . . . , θ m−1 ) are defined on specific spaces of nested contours defined as follows: Let L and R be two simply connected regions on the complex plane such that Definition 2. 16 We define where two operators are defined by the kernels Here the function f i (w) is the same as the one appearing in the kernels for periodic case defined in (2.35) and we define if j = m is even or j = 1.

Transition probability
In this section we give an explicit integral formula for the transition probability of discrete time parallel TASEP in the configuration space X N . This is the starting point for deriving the finite-time joint distribution formulas.
Here is any simple closed contour with 0 inside and S z consists of all the roots of the degree L polynomial q z (w) The functions F i, j (w; t) and J (w) are given by

3)
and J (w) := w(w + 1)(1 + pw) Remark 3. 2 We remark that a different formula for the transition probability of discrete time parallel TASEP on a ring was obtained in [18]. The key difference is that the formula in [18] is expressed as an infinite sum of determinants while our formula is a single contour integral of determinants. The main reason for this is that in [18] the authors do not distinguish particle configurations which differ by a translation of an integer multiple of the period, so in our language their transition probability really is We believe our formula is simpler and more suitable for deriving finite-time joint distributions.

Remark 3.3
If we take the continuous time limit by setting p = , t = T / and send → 0, the dynamics then becomes continuous time periodic TASEP considered in [3] and our formula (3.1) reduces to equation (5.4) in [3] for the transition probability of continuous time periodic TASEP up to an index reversing.
We first list a few elementary properties of the transition probability formula (3.1) before discussing the proof of Proposition 3.1.

Proposition 3.4 (Properties of the transition probability formula)
The right hand side of (3.1) satisfies the following properties: (i) The right hand side of (3.1) can also be written as: has the same roots as q z (w) for any |z| > 0 and S z is any simple closed contour with all the roots in S z inside and −1 and −1/ p outside. (ii) The outer integral with respect to z in (3.1) (and also (3.5)) does not depend on the contour . (iii) Assume further that L ≥ x 1 − y N + 2. Then the right-hand side of (3.1) can be further written as

6)
where 0,−1 is any simple closed contour enclosing 0 and −1 as the only possible poles for the integrand. Note that (3.6) agrees with the transition probability for discrete time parallel TASEP on Z, see for example equation (3.21) of [9]. (iv) The right-hand side of (3.1) is invariant under cyclic translation and re-indexing.
. Hence by the residue theorem we have Note that F i, j (w, x, y, t)q z (w)+z L q z (w) is analytic at w = 0 for any 1 ≤ i, j ≤ N and the only possible poles of the integrand besides S z are w = −1 and w = −1/ p. (ii) Choose R > 0 large enough and > 0 small enough so that all the roots in S z are inside the region (3.7) Since R and can be arbitrarily large or small, the right hand side of (3.7) as a function in z is analytic for any |z| > 0. Hence the integral with respect to z in (3.1) is independent of since the integrand has an analytic continuation to Now for fixed R large enough and small enough, the right hand side of (3.8) is an analytic function in z for |z| sufficiently small such that all the roots ofq z are in the region {|w + 1| ≤ R}\{| pw + 1| ≤ }. Now by the residue theorem the outer contour integral with respect to z in (3.5) equals the integrand evaluated at z = 0, which equals .
(iv) We remark that this property can be easily understood if we use the probabilistic interpretation since x and x (and also y and y ) actually represent the same particle configuration on Z (we just use particles in different period as representatives), hence the transition probability between y and x and y and x should be the same. Here however we can not directly use this since we have not proven (3.1). In fact, we will need this fact in our proof of (3.1) so we give an independent algebraic proof here. It suffices to assume k = 2.
Here again Then the last line in the above equation equals We claim that for any fixedσ ∈ S N , This then implies To see the claim note that for i = 1 we have τ (i) = i − 1 and forσ For the other situations we split into two cases: Similarly if we denoteσ −1 (1) = j, then the term involvingw j among the sum in (3.9) equals so the claim follows and this completes the proof of the cyclic invariance.

Proof of the transition probability formula
Now we turn to the proof of formula (3.1). The proof basically follows the idea of [3,9] where we replace the Kolmogorov forward equation by a free evolution equation with additional boundary conditions. The main extra difficulty here is that due to parallel update rule, the stationary distribution for the dynamics is non-uniform. In fact one can satisfies a relatively simpler dynamics than P t ( y → x).
More precisely we have: Lemma 3.5 Given x, y ∈ Z N . Let G( x, t; y, 0) be the (unique) solution of the following free evolution equation (3.10) together with the boundary conditions: (3.11) and the initial condition: .
Note that we have used the convention x 0 = x N + L so when i = 1, (3.11) should be interpreted as (3.14) Proof Set . (3.15) It suffices to show that for x, y ∈ X (L) N , H ( x, t; y, 0) and G( x, t; y, 0) satisfy the same evolution equation. This then implies that H ( x, t; y, 0) = G( x, t; y, 0) for all x, y ∈ X (L) N since they satisfy the same initial condition. To better describe the evolution equation for H ( x, t; y, 0), it is convenient to introduce the notion of clusters of a particle configuration x.
Namely particle i through i − k + 1 are right next to each other while there are at least one empty site to the left of x i and right of x i−k+1 . Here we abuse notation by allowing the index to exceed {1, . . . , N } and this should be understood with the convention forms a cluster so that all the indices appearing will be between 1 and N .
Let N c ( x) be the number of clusters in configuration x and let x c j , 1 ≤ j ≤ N c ( x) be the locations of the left-most particles in each cluster. Then it is straightforward to check that H ( x, t; y, 0) satisfies where e c j ∈ Z N has 1 in the c j -th coordinate and 0 in the other coordinates.
We claim that for x, y ∈ X (L) N , (3.16) and (3.10) take the same form (with H replaced by G) provided that G( x, t; y, 0) satisfies boundary conditions (3.11). Due to the sum of products form of (3.16) and (3.10) it suffices to check of size m. We will show the stronger statement: (3.17) actually holds for any We do not require empty sites at the left and right ends so they may not form a cluster.
We prove this by induction on m. For m = 1 this is trivial. Assume the claim is true for any clusters of size ≤ m.
Here without loss of generality we can assume i − m ≥ 1, otherwise replace j by j + N for indices j ≤ 0. Then (3.18) where we used the induction hypothesis in the second equality of (3.18) for the sum inside the brackets. For notational conveniece we have suppressed the dependence on y in the last line to save space. Now by the boundary conditions (3.11) (possibly (3.14)) we have Inserting (3.19) into (3.18) we see the last line of (3.18) simplifies to which is precisely the left hand side of (3.17) and this completes the proof of Lemma 3.5.
Now we discuss the proof of Proposition 3.1.

Proof of Proposition 3.1 By Lemma 3.5 it suffices to prove
To see this, we check that the right-hand side of (3.21) satisfies the free evolution equation (3.10), the boundary conditions (3.11) and the initial condition (3.12). For the free evolution equation (3.10), note first that it is straightforward to check Here e i ∈ Z N has 1 in the i-th entry and 0 for the others. By multi-linearity of .
Here in the last equality above we used the fact that F i, j (w; x, y, t) only depends on the (N − i + 1)-th entry of x. Now (3.10) follows from linearity of the integral. Next we check the boundary conditions (3.11).
Then the boundary conditions (3.11) can be expressed as for k = 1. We prove (3.23) first. Note that for 2 ≤ k ≤ N and x k−1 = x k + 1 we have Hence by multi-linearity we have The determinant is 0 since there are two proportional rows. Now (3.23) follows from linearity of the integral.
The proof of (3.24) is similar. The only thing changes is that when k = 1, we have M i, j ( x, y, t)] = 0 since row 1 and N are proportional. Note that in the last equality of (3.25) we used the fact that w ∈ S z .
Finally we check the initial condition (3.12). We need to show . (3.26) Thanks to the cyclic-shift invariance of both sides of (3.26) (see (iv) of Proposition 3.4) we can assume without loss of generality that Then the two sides of (3.26) remain the same and where I 1 (i, j), I 2 (i, j), I 3 (i, j) are the integrals over the three contours, respectively. Here we recall that R and are chosen to be large (small) enough so that S z is contained in the region { < |w + 1| < R}\{|1 + pw| ≤ }. For I 3 (i, j) we note that For the other parts we write and (3.28) Depending on the analytical properties of the integrands in I 1 and I 2 we split into two cases: Case 1 x 1 < y 1 . First note the basic fact that for any x, y ∈ X (L) N we have Here in the last equality of (3.30) we take the outer integral contour to be |z| = r and let r → 0. 0,−1 is any simple closed contour with 0 and −1 inside and −1/ p outside. (3.31) Again by (3.29) we have We claim that the first inequality in (3.32) is strict and hence − This is due to our original assumption that But since it is independent of large enough R , we have I 1 (i, j) = 0 for all R large enough. Now a similar argument as in (3.30) with instead be a large circle |z| = r with r → ∞ implies: . (3.33) In conclusion we have reduced checking (3.26) to checking the following: for any x, y ∈ X But this is precisely equation (3.30) of [9] which appears in checking that the determinantal formula for the transition probability of discrete parallel TASEP on Z satisfies the initial condition. We will not repeat the proof here but point out that due to the assumption

A Toeplitz-like determinant formula
In this section we derive a formula for the finite-time multi-point joint distributions of discrete time parallel periodic TASEP under arbitrary initial condition. The proof basically follows the strategy of [4] by performing a multiple sum of transition probabilities over suitable particle configurations. The main technical part is a Cauchy-type identity for summation of left and right eigenfunctions (see Proposition 5.4) which generalizes Proposition 3.4 of [4] and some new difficulties appear.
where the contours for the integrals are nested circles 0 < |z m | < · · · < |z 1 |. Here Here k 0 = t 0 = a 0 := 0 and we suppress the dependence on a i , k i and t i 's in C (L) ( z) and D (L) y ( z).
Proof We start with the case m = 1 for a warm-up. First by Cauchy-Binet formula we can rewrite the transition probability (3.1) as Now to get the one-point distribution P N with x k ≥ a of the transition probability and interchange the order of integration and summation: (4.8) Now by Corollary 5.3 in Sect. 5 we have x∈X (L) Inserting (4.9) back to (4.8) and use Cauchy-Binet formula backwards we conclude that . (4.10) Here in order to interchange the order of summation and integration as in (4.8) we need to make sure that the summation over x ∈ X (L) N ∩ {x k ≥ a} converges absolutely. This is guaranteed if we assume |w + 1| > 1 for all w ∈ S z (see the discussion in Proposition 5.2). By choosing the contour to be a circle with large radius R we can make sure that |w + 1| > 1 for all w ∈ S z since for |z| = R and all w satisfying Finally a similar argument as in (ii) of Proposition 3.4 shows that the right-hand side of (4.10) does not depend on the choice of , so we can deform to be any simple closed contour containing 0. Now assume m ≥ 2. Then where t 0 := 0. As in the m = 1 case we rewrite the transition probability using Cauchy-Binet formula and interchange the order of summation and integration (will be justified later) so that (4.14) Now similar as in the m = 1 case we evaluate the sums appearing in (4.13) and (4.14) using Corollaries 5.3 and 5.6 in Sect. 5 and apply Cauchy-Binet formula backwards to conclude that (4.15) for C (L) ( z) and D (L) y ( z) defined in (4.2) and (4.3). Similar as the discussion before, in order to interchange summation and integration we need the absolute convergence of all the infinite sums which is guaranteed if we assume (4.16) By the same reasoning as in m = 1 case this can be achieved assuming the integral contours for z i 's are large nested contours |z | = r with r − r +1 also large enough for all 1 ≤ ≤ m. Finally we can deform the integral contours in (4.1) into arbitrary nested contours with 0 inside, not necessarily with large radius. This is due to the analyticity of C (L) ( z) and D (L) y ( z) in z for any z i 's nonzero and distinct, which can be shown in a similar way as (ii) of Proposition 3.4.

Toeplitz-like determinant to Fredholm determinant
The finite time formula obtained in Theorem 4.1 contains a factor D (L) y ( z) inside the integrals which is a Toeplitz-like determinant and is hard to take large-time limits. In this section we rewrite the formula based on a remarkable finite determinantal identity obtained in [5]. It is an identity between a Toeplitz-like determinant with symbol supported on a finite set and a Fredholm determinant with kernel acting on the 2 space supported on the same finite set.

Proof of Theorem 2.2
The proof is almost verbatim to the proof of Theorem 3.1 in [5] so we omit most of the details. We apply Proposition 4.1 of [5] with the functions p i (w), q i (w) and h i (w) given by for 1 ≤ i ≤ N and

Summation identities of eigenfunctions
In this section we state and prove the summation identities used in computing the multi-point joint distribution in Sect. 4. First we recall the left and right eigenfunctions defined in (4.5) and (4.6). They are certain (anti)-symmetric functions appearing naturally in the transition probability formula (3.1). Our summation identities should be viewed as Littlewood or Cauchy type identities for these symmetric functions, but over configuration spaces x N +L > x 1 > · · · > x N (partitions wrapped on a cylinder).

Definition 5.1 (Left and right eigenfunctions)
and p ∈ C, we define the functions x ( w) and r x ( w) for w ∈ C N as follows: The rest of this section is organized as follows: we list all the summation identities we need in Sect. 5.1 and postpone the proofs to the next few subsections. The proof of the summation identity over a single eigenfunction (Proposition 5.2) consists of straightforward manipulations of determinants and will be discussed in Sect. 5.2. The summation identity over left and right eigenfunctions (Proposition 5.4) is more involved and the proof will be splitted into three steps discussed in Sects. 5.3, 5.4, 5.5 and 5.6. Throughout the proof several rank-one perturbation formulas for Cauchy determinants are frequently used so we collect all these elementary formulas in a separate section for convenience, see Sect. 5.7 for details. Finally in Sect. 5.8 we discuss the proof of Corollaries 5.3 and 5.6 which are simple consequences of the previous propositions and the periodic nature of the identities.

Summation identities over eigenfunctions
We start with a summation identity for r x ( w): Proposition 5.2 (Summation over a single eigenfunction) Let z ∈ C be nonzero. Let r x ( w) be as in (4.6) where w = (w 1 , . . . , w N ) ∈ (S z ) N satisfies N j=1 |w j + 1| > 1. Then 3) The following corollary which allows slightly more general constraints on the summation is a simple consequence of Proposition 5.2 and the periodic nature of X for all 0 ≤ k ≤ N − 1 and a ∈ Z.
The above summation identities over a single eigenfunction are sufficient for computing the one-point distribution. To get the multi-point distribution we also need the following Cauchy-type summation identities over products of left and right eigenfunctions.

Remark 5.5
It is interesting to note that the right-hand side of (5.5) does not depend on p explicitly (of course the w i ' should satisfy certain algebraic equations which depend on p). Taking p → 0, Proposition 5.4 degenerates to Proposition 3.4 of [4], which can be understood as the Cauchy identity (a periodic version) for the Grothendieck polynomial (and its dual), and can be derived from the finite-sum Cauchy identity for Grothendieck polynomials obtained in Theorem 5.3 of [17]. For 0 < p < 1, to the best of our knowledge the corresponding Cauchy-type identity (5.5) has not been discussed in the existing literature, at least for the periodic case. The key point here is instead of summing over all configuration x with 0 ≤ x N < x N −1 < · · · < x 1 as in the usual Cauchy identity, we are only summing over those configurations satisfying the extra constraint x 1 < x N + L. For general spectral parameters w and w this sum only gives a deformed or generalized Cauchy determinant. It further reduces to a genuine Cauchy determinant when we impose the conditions as in Proposition 5.4 that the spectral parameters satisfy suitable Bethe equations.
for all 0 ≤ k ≤ N − 1 and a ∈ Z.

Proof of Proposition 5.2
First we write − (a, . . . , a). The summation over a converges absolutely provided that N j=1 |1 + w j | > 1. Thus we have reduced the computation to a summation with constraint x N = 0 which is a finite sum so there is no convergence issue. We will see the telescoping nature of the summation over a.
The main difficulty for computing the sum comes from the N i=1 (1 − p1 x i−1 −x i =1 ) factor in r x ( w) which gives different weights to different terms in the sum. To handle this we first split the sum according to the position of the first empty site to the right of x N in the configuration x. More precisely we set We will first compute the summations over X (L) N ,k for each k and then sum over k. The results for each of these steps are summarized in the following two lemmas, whose proofs will be at the end of this subsection.

Lemma 5.7 We have
Now by Lemma 5.8 and (5.7) we see that This completes the proof of Proposition 5.2.

Proof of Lemma 5.7 Note that configurations
We sum over configurations in X (L) N ,k in the following order: Note that Adding the (k + 3)-th column to the (k + 2)-th column we get Now we repeat this procedure and perform the sum over x N −k−2 , . . . , x 2 to get Finally note that L−1 where we used the fact that w i ∈ S z . Thus we conclude that

Proof of Lemma 5.8
For the purpose of summing over k we further rewrite R . Now we perform the sum over k in the order k = N − 1, N − 2, . . . , 0. Note that for each 0 ≤ k ≤ N − 1,R ( w) except for j = k + 1. Hence by multi-linearity of the determinants we have Here to simplify the expression we have multiplied the (N −1)-th column by − p(1− p) and added it to the sum of the N -th column ofR (N −1) ( w) andR (N −2) ( w), using the simple identity that Repeating this procedure we get Multiplying the j-th column of T (0) ( w) by − p and adding to the ( j + 1)-th column for j = N − 1, . . . , 2, we get This completes the proof of Lemma 5.8.

Proof of Proposition 5.4: Strategy
The proof of Proposition 5.4 is rather lengthy so we divide it into three steps and discuss them one by one in the next few subsections. The proof mainly follows the proof strategy of Proposition 3.4 of [4] but there are several new technicalities. The non-uniform term leads to extra difficulty for the sum. In Step 1 we overcome this by introducing a different way (and slightly more convenient way in our opinion) of expressing the sum in (5.5) comparing to the proof in [4], see Lemma 5.9 for details. In Step 2 we establish a key summation identity (see Lemma 5.12) which generalizes Lemma 5.4 in [4] while the computation is more delicate. Finally in Step 3 we combine the formula obtained in Step 1 and the summation identity obtained in Step 2 to conclude the final result.

Proof of Proposition 5.4: Step 1
Similar as in the proof of Proposition 5.2 we first write so that it suffices to compute the sum over x ∈ X (L) N ∩ {x N = 0}, which is a finite sum so there is no convergence issue. The summation over a converges absolutely by our assumption that N j=1 |w j + 1| > N j=1 |w j + 1|. Expanding the determinants in (5.9) Lemma 5.9 Let N ≥ 1 be an integer and w 1 , . . . , w N and w 1 , . . . , w N be distinct complex numbers. Then for any integer L > N we have (5.10) Here any empty product is set to be 1.

Proof
We use an induction on N . For N = 1 the identity is obvious. Assume now N ≥ 2 and the identity holds for all indices less than N . We split the sum into two sums depending on whether x 1 = L − 1 or not: For T 1 we first relabel the indices so that (x 0 , x 1 , . . . , x N −1 ) := (x 1 , x 2 , . . . , x N ) and L := L − 1. Then by induction hypothesis we have 1 ( w, w ) := 0 for convenience. For T 2 we calculate the sum directly in the order L − 1 > x 1 > · · · > x N = 0 using Lemma 5.10 below which gives: for all 1 ≤ k ≤ N . Using the simple identity which follows from Lemma 5.11 below by taking z j = w j +1 w j +1 and properly shifting the indices.
Proof This lemma is a slightly modified version of Lemma 5.3 of [4]. The proof is an elementary computation of geometric sums and almost identical to the proof in [4] so we omit it.

Lemma 5.11
Let n ≥ 1 be an integer and z 1 , . . . , z n be complex numbers such that k j= z j = 1 for all 1 ≤ ≤ k ≤ n. Then Proof Fixing an integer M > n, consider the following sum: We calculate the sum in two different ways: from right to left or from left to right. Namely we set Where the coefficients C r → k (z 1 , . . . , z n )'s and C →r k (z 1 , . . . , z n )'s are some very explicit functions in z 1 , . . . , z n independent of M which are analytic for all z i 's satisfying the assumption that k j= z j = 1 for all 1 ≤ k ≤ ≤ n. In particular it is straightforward to check (see Lemma 5.10 for example) that LHS of (5.18) = C r → n (z 1 , . . . , z n ), RHS of (5.18) = C →r n (z 1 , . . . , z n ).
We claim that C r → k (z 1 , . . . , z n ) = C →r k (z 1 , . . . , z n ) for all 0 ≤ k ≤ n. This in particular implies (5.18). Due to analyticity it suffices to check this for the z j 's satisfying |z j | < 1 for all 1 ≤ j ≤ n. In this case by letting M → ∞ in the equality Repeat this procedure we see C r → k = C →r k for all 0 ≤ k ≤ n.

Proof of Proposition 5.4: Step 2
In this section we simplify the sum (5.9). We rewrite the sum further by first choosing two index sets J and J with |J | = |J | = k and then expressing the sum in terms of summation over index sets J , J : (5.21) By Lemma 5.12 below (which is of interest on its own) we have To simplify H 2 (J c , (J ) c ) we use the assumption that w ∈ (S z ) N and w ∈ (S z ) N . Namely for all 1 ≤ i, i ≤ N we have (5.24) Now apply Lemma 5.12 again with the role of w and w exchanged we have Here the extra factor in (5.25) comes from the fact that in (5.24) the product starts from i = 1 instead of i = 2 and the exponent is i instead of i − 1 comparing to (5.26). This completes Step 2 of the proof.

Lemma 5.12
Let n ∈ N. Given any complex numbers w i and w i , i = 1, . . . , n such that w i = w i , for all 1 ≤ i, i ≤ n. Then for any 0 < p < 1 we have = LHS of (5.26), (5.27) assuming all the infinite geometric series converge absolutely.

Proof of Lemma 5.12
The proof is based on induction on n. The main tools are several rank-one perturbation formulas for the Cauchy determinants. Since we will use them several times we will collect them in a separate section, see Sect. 5.7. For n = 1, (5.26) is trivial. Let n ≥ 2 and assume (5.26) is true for all indices ≤ n − 1. Given σ, σ ∈ S n we first choose two indices = σ (1) and = σ (1) and shift the restriction of σ and σ on {2, . . . , n} by 1 but still denote them by σ and σ . Then LHS of (5.26) = n , =1 By the induction hypothesis the term inside the big bracket above equals ⎡ Hence LHS of (5.26) (5.28) Here we set A(z) := n i=1 (z − w i ) and B(z) := n i=1 (z − w i ) and D 1 , D 2 represent the first and second sum over , respectively. Now by (5.44) and Lemma 5.17 we have , (5.29) and (5.30) Inserting (5.29) and (5.30) into (5.28), after necessary cancellation we obtain This completes the proof of Lemma 5.12.
Finally inserting (5.22) and (5.25) into (5.19) and apply Lemma 5.14 below we obtain Note that in the first equality of (5.32) we have added an extra term corresponding to |J | = |J | = 0 comparing to (5.19) which is harmless since the summand is 0 in this case.

Proof of Proposition 5.4: Step 3
In this section we further simplify equation (5.32) to conclude the proof of Proposition 5.4. We compute det[Ĉ(i, i )] first. For notational convenience we set (z /z) L := μ. Note thatĈ Hence by Lemma 5.16 and Lemma 5.17 we have whereĈ ,k is the matrix obtained by removing row and column k fromĈ. Similarly (5.37) Now sinceĈ ,k has the same entries asĈ only omitting row and column k, by (5.35) we have Where C ,k is obtained from removing row and column k from the N × N Cauchy . .
Inserting (5.41) into (5.39) and (5.40) and combine with (5.32), after some tedious simplification we conclude that This completes the proof of Proposition 5.4.

Perturbation formulas for Cauchy determinants
In this section we collect all the elementary linear algebra facts needed in the proof of Proposition 5.4. Some of them have already been discussed in [4]. First we state a general linear algebra lemma on rank-one perturbations: i, j=1 be an n × n matrix. Then for any function f , g : C → C and complex numbers x 1 , . . . , x n , y 1 , . . . , y n we have where D ,k is obtained by removing row and column k from D.
Proof For D invertible, by the rank-one property and the Cramer's rule we have For general matrix D we pick { k } ∞ k=1 such that k → 0 as k → ∞ and D + k I n are invertible for all k . Now apply the above argument for D + k I n and let k → ∞.
Next we specialize to the case of C being a Cauchy matrix when the minors can be explicitly calculated: 1104 Y. Liao 1 x i −y j for distinct complex numbers x 1 , . . . , x n and y 1 , . . . , y n . Then we further have

monic polynomials with roots at x i 's and y i 's.
Proof For Cauchy matrix C we have .
Note that C ,k is also a Cauchy matrix so we have Hence . For special choices of f and g, (5.44) can be further simplified using the residue theorem. We list here all the situations encountered in the proof of Proposition 5.4.

Lemma 5.17
Given distinct complex numbers x 1 , . . . , x n and y 1 , . . . , y n . Let C be the Cauchy matrix with (i, j)-th entry 1 (5.46) (2) For f (x) = x 1+ px and g(y) = 1+ py y(1+y) we have n Proof We will only prove part (1) since the arguments for the other parts are similar.
we consider the following double contour integral: Where R > r are both large enough so that all the possible poles of the integrand are inside the integral contours. Now since for fixed r the integrand is of order O(R −2 ), the double integral goes to 0 as R → ∞. Thus for all R large enough the double contour integral equals 0. On the other hand by the residue theorem we have .
Where in the first equality the first contour integral is O(R −2 ) hence 0 for R large enough. The single sum over 1 ≤ k ≤ n can be obtained as the residue terms for the following single contour integral: Here r is large enough so that |ξ | = r contains all possible poles of the integrand inside. On the other hand by considering the residue at ∞ we have |ξ |=r dξ 2πi Similar residue analysis on the contour integral |ξ |=r Inserting (5.50) and (5.51) into (5.49) we obtain This completes the proof of part (1).

Proof of Corollary 5.3 and 5.6
Finally we discuss how Corollaries 5.
We move the first k columns of the matrix to the end. The resulting determinant equals (−1) k(N −1) times the determinant of the matrix whose (i, j)-th entry has the form After factoring out the common factor (1 + w i ) k from each row and z L from the last k columns we conclude that (5.54) Corollary 5.6 follows from Proposition 5.4 in a similar way.

Large-time asymptotics under relaxation time scale
In this section we discuss the large time limit of the multi-point distribution of dpTASEP(L, N , y) under the relaxation time scale t = O(L 3/2 ). In Theorem 2.4 we state the limit theorem for general initial condition satisfying certain assumptions.
Below are the precise assumptions on the initial conditions we need:

Assumptions on the initial condition
We now state the assumptions on the sequence of the initial conditions y(L) under which we prove the limit theorem. The conditions are in terms of the global energy function and the characteristic function defined in Definition 2.12.
Recall that the finite-time formula (2.4) involves a m-fold contour integral with respect to z = (z 1 , . . . , z m ) for 0 < |z m | < · · · < |z 1 | < r c where the critical radius r c is defined in (2.16). It turns out that we need to rescale the parameters z i to be close to the critical value r c in order to make the large time limits converge. Notation 6.1 For given complex parameter z with 0 < |z| < r c , we introduce the rescaled parameter z defined by The constraint 0 < |z| < r c then translates to 0 < |z| < 1. We introduce a similar rescaling for the parameters z = (z 1 , . . . , z m ). Then for 0 < |z m | < · · · < |z 1 | < r c , the rescaled parameters satisfy 0 < |z m | < · · · < |z 1 | < 1. Throughout the rest of the paper we will always use z i to represent the unscaled parameters and z i to represent the rescaled ones satisfying (6.1).

Assumption 6.2
We assume that the sequence of the initial profiles y = y(L) satisfies the following three conditions as L → ∞.

Step and flat initial conditions
It turns out that Assumption 6.2 is not easy to check in general. Nevertheless we are able to verify them for at least the classical step and flat initial conditions. The following proposition combined with Theorem 2.4 gives the corresponding limit theorems for dpTASEP starting with step and flat initial conditions. for 0 < |z| < 1, where B(z) is defined in (6.14), h(ζ, z) is defined in (6.15) and (6.16).
The step case is trivial since by the discussion in Remark 2.13 we have E step (z) = χ step (v, u; z) ≡ 1. The calculation for the flat case is a bit more involved and we postpone the proof to Sect. 7.5.

Formula for the limiting distribution
The following formula for the relaxation-time limiting distribution was first obtained in [4] for the step initial condition and [5] for more general initial conditions. The formula involves C and where L z and R z are the sets defined in Definition 6.4. We express the limiting distribution function F per ic in terms of the above terms. Definition 6.4 Given 0 < |z| < 1, we define the discrete sets S z := L z ∪ R z where dz m 2π iz m · · · dz 1 2π iz 1 , (6.8) where z = (z 1 , . . . , z m ) and the contours are nested circles satisfying 0 < |z m | < · · · < |z 1 | < 1 and also, r 1 < |z 1 | < r 2 with r 1 , r 2 being the constants in Assumption 6.2 (B). The first function in the integrand is given by The second function is given by The function C per step ( z) and kernels K per step,1 and K per step,2 are first obtained in [4] and the definitions will be recalled in the next two sections for completeness.

The factor C per step ( z)
Let Li s (z) be the polylogarithm function which is defined by dt, (6.11) for |z| < 1 and s ∈ C. Set Let log z denote the principal branch of the logarithm function with branch cut R ≤0 . Set √ kk (6.13) for 0 < |z|, |z | < 1 where the integral contours are the vertical lines Re(ξ ) = a and Re(η) = b with constants a and b satisfying − √ − log |z| < a < 0 < b < − log |z |.
Note that C per step ( z), and hence C per ic ( z), depend on x i and τ i , but not the spatial parameters γ i . For each i, define (6.17) where we set τ 0 = γ 0 = x 0 = 0. We also define where we set z 0 = z m+1 = 0.
Definition 6.7 Let S 1 and S 2 be the discrete sets defined in (6.5) and (6.6). Let K per step,1 : 2 (S 2 ) → 2 (S 1 ) and K per step,2 : 2 (S 1 ) → 2 (S 2 ) denote the operators defined by their kernels Here L z i and R z i are again defined in Definition 6.4.

Proof of Theorem 2.4
In this section we discuss the proof of Theorem 2.4. The ideas are similar to the one in [4,5] so we omit some technical details. Clearly the theorem follows immediately from the following two lemmas, dealing with the asymptotics of C The rest of the section is devoted to proving Lemmas 7.1 and 7.2. We start with a discussion on the asymptotic behaviors of the roots of the Bethe polynomial q z (w) = w N (1 + w) L−N − z L (1 + pw) N under the critical re-scaling in Sect. 7.1. Then in Sect. 7.2 we list a few lemmas discussing the asymptotics of several products involving these roots under the critical re-scaling. With these preparations we prove Lemmas 7.1 and 7.2 in Sects. 7.3 and 7.4, respectively. Finally in Sect. 7.5 we verify the Assumption 6.2 for the classical step and flat initial conditions.

Asympotics of the Bethe roots
We assume that the particle density ρ := N /L stays within a compact subset of (0, 1) for all L. From the discussion in Sect. 2.3 we know the level set {w ∈ C : |w N (1 + w) L−N | = |z L (1 + pw) N |} consists of two disjoint closed contour for |z| < r c so we can define: Definition 7.3 Given |z| < r c , we define two closed contours L and R by where D(a, r ) is a disc centered at a with radius r and c 0 is defined in (7.5). Then for the re-scaled parameter Similar results hold if we replace L z and L z by R z and R z .
Proof This lemma is a minor generalization of Lemma 8.1 of [3] by allowing one extra parameter p in the Bethe equation. The proof is almost identical to the one in [3] so we omit the details.

Asymptotics of various products over Bethe roots
In this section we collect all the results involving limits of products of the Bethe roots appearing in the finite-time formula that are independent of the parameters y, a i , t i and k i . The starting point is the following simple integral formula for the sums of functions evaluated at the left or right Bethe roots.

Similarly if φ(w) be a function analytic in the interior and a neighborhood of L , then
J (w) . (7.10) Here L and R are simple closed contours lie in the half-plane {w : Re(w) < w c } (respectively {w : Re(w) > w c }) with L (respectively R ) inside. Taking the function φ to be the constant function 1 implies in particular that L dw 2πi Hence by the residue theorem we know The proof for (7.10) is similar.
As taking logarithm transforms products into sums, the following lemma is a direct consequence of Lemma 7.5 and the method of steepest descent.

Lemma 7.6
Given |z| < r c and z = (−1) N z L r −L c . Suppose ρ = N /L stays in a compact subset of (0, 1), then for every 0 < < 1/2 the following holds for all large enough L.
There is a constant C > 0 such that for every w satisfying |w − w c | ≥ L −1/2 , , which comes from the larger root of the quadratic equation p(L − N )w 2 + Lw + N = 0, as opposed to w = −ρ for the p = 0 degeneration discussed in [3]. A standard steepest descent analysis using integral representations obtained in Lemma 7.5 with critical point w c yields all the estimates. We omit the details.

Asymptotics of C (L) y ( z)
Now we are ready to prove Lemma 7.1. Recall that (see Definition 2.14), C (L) as L → ∞ due to Assumption 6.2. Hence it suffices to prove where C per step ( z) is defined in Definition 6.6 and C (L) step ( z) is defined as follows (see also Definition 2.14): where a, k ∈ Z and t ∈ N are given parameters. Then for z L = (−1) N r L c z with 0 < |z| < 1 and the parameters satisfying (7.19) we have for L large enough and fixed 0 < < 1/2 where A 1 (z) and A 2 (z) are scaled polylogarithm functions as in (6.12). The constants c i are the same as in (2.8) and the re-scaled parameters are chosen such that τ > 0, γ ∈ [0, 1] and x ∈ R.
. Apply Lemma 7.5 to the two sums over left and right Bethe roots and deform both contours L and R to the vertical line with real part where G(w) = (−k) log(−w) + (a + k) log(w + 1) + (−t + k) log(1 + pw). Note that in (7.21) we have used the fact discussed in Lemma 7.5 that L dw 2π i Splitting the integral representation for log E(z) into two parts with |ζ | ≤ L /4 and |ζ | > L /4 and using the estimates for z L (1+ pw) N q z (w) and 1 J (w) obtained in Lemma 7.6 we see for some constants c, α > 0. Now (7.20) follows from integral representations of polylogarithm (6.11).

Asymptotics of D (L) y ( z)
Next we discuss the asymptotics of the Fredholm determinant part D y (z). Note first that by a standard series expansion of Fredholm determinants we have y, n ( z), (7.24) where n = (n 1 , . . . , n m ) and if k = n 1 + · · · + n −1 + k for integer k ≤ n with odd , v ( ) k if k = n 1 + · · · + n −1 + k for integer k ≤ n with even , (7.26) and if k = n 1 + · · · + n −1 + k for integer k ≤ n with odd , u ( ) k if k = n 1 + · · · + n −1 + k for integer k ≤ n with even . (7.27) A similar series expansion holds for the limiting Fredholm determinant D per ic ( z) with L z and R z replaced by the limiting roots L z and R z and the kernels replaced by the limiting kernels K It is clear that Lemma 7.2 follows immediately from Lemma 7.8 by dominated convergence theorem. To prove Lemma 7.8 we will prove the convergence of the kernels after proper conjugation for points inside the critical region as well as exponential decay estimates for the kernel at points outside the critical region. The conjugation is as follows: we replace K We define the square root to be √ w = r 1/2 e iθ/2 for w = re iθ with −π < θ ≤ π . Note that the product of determinants will always be continuous even though the square root function is not since every (f i ) 1/2 is multiplied twice. We change the limiting kernels in a similar way by replacing f i or f j in the kernels with √ f i f j and denote them as K per 1 andK per ic . Then we have the following asymptotics for the conjugated kernels, which easily implies Lemma 7.8. Lemma 7.9 Fix 0 < < 1/(1 + 2m). Let Proof The proof is similar to the proof of Lemma 10.2 and Lemma 10.3 in [5]. The key observation is the existence of a d − 1 to 1+ pv . Using this relation we can express the global energy and characteristic functions in terms of products over the Bethe roots. We omit the details. Combining Lemma 7.11 with the asymptotics obtained in Lemma 7.6 we can now prove Theorem 6.3.
Proof of Theorem 6.3 For fixed 0 < < 1/2 by Lemma 7.6 (viii) we have Similarly by Lemma 7.6 (v) we have On the other hand we have p(d − 1)v 2 + dv Hence by Lemma 7.6 (iv) and (v) we have where we used the fact that w c · w * c = 1 p(d −1) and also e h(0,z) = (1 − z) 1/2 which follows from (6.15). Combining (7.37), (7.38) and (7.39) we conclude that The argument for the characteristic function part is quite similar. To verify part (B) of Assumption 6.2 we note that given 0 < < 1/8, for u ∈ L ( ) This then implies that Finally by Lemma 7.6 (iii) we have Finally part (C) in Assumption 6.2 is clearly true since by Lemma 7.6 every factor in for . Thus |χ flat (v, u; z)| ≤ C · L and part (C) of Assumption 6.2 is satisfied.

Multi-time distribution of discrete time parallel TASEP on Z
So far we have been focusing on the so-called relaxation time scale when the period L is comparable to t 2/3 and the height fluctuation is critically affected by the boundary conditions. In this section we discuss another important regime, the sub-relaxation time scale L t 2/3 when the particles are not affected by the boundary and are expected to behave the same as when the underlying space has infinite-volume.

Discrete time parallel TASEP with large period
Instead of taking the re-scaled time parameter τ → 0 in the relaxation-time limiting distribution formula, we will start with the pre-limit joint distribution formula (2.4). The following elementary proposition confirms the intuition that for fixed spatial and time parameters a i , k i and t i , when L is large enough, the joint distribution P y is the probability with respect to discrete time parallel TASEP on a periodic domain with period L and P (∞) y is the probability with respect to an independent discrete time parallel TASEP on the infinite lattice Z, both starting with the initial condition satisfying x i (0) = y i for 1 ≤ i ≤ N .
Proof The proof is almost identical to the similar statement for continuous time TASEP as in Lemma 8.1 of [4] where the particle ordering is reversed. We omit the details.
In principle for every L satisfying (8.1), the left-hand side of (8.2) which has an expression as contour integrals of Fredholm determinant given by (2.4) gives a formula for the multi-point joint distribution of discrete time parallel TASEP on Z for fixed parameters a i , k i and t i . The main obstacle here is the dependence on the extra parameter L for the left-hand side of (8.2). The key idea to get rid of the dependence on L was already illustrated in Proposition 3.4 (iii) and can be summarized as follows: (1) By rewriting the sum over the Bethe roots as contour integrals using the residue theorem we get an analytic continuation of the integrand (for the outer integral with respect to z) to {|z| > 0}.
(2) Under proper assumption on the parameter L (in fact exactly (8.1)), one can further show the integrand is analytic at z = 0, thus the outer integral with respect to z equals the evaluation of the integrand at z = 0, which turns out to be independent of L.
Now as in the above discussion, we would like to send z 1 , . . . , z m all to 0 in (2.4) while still keep the nested relation 0 < |z m | < · · · < |z 1 | < r c . For this we change the integral variables as follows: for given L satisfying (8.1), we set where z 0 := 1 and 0 < |z m | < · · · < |z 1 | < r c . Then the corresponding relations between θ i 's are precisely 0 < |θ 0 | < r L c := r c and 0 < |θ i | < 1 for 1 ≤ i ≤ m − 1. After this change of variable, the finite-time multi-point joint distribution formula (2.4) can be rewritten as where the integral contours are circle centered at origin with radius smaller than r c for θ 0 and circles centered at origin with radii less than 1 for θ i , 1 ≤ i ≤ m − 1. Here the integrand is defined such that

Proof of Theorem 2.5
As discussed above, Theorem 2.5 follows immediately once we establish the analyticity ofĈ y (θ 0 , . . . , θ m−1 ) at θ 0 = 0 so that we can get rid of the integral with respect to θ 0 whose value equals the integrand evaluated at θ 0 = 0. The precise statement needed are summarized in the following two lemmas which will be proved in the next two subsections. Finally using (8.2) we complete the proof of Theorem 2.5.

Proof of Lemma 8.2
By the definition of the function C  . . . , z m )A 2 (z 1 , . . . , z m )A 3 (z 1 , . . . , z m ). (8.8) where E (z) = u∈L z (−u) −k v∈R z (v + 1) −a −k ( pv + 1) t −k and the functions A i represents the products inside each brackets. The parameters θ i and z i are related by the equations z L i = i−1 j=0 θ j for 1 ≤ i ≤ m. Now if we send θ 0 → 0, then all the z i 's will converge to 0. As a result, all the roots in L z i converge to −1 while all the roots in R z i converge to 0, for all 1 ≤ i ≤ m. It is then easy to see that for j = 1, 2, 3

Proof of Lemma 8.3
This Lemma is an analogue of Lemma 5.2 of [15] and the proof is similar so we omit most of the details. What we want to prove is an identity between two Fredholm determinants, with the first one acting on discrete sets consisting of roots of certain polynomials and the other one acting on continuous contours. By looking at the standard series expansions of both Fredholm determinants, it suffices to show that every term in the two series expansions matches. These type of identities between multiple sums over roots of certain algebraic equations and multiple contour integrals can be understood as a highly nontrivial consequence of the residue theorem, which is stated and proved with enough generality for our purpose in Proposition 4.4 of [15]. We take q(w) = w N (w + 1) L−N (1 + pw) −N in the proposition and follow the same argument as in Lemma 5.2 of [15].

KPZ scaling limit
In this final section we briefly discuss the large-time asymptotics for the multi-time distribution of discrete time parallel TASEP on Z under the 1 : 2 : 3 KPZ scaling. For completeness we first recall the definition of the limiting distribution function F kpz step which essentially agrees with the function F step appearing in Theorem 2.20 of [15] with a slight different choices of the parameters. The structure of F kpz step is similar to the finite-time distribution as contour integrals of a Fredholm determinant D kpz step (θ 1 , . . . , θ m−1 ).

Spaces of the operators
Given integer m ≥ 1. We first define the contour in the complex plane Now we deform the contours j to be sufficiently close to the critical points w c = − 1 1+ √ q such that locally they look like the limiting contours j . More concretely we deform + j,R such that (2) There exists constant C > 0 independent of T and n such that for all n ∈ N dw 1 2π i · · · dw n 2π i ≤ C n .
Acknowledgements I would like to thank Jinho Baik for suggesting this problem and for useful discussions. I would also like to thank Zhipeng Liu and Mustazee Rahman for helpful discussions. I am grateful to the two anonymous referees for catching up several typos and to their valuable suggestions which help improve the quality of the paper. This work is supported in part through Jinho Baik's NSF grants DMS-1664531.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.