PushTASEP in inhomogeneous space

We consider the PushTASEP (pushing totally asymmetric simple exclusion process, also sometimes called long-range TASEP) with the step initial configuration evolving in an inhomogeneous space. That is, the rate of each particle's jump depends on the location of this particle. We match the distribution of the height function of this PushTASEP with Schur processes. Using this matching and determinantal structure of Schur processes, we obtain limit shape and fluctuation results which are typical for stochastic particle systems in the Kardar-Parisi-Zhang universality class. PushTASEP is a close relative of the usual TASEP. In inhomogeneous space the former is integrable, while the integrability of the latter is not known.


Introduction and main results
1.1. Overview. The totally asymmetric simple exclusion process (TASEP) was introduced about 50 years ago, independently in biology [MGP68], [MG69] and probability theory [Spi70]. The latter paper also introduced zero-range processes and long-range TASEP. The long-range TASEP (which we call PushTASEP following more recent works) is the focus of the present paper.
Since early works, TASEP and PushTASEP were often studied in parallel. Once a result (such as description of hydrodynamics and local equilibria [Lig05], limiting density [Ros81], or asymptotic fluctuations [Joh00]) for TASEP is established, it can often be generalized to PushTASEP using similar tools. See [DLSS91] for hydrodynamics and related results for the PushTASEP (viewed as a special case of the Tooms model), and, e.g., [DW08], [BF08] for fluctuation results. Borodin and Ferrari [BF14] introduced a two-dimensional stochastic particle system whose two different one-dimensional (marginally Markovian) projections of are TASEP and PushTASEP. This coupling works best for special examples of initial data, most notably, for step initial configurations. It is worth pointing out that most known asymptotic fluctuation results for TASEP, PushTASEP, and related systems (in the Kardar-Parisi-Zhang universality class, cf. [Cor12], [QS15]) require integrability, that is, the presence of some exact formulas in the pre-limit system.
Running either TASEP or PushTASEP in inhomogeneous space (such that the particles' jump rates depend on their locations) is a natural generalization. Hydrodynamic approach works well in macroscopically inhomogeneous systems, and allows to write down PDEs for limiting densities [Lan96], [Sep99], [RT08], [GKS10], [Cal15]. This leads to law of large numbers type results for the height function (in particular, of the inhomogeneous TASEP). However, when the disorder is microscopic (such as just one slow bond), this affects the local equilibria, and makes the analysis of both limit shape and asymptotic fluctuations of TASEP much harder [BSS14], [BSS17]. Overall, putting TASEP on inhomogeneous space breaks its integrability.
On the other hand, considering particle-dependent inhomogeneity (when the jump rate depends on the particle's number, but not its location) in TASEP preserves its integrability, and allows to extract the corresponding fluctuation results, cf. [Bai06], [BFS09], [Dui13].
The main goal of this paper is to show that, in contrast with TASEP, the PushTASEP in inhomogeneous space started from the step initial configuration retains the integrability for arbitrary inhomogeneities. Namely, we obtain a matching of the PushTASEP to a certain Schur process, which follows by taking a third marginally Markovian projection of the two-dimensional dynamics of Borodin-Ferrari [BF14], which was not observed previously. This coupling is also present in the Robinson-Schensted-Knuth insertion -a mechanism originally employed to obtain TASEP fluctuations in [Joh00]. The coupling of inhomogeneous PushTASEP to Schur processes, together with their determinantal structure [Oko01], [OR03] leads to exact formulas for the PushTASEP. We illustrate the integrability by obtaining limit shape and fluctuation results for PushTASEP with arbitrary macroscopic inhomogeneity.
Remark 1.1. Based on the tools employed in the present work, one can even say that our results could have been observed already in the mid-2000s. However, it is the much more recent development of stochastic vertex models, especially their couplings in [BBW18], [BM18], [BP17] to Hall-Littlewood processes, that prompted the present work (as a t = 0 degeneration of the Hall-Littlewood situation). Asymptotic behavior of the Hall-Littlewood deformation of the PushTASEP (in a homogeneous case) was studied in [Gho17].
Other examples of integrable stochastic particle systems in one-dimensional inhomogeneous space have been recently studied in [BP18], [KPS19]. These systems may be viewed as analogues of q-TASEP or TASEP, respectively in continuous space. (The q-TASEP is a certain integrable qdeformation of TASEP [BC14].) In those inhomogeneous systems, a certain choice of inhomogeneity leads to interesting phase transitions corresponding to formation of "traffic jams", when the density goes to infinity. In PushTASEP the density is bounded by one, and so we do not expect this type of phase transitions to appear. A two-dimensional stochastic particle system in inhomogeneous space unifying both the inhomogeneous PushTASEP considered in the present paper, and a process like the one in [KPS19], is studied in the forthcoming work [Ass19].
In the rest of the introduction we give the main definitions and formulate the results.
1.2. PushTASEP in inhomogeneous space. Fix a positive speed function ξ • = {ξ x } x∈Z ≥1 , uniformly bounded away from zero and infinity. By definition, the PushTASEP is a continuous time Markov process on particle configurations (1.1) on Z ≥1 (at most one particle per site is allowed). We consider only the step initial configuration x i (0) = i for all i ≥ 1, so at all times the particle configuration has the leftmost particle.
The system evolves as follows. At each site x ∈ Z ≥1 there is an independent exponential clock with rate ξ x (i.e., the mean waiting time till the clock rings is 1/ξ x ). When the clock at site x ∈ Z ≥1 rings and there is no particle at x, nothing happens. Otherwise, let some particle x i be at y. When the clock rings, x i jumps to the right by one. If the destination x + 1 is occupied by x i+1 , then x i+1 is pushed to the right by one, which may trigger subsequent instantaneous pushes. That is, if there is a packed cluster of particles to the right of x i , i.e., x i+m − m = . . . = x i+1 − 1 = x i and x i+m+1 − 1 > x i+m for some m ∈ {1, 2, . . .} ∪ {+∞}, then each of the particles x i+ , = 1, . . . , m, is instantaneously pushed to the right by one. The case m = +∞ corresponds to pushing the whole right-infinite densely packed cluster of particles to the right by one. Clearly, the evolution preserves the order (1.1) of particles. See Figure 1 for an illustration.
Thus described Markov process on particle configurations in Z ≥1 is well-defined. Indeed, consider its restriction to {1, . . . , N } ⊂ Z ≥1 for any N . This is a continuous time Markov process on a finite space in which at most one exponential clock can ring at a given time moment. For different N , such restrictions are compatible, and thus the process on configurations in Z ≥1 exists by the Kolmogorov extension theorem. Note however that the number of jumps in the process on Z ≥1 during each initial time interval [0, ε), ε > 0, is infinite. 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 rate= ξ 3 rate= ξ 12 Figure 1. Two possible transitions in the inhomogeneous PushTASEP. The first one corresponds to activating the particle at 3 which happens at rate ξ 3 . This particle then pushes two other particles (located at 4 and 5) to the right by one. The second transition corresponds to activating the particle at 12 at rate ξ 12 .
1.3. Determinantal structure. The height function of the PushTASEP is defined as (1. 2) The step initial condition corresponds to h(0, N ) = N for all N ≥ 0.
and the points (t i , N i ) are pairwise distinct.
Remark 1.3. Down-right paths are also called space-like (as opposed to time-like, when both t i and N i increase). These names come from a growth model reformulation, cf. [DLSS91], [Fer08].
Define a kernel depending on the speed function ξ • by The integration contours are positively oriented simple closed curves around 0, the w contour additionally encircles {ξ x } x∈Z ≥1 , and the contours satisfy |z| > |w| for t ≤ t and |z| < |w| for t > t . (Throughout the text 1 A stands for the indicator of an event A. By 1 without subscripts we will also mean the identity operator.) Fix a down-right path p = {(t i , N i )} r i=1 , and define the space For y ∈ X i set t(y) = t i , N (y) = N i .
Definition 1.4. We define a determinantal random point process 1 L p on X with the correlation kernel expressed through K (1.4). Namely, for any m ≥ 1 and any pairwise distinct y 1 , . . . y m ∈ X , let the corresponding correlation function of L p be given by The process L p exists because it corresponds to column lengths in a certain specific Schur process, see Section 2 for details.
On each X i = Z the random point configuration L p almost surely has a leftmost point, denote it by (t i , N i ). The joint distribution of these leftmost points is identified with the inhomogeneous PushTASEP. The following theorem is the main structural result of the present paper.
Theorem 1.5. Fix an arbitrary down-right path p = {(t i , N i )} r i=1 . The joint distribution of the PushTASEP height function along this down-right path is related to the determinantal process L p defined above as Corollary 1.6. For any t ≥ 0, N ≥ 1, and y ≥ 0 we have where the second equality is the series expansion of the Fredholm determinant given in the first equality. Similar Fredholm determinantal formulas are available for joint distributions of the PushTASEP height function along down-right paths.
Theorem 1.5 is a restatement of a known result on how Schur processes appear in stochastic interacting particle systems in (2 + 1) dimensions. Via Robinson-Schensted-Knuth (RSK) correspondences, such connections can be traced to [VK86], and they were heavily utilized in probabilistic context starting from [BDJ99], [Joh00], [PS02]. Markov dynamics on particle configurations coming from RSK correspondences were studied in [Bar01], [O'C03b], [O'C03a]. Another type of Markov dynamics whose fixed-time distributions are given by Schur processes was introduced in [BF14], and it, too, can be utilized to obtain Theorem 1.5. A self-contained exposition of the proof of this theorem following the latter approach is presented in Section 2.
Remark 1.7 (Connection to vertex models). Yet another alternative way of getting Theorem 1.5 is to view the PushTASEP as a degeneration of the stochastic six vertex model [GS92], [BCG16]. The latter was recently connected to Hall-Littlewood processes [Bor18], [BBW18], [BM18]. Setting the Hall-Littlewood parameter t to zero leads to a distributional mapping between our inhomogeneous PushTASEP and Schur processes.
1.4. Hydrodynamics. Let the space and time in the PushTASEP, as well as the speed function scale as follows: (1.7) where L is the large parameter going to infinity, and ξ(·) is a fixed positive limiting speed function bounded away from zero and infinity. Under (1.7), one expects that the height function h(t, N ) admits a limit shape (i.e., law of large numbers type) behavior of the form h(τ L, ηL ) L → h(τ, η), in probability as L → +∞. (1.8) Let us first write down a partial differential equation for the limiting density using hydrodynamic arguments as in [AK84], [Rez91], [Lan96], [GKS10]. Because of our scaling (1.7), locally around every scaled location ηL the behavior of the Push-TASEP (when we zoom at the lattice level) is homogeneous with constant speed ξ = ξ(η). Thus, locally around ηL the PushTASEP configuration should have a particle distribution on Z which is invariant under shifts of Z and is stationary under the speed ξ homogeneous PushTASEP dynamics on the whole line.
A classification of translation invariant stationary distributions for the PushTASEP is available [Gui97], [AG05]. Namely, ergodic (= extreme) such measures are precisely the Bernoulli product measures. For the Bernoulli product measure of density ρ ∈ [0, 1], the flux (= current) of particles in the PushTASEP (i.e., the expected number of particles crossing a given bond in unit time interval) is readily seen to be j(ρ) = ξρ 1 − ρ . Therefore, the partial differential equation for the limiting density has the form: The singularity at τ = 0 coming from the initial data corresponds to the fact that the PushTASEP makes an infinite number of jumps during every time interval [0, t], t > 0. That is, from t = 0 to t = τ L (for every τ > 0 in the regime L → +∞), the density of the particles drops below 1 everywhere.
Remark 1.8. One sees from, e.g., [BF14] that in the homogeneous case ξ(η) ≡ 1, a solution to (1.9) has the form (1.10) The condition η ≥ τ for nonzero density comes from the behavior of the leftmost particle in the PushTASEP which performs a simple random walk. Integrating this density in η gives the limiting Next, we present a solution to (1.9) for general ξ(·).
1.5. Limit shape. For any η > 0, set this is the rescaled time when the leftmost particle in the PushTASEP reaches ηL . Consider the following equation in z: (1.12) Lemma 1.9. For any η > 0 and τ ∈ (0, τ e (η)) equation (1.12) has a unique root z on the negative real line.
Proof of Lemma 1.9. This is evident due to the strict increasing of the right-hand side of (1.12) in z ∈ (−∞, 0), and the fact that at z = 0 this right-hand side is equal to τ e (η) given by (1.11).
where z(τ, η) comes from Lemma 1.9. We call h the limiting height function.
Remark 1.11. 1. Since the right-hand side of (1.11) depends on z and η in a continuous way when z ≤ 0, the function η → z(τ, η) is continuous for each fixed τ . Thus, the height function (1.13) is also continuous in η. (This continuity extends to the unique η e such that τ e (η e ) = τ because both cases in (1.13) give zero.)

Equivalently, the function
One can check that the limiting density corresponding to h(τ, η) (defined as ρ(τ, η) = ∂ ∂η h(τ, η) when this derivative exists) is expressed through z as (1.14) One can also verify that ρ(τ, η) formally satisfies the hydrodynamic equation (1.9). This is done in Appendix A. See Figures 2 and 3 for illustrations of limit shapes of the height function and the density. Figure 2. Limiting density and height function for three cases of piecewise linear ξ(·).
Theorem 1.12 (Limit shape). Fix arbitrary τ, η > 0. If the limiting speed function ξ(·) is piecewise continuously differentiable on [0, η], then in the regime (1.7) we have the convergence Here h is the random height function of our PushTASEP, and h is defined by (1.13).
Theorem 1.12 follows from the fluctuation result (Theorem 1.13) which is formulated next.
1.6. Asymptotic fluctuations. Using the notation of Section 1.5, define for 0 < τ < τ e (η): (1.15) Note that this quantity is also continuous in η similarly to Remark 1.11.1. The following is the main asymptotic fluctuation result of the present paper: Theorem 1.13. Fix arbitrary η > 0 and τ ∈ (0, τ e (η)). If the limiting speed function ξ(·) is piecewise continuously differentiable on [0, η], then in the regime (1.7) we have the convergence This result implies the law of large numbers (Theorem 1.12).
Using the determinantal structure of Theorem 1.5, it is possible to also obtain (under slightly more restrictive smoothness conditions on ξ(·)) multipoint asymptotic fluctuation results along space-like paths. These fluctuations are governed by the top line of the Airy 2 line ensemble. We will not focus on this result as it is a standard extension of Theorem 1.13 which readily follows from the determinantal structure. 1.7. Outline. In Section 2 we describe the connection between the inhomogeneous PushTASEP and Schur processes, and establish Theorem 1.5. In Section 3 we perform asymptotic analysis and establish Theorems 1.12 and 1.13 on the limit shape and asymptotic fluctuations. In Appendix A we check that the limiting density ρ defined in Section 1.5 formally satisfies the hydrodynamic equation coming from the inhomogeneous PushTASEP.
1.8. Acknowledgments. I am grateful to Konstantin Matveev for an insightful remark on connections with RSK column insertion, and to Alexei Borodin, Alexey Bufetov, Patrik Ferrari, Alisa Knizel, and Axel Saenz for helpful discussions. The work was partially supported by the NSF grant DMS-1664617.

Schur processes and inhomogeneous PushTASEP
Here we present a self-contained proof of Theorem 1.5 which follows from results on Schur processes [OR03] and the two-dimensional stochastic particle dynamics introduced in [BF14].
2.1. Young diagrams. A partition is a nonincreasing integer sequence of the form λ = (λ 1 ≥ . . . ≥ λ (λ) > 0). The number of nonzero parts (λ) (which must be finite) is called the length of a partition. Partitions are represented by Young diagrams, such that λ 1 , λ 2 , . . . are lengths of the successive rows. The column lengths of a Young diagram are denoted by λ 1 ≥ λ 2 ≥ . . .. They form the transposed Young diagram λ . See Figure 4 for an illustration. . (2.1) If N < (λ), then s λ (u 1 , . . . , u N ) = 0 by definition. When all u i ≥ 0, the value s λ (u 1 , . . . , u N ) is also nonnegative. Along with evaluating Schur functions at finitely many variables, we also need their Plancherel specializations defined as This limit exists for every λ. It can be expressed through the number of standard Young tableaux of shape λ, which is the same as the dimension of the corresponding irreducible representation of the symmetric group of order λ 1 + λ 2 + . . .. The values s λ (Pl t ) are nonnegative for all t ≥ 0. When t = 0, we have s λ (Pl 0 ) = 1 λ=∅ . The Schur functions satisfy Cauchy summation identities. We will need the following version: where the sum runs over all Young diagrams. However, summands corresponding to (λ) > N vanish. There are also skew Schur symmetric functions s λ/µ which are defined through The function s λ/µ vanishes unless the Young diagram λ contains µ (notation: λ ⊃ µ). Skew Schur functions satisfy skew modifications of the Cauchy summation identity. They also admit Plancherel specializations, and, moreover, s λ/µ (Pl t ) is expressed through the number of standard tableaux of the skew shape λ/µ. We refer to, e.g., [Mac95, Ch I.5] for details.

Schur processes.
Here we recall the definition (at appropriate level of generality) of Schur processes introduced in [OR03]. Let ξ • be a speed function as in Section 1.2, and take a down- (Definition 1.2). A Schur process associated with this data is a probability distribution on sequences (λ; µ) of Young diagrams (see Figure 5 for an illustration) with probability weights Here ξ (a,b] for a ≤ b means the string (ξ a+1 , . . . , ξ b ). Note that some of the specializations above can be empty. The normalizing constant in (2.3) is which is computed using the skew Cauchy identity.
µ (4) Figure 5. An illustration of the Schur process (2.3) corresponding to a down-right path with r = 4. For convenience we take t 1 = N r = 0 so that the corresponding Young diagrams are almost sure empty under the Schur process.
The marginal distribution of any λ (i) under the Schur process (2.3) is a Schur measure [Oko01] whose probability weights are (2.4) 2.4. Correlation kernel. As shown in [OR03], the Schur process such as (2.3) can be interpreted as a determinantal random point process, and its correlation kernel is expressed as a double contour integral. To recall this result, consider the particle configuration corresponding to a sequence (2.2) (where we sum over all the µ (j) 's). The configurations λ (i) j − j, j ≥ 1, are infinite and are densely packed at −∞ (i.e., we append each λ (i) by infinitely many zeroes). Then for any m and any pairwise distinct locations (a p , x p ), p = 1, . . . , m, where 1 ≤ a p ≤ r − 1 and x p ∈ Z, we have P there are points of the configuration (2.5) at each of the locations (a p , x p ) = det [K SP (a p , x p ; a q , x q )] m p,q=1 . The kernel K SP has the form The integration contours in (2.6) are positively oriented simple closed curves around 0, the contour w in addition encircles {ξ x } x∈Z ≥1 , and on these contours |z| > |w| for i ≤ j and |z| < |w| for i > j.
2.5. Coupling PushTASEP and Schur processes. Fix a speed function ξ • as above. We will consider (half continuous Schur) random fields of Young diagrams satisfying the following properties: (1) (Schur field property) For any down-right path {(t i , N i )} r i=1 , the joint distribution of the Young diagrams λ (i) = λ (t i+1 ,N i ) is described by the Schur process corresponding to this down-right path. Note that this almost surely enforces the boundary conditions λ (0,N ) = λ (t,0) ≡ ∅, and also forces each diagram λ (t,N ) to have at most N parts.
(2) (PushTASEP coupling property) The collection of random variables is the length of the first column of λ (t,N ) ) has the same distribution as the values of the height function in the inhomogeneous PushTASEP having the speed function ξ • and started from the empty initial configuration. The first property states that a field couples together Schur processes with different parameters in a particular way, and the second property requires a field to possess additional structure relating it to the PushTASEP. The random field point of view was recently useful in [BBW18], [BM18], [BP17], [BMP19] in discovering and studying particle systems powered by generalizations of Schur processes.
The above two properties do not determine a field uniquely. In fact, there exist several constructions of fields satisfying these properties. They lead to different joint distributions of all the diagrams {λ (t,N ) }. However, due to the Schur field property, along down-right paths the joint distributions of diagrams are the same.
The oldest known such construction is based on the column RSK insertion. Connections between RSK, random Young diagrams, and stochastic particle systems can be traced [VK86], see also Another field coupling Schur processes and the PushTASEP was suggested in [BF14] based on an idea [DF90] of stitching together Markov processes connected by a Markov projection operator. A unified treatment of these two approaches was performed in [BP16], see also [OP13]. A variation of the field of [BF14] based on the Yang-Baxter equation was suggested recently in [BP17] (for Schur processes, as well as for their certain two-parameter generalizations), and further extended in [BMP19]. Since either of these approaches suffices for our purposes, let us outline the simplest one from [BF14].
Fix K ≥ 1, and consider the restriction of the field to the first K horizontal levels. Interpret t as continuous time, and the integers {λ (t,N ) i : 1 ≤ N ≤ K, 1 ≤ i ≤ N } as a two-dimensional timedependent array. 2 We will describe a Markov evolution of this array. Throughout the evolution, the integers will almost surely satisfy the following interlacing constraints: The array evolves as follows. Each of the integers at each level 1 ≤ N ≤ K has an independent exponential clock with rate ξ N . When the clock of λ (t,N ) j rings (almost surely, at most one clock can ring at a given time moment since the number of clocks is finite), its value is generically incremented by one. In addition, the following mechanisms are at play to preserve interlacing in the course of the evolution: for some m ≥ 1 before the increment of λ Thus described Markov processes are compatible for various K, and so they define a random field λ (t,N ) , t ∈ R ≥0 , N ∈ Z ≥1 . From [BF14] it follows that the collection of random Young diagrams {λ (t,N ) } satisfies the Schur field property, i.e., its distributions along down-right paths are given by Schur processes.
Proof of the PushTASEP coupling property. Let us now prove that the just constructed collection {λ (t,N ) } of Young diagrams satisfies the PushTASEP coupling property. Observe that N − λ 1 (t, N ) is the number of zeroes in the N -th row in the array (2.9). Due to interlacing, for each fixed t we can interpreth(t, N ) := N − λ 1 (t, N ) as the height function of a particle configurationx(t) = {x i (t)} i≥1 in Z ≥1 , with at most one particle per site. The initial condition isx i (0) = i, i ≥ 1. That is, we can determinex fromh using (1.2).
The time evolution of the particle configurationx(t) is recovered from the field λ (t,N ) . First, observe that any change inx can come only from the exponential clocks ringing at the rightmost zero elements of the interlacing array. There are two cases. Ifh(t, N ) =h(t, N − 1), then the rightmost clock at zero on level N corresponds to a blocked increment, which agrees with the fact thatx has no particle at location N . If, on the other hand,h(t, N ) =h(t, N − 1) + 1, then there is a particle inx at N which can jump to the right by one. This happens at rate ξ N . If this particle at N jumps and, moreover,h(t, N + 1) =h(t, N ) + 1, then the particle at N + 1 which is also present inx is pushed by one to the right, and so on. See Figure 6 for illustration. Figure 6. Left: in the interlacing array the framed zero is blocked and cannot increase. This corresponds to no particle inx(t) at N . Right: the circled zero at level N decides to increase at rate ξ N , and forces the circled zero at level N + 1 to increase, too. Inx(t) this corresponds to a jump of the particle at N which then pushes a particle at N + 1.
We see that the Markov processx(t) coincides with the PushTASEP in inhomogeneous space x(t) introduced in Section 1.2.
Remark 2.1. The field λ (t,N ) from [BF14] described above has another Markov projection to a particle system in Z which coincides with the PushTASEP with particle-dependent inhomogeneity. Namely, start the PushTASEP from the step initial configuration x i (t) = i, i ≥ 1, and let the particle x i have jump rate ξ i . The space is assumed homogeneous, so now variable jump rates are attached to particles. Then the joint distribution of the random variables {x i (t)} for all t ≥ 0, i ≥ 1, coincides with the joint distribution of {λ (t,i) 1 + i}. In particular, each x N (t) has the same distribution as λ 1 + N under the Schur measure ∝ s λ (ξ 1 , . . . , ξ N )s λ (Pl t ) (this is the same Schur measure as in (2.4)). Asymptotic behavior of PushTASEP with particle-dependent jump rates was studied in [BG13] by means of Rákos-Schütz type determinantal formulas [RS05], [BF08].
A third Markov projection of the field λ (t,N ) onto {λ (t,N ) N − N } N ≥1 recovers TASEP on Z with particle-dependent speeds. We refer to [BF14] for details on these other two Markov projections.
2.6. From coupling to determinantal structure. For any random field λ (t,N ) satisfying the Schur field property, the determinantal structure result of [OR03] recalled in Section 2.4 can be restated as follows: Theorem 2.2. For any m ∈ Z ≥1 and any collection of pairwise distinct locations (2.10) The integration contours are positively oriented simple closed curves around 0, the w contour additionally encircles {ξ x } x∈Z ≥1 , and the contours satisfy |z| > |w| for t ≤ s and |z| < |w| for t > s.
In particular, this theorem applies to the field from [BF14] recalled in Section 2.5 whose first columns are related to the PushTASEP as in (2.7)-(2.8).

2.7.
Kernel for column lengths. Let us restate Theorem 2.2 in terms of column lengths so that we can apply it to PushTASEP. The correlation kernel for the complement configuration is given by K := 1 − K F , where 1 is the identity operator whose kernel is the delta function. This follows from an observation of S. Kerov based on the inclusion-exclusion principle see [BOO00,Appendix A.3]. This leads to: First, recall Fredholm determinants on an abstract discrete space X. Let K(x, y), x, y ∈ X be a kernel on this space. We say that the Fredholm determinant of 1 + zK, z ∈ C, is an infinite series (2.11) One may view (2.11) as a formal series, but in our setting this series will converge numerically. Details on Fredholm determinants may be found in [Sim05] or [Bor10]. Fix a down-right path p = {(t i , N i )} r i=1 and consider the space ..,r as a determinantal process L p on X with kernel K (1.4) in the sense of Corollary 2.4.
Fix an arbitrary r-tuple y = (y 1 , . . . , y r ) ∈ Z r . We can interpret as the probability of the event that there are no points in the random point configuration L p in the subset X y := r i=1 {. . . , −y i − 2, −y i − 1, −y i } of X . This probability can be written (e.g., see [Sos00]) as the Fredholm determinant det(1 − χ y Kχ y ) X , (2.12) where χ y (x) = 1 x≤−y i for x ∈ X i is the indicator of X y ⊂ X viewed as a projection operator acting on functions. In particular, for r = 1 this implies Corollary 1.6 from the Introduction.
Remark 2.5. One can check that the sums in the Fredholm determinant (2.12) (as well as in (1.6) in the Introduction) are actually finite due to vanishing of K far to the left.

Asymptotic analysis
In this section we study asymptotic fluctuations of the random height function of the inhomogeneous PushTASEP at a single space-time point, and prove Theorems 1.12 and 1.13. We also establish more general results on approximating the kernel K (1.4) by the Airy kernel under weaker assumptions on ξ(·).
3.1. Rewriting the kernel. Let us rewrite K given by (1.4) to make the integration contours suitable for asymptotic analysis via steepest descent method.

(3.1)
Here the z contour in both integrals is a small positively oriented circle around 0, and the w contour in the double integral is a vertical line traversed downwards and located to the left of the z contour.
Proof. We start from formula (1.4) for the kernel. For t ≤ t (thus necessarily N ≥ N because we consider correlations only along down-right paths, cf. Definition 1.2) the z contour encircles the w contour. Note that the integrand does not have poles in z at the ξ a 's. Thus, exchanging for t ≥ t the z contour with the w contour at a cost of an additional residue, we see that the new contours in the double integral in (1.4) can be taken as follows: • the z contour is a small positive circle around 0; • the w contour is a large positive circle around 0 and {ξ a } a≥1 . The additional residue arising for t ≥ t is equal to the integral of the residue at z = w of the integrand over the single w contour. Because t ≤ t and N ≥ N , this residue does not have poles at the ξ a 's, and so the integration can be performed over a small contour around 0. Renaming w to z we arrive at the single integral in (3.1).
Finally, in the double integral the w integration contour can be replaced by a vertical line because: • the exponent e −t w ensures rapid decay of the absolute value of the integrand sufficiently far in the right half plane; • the polynomial factors w x +N z−w N b=1 (w − ξ b ) −1 for x < 0 ensure at least quadratic decay of the absolute value of the integrand for sufficiently large | Im w|.
This completes the proof.
Remark 3.2. The assumption x < 0 made in Proposition 3.1 agrees with the fact that we are looking at the leftmost points in the determinantal point process L p (Theorem 1.5), and these leftmost points almost surely belong to Z ≤0 . At the level of the PushTASEP this corresponds to h(t, N ) ≤ N . The event h(t, N ) = N (i.e., for which it would be x = 0) can be excluded, too, since it corresponds to no particles ≤ N jumping till time t. Since t goes to infinity, this is almost surely impossible. We thus assume that x < 0 throughout the text.
3.2. Critical points and estimation on contours. Rewrite the integrand in the double contour integral in (3.1) as where h := x + N , h := x + N , and the function S L has the form The signs inside logarithms are inserted for future convenience. The branches of the logarithms are assumed standard, i.e., they have cuts along the negative real axis. We apply the steepest descent approach (as outlined in, e.g., [Oko02, Section 3], in a stochastic probabilistic setting) to analyze the asymptotic behavior of the leftmost points of the determinantal process L p . To this end, we consider double critical points of S L which satisfy the following system of equations: (3.5) In the rest of this subsection we assume that N ≥ 1 and 0 < t < t e (N ) are fixed. Proof. Follows by monotonicity similarly to Lemma 1.9.
Denote the solution afforded by Lemma 3.4 by z L = z L (t, N ). Also denote by h L = h L (t, N ) the result of substitution of z L (t, N ) into the right-hand side of (3.5).
Lemma 3.5. The function z → S L (z; t, N, h L (t, N )) has a double critical point at z L (t, N ) which is its only critical point on the negative real half-line. All other critical points (of any order) of this function are real and positive.
Proof. The fact that z L (t, N ) is a double critical point of S L (·; t, N, h L (t, N )) follows from the above definitions. It remains to check that all other critical points of S L are real and positive. Let 0 < b 1 < . . . < b k be all of the distinct values of ξ 1 , . . . , ξ N . Then equation S L (z) = 0, that is, is equivalent to a polynomial equation of degree k + 1 with real coefficients. The right-hand side of (3.6) takes all values from −∞ to +∞ on each of the k − 1 intervals of the form (b i , b i+1 ). Therefore, (3.6) has at least k − 1 positive real roots. Since z L is a double root when h = h L , we have described at least k + 1 real roots to the equation S L (z) = 0, i.e., all of its roots. This completes the proof.
Keeping t, N fixed, plug h = h L (t, N ) into S L , and using (3.4)-(3.5) rewrite the result in terms of z L : Denote the expression inside the sum by R(z; ξ a ).
Lemma 3.6. On a circle through z L centered at the origin, Re S L (z; t, N, h L (t, N )) viewed as a function of z attains its maximum at z = z L .
for ϕ ∈ [0, π] (by symmetry, it suffices to consider only the upper half plane), and this derivative is equal to zero only for ϕ = 0. This implies the claim. Proof. For w = z L + ir, r > 0, we have (recall that z L < 0). This implies the claim.
We need one more statement on derivatives of the real part at the double critical point: Lemma 3.8. Along the w and z contours in Lemmas 3.6 and 3.7 the first three derivatives of Re S L vanish at z L , while the fourth derivative is nonzero.
Let us now deform the integration contours in the double contour integral in (3.1) so that they are as in Lemmas 3.6 and 3.7 (but locally do not intersect at z L ). We can perform this deformation without picking any residues in particular because the integrand is regular in z at all the ξ a 's. Lemmas 3.6 and 3.7 then imply that the asymptotic behavior of the integral for large L is determined by the contribution coming from the neighborhood of the double critical point z L . In Section 3.4 we make precise estimates.
3.3. Airy kernel. Before we proceed, let us recall the Airy kernel [TW93], [TW94] A(x; x ) := A ext (s, x; s, x ) = 1 (2πi) 2 x, x ∈ R. (3.7) In the contour integral expression, the v integration contour goes from e −i 2π 3 ∞ through 0 to e i 2π 3 ∞, and the u contour goes from e −i π 3 ∞ through 0 to e i π 3 ∞, and the integration contours do not intersect. The GUE Tracy-Widom distribution function is the following Fredholm determinant of (3.7): (3.8) Its expansion is defined analogously to (2.11) but with sums replaced with integrals over (r, +∞).
3.4. Approximation and convergence. Our first estimate is a standard approximation of the kernel K(t, N, x; t, N, x ) by the Airy kernel A (3.7) when both x, x are close to h L (t, N ) − N . Denote (3.9) In this subsection we assume that t = t(L) and N = N (L) depend on L such that for all sufficiently large L: • 0 < t < t e (N ) − cL for some c > 0; • for some m, M > 0 we have m < t(L) L < M and m < N (L) L < M . Lemma 3.9. Under our assumptions on (t(L), N (L)), as L → +∞ we have Proof. When (t, N ) = (t , N ), the indicator and the single contour integral in (3.1) cancel out, and so we have where h = x + N , h = x + N . Let the z and w integration contours pass near z L (without intersecting each other) and be as in Lemmas 3.6 and 3.7. We have For large L the main contribution to the double integral comes from a small neighborhood of the critical point z L = z L (t, N ). Indeed, fix a neighborhood of z L of size L −1/6 . By Lemma 3.8, if w or z or both are outside the neighborhood of size L −1/6 of z L , we can estimate Re(S L (z; t, N, h L ) − S L (w; t, N, h L )) < −cL 1/3 for some c > 0. This means that the contribution coming from outside the neighborhood of z L is asymptotically negligible compared to (−z L ) (h −h)L 1/3 in (3.10).
Inside the neighborhood of z L make a change of variables Herez,w are the scaled integration variables which are integrated over the contours in Figure 8. More precisely, |z|, |w| go up to order L 1/6 , and the contribution to the Airy kernel A coming from the parts of the contours in (3.7) outside this large neighborhood of zero is bounded from above by e −cL 1/2 for some c > 0, and so is asymptotically negligible. z w Figure 8. The integration contours forz andw in (3.12) leading to the Airy kernel approximation. Shaded are the regions where Re(z 3 ) < 0.
Using (3.12) and Taylor expanding as L → +∞ we have and similarly for the other term in the exponent in (3.11). Therefore, we have withz,w contours as in Figure 8. This completes the proof.
Lemma 3.10. Under our assumptions on (t(L), N (L)), let h = x + N and h = x + N be such that h − h L (t, N ) ≤ −sL 1/3 for some s > 0. Then for some C, c 1 , c 2 > 0 and L large enough we have Proof. First, observe that the assumptions imply that the double critical point |z L (t, N )| is uniformly bounded away from zero and infinity. Write the kernel as (3.11) with integration contours described in Lemmas 3.6 and 3.7 and locally around z L in the proof of Lemma 3.9. In the exponent we have If either z or w or both are outside of a L −1/6 -neighborhood of z L , we estimate (3.14) ≤ (h − h) log |z L | + (h − h L ) log |w/z L | − cL 1/3 (3.15) as in the proof of Lemma 3.9. The part in the exponent containing w is integrable over the vertical w contour, which leads to the first term in the estimate for |K|. Now, if both z, w are close to z L , make the change of variables (3.12) and write The part containingz,w is integrable over the scaled contours in Figure 8. Since the coefficient by Rew is positive and Rew ≤ −1 on our contour, we estimate this integral using the exponential integral ∞ 1 e Au du = e −A /A, where u corresponds to Rew, and A is the coefficient by Rew). This produces the second term in the estimate for |K|. This completes the proof. Fix s > 0 (to be taken large later) and separate the terms in the above Fredholm expansion where all x i > y − N − sL 1/3 , plus the remainder. In the former terms we use Lemma 3.9, and the latter terms are smaller due to Lemma 3.10.
When all x i > y − N − sL 1/3 ∼ h L − N + yL 1/3 − sL 1/3 , let us reparametrize the summation variables as x i = h L − N + u i L 1/3 , with u i ∈ R going from y − s to y (in increments of L −1/3 ). From Lemma 3.9 we have 3 n det i,j=1 [K(t, N, x i ; t, N, x j )] = 1 + O(L −1/3 ) and each n-fold sum over x i > h L − N + (y − s)L 1/3 can be approximated (within O(L −1/3 ) error) by the n-fold integral of the Airy kernel A from −y/d L to (s − y)/d L . Taking s sufficiently large and using the decay of the Airy kernel (e.g., see [TW94]) leads to the Tracy-Widom GUE distribution function at −y/d L .
Consider now the remaining terms. Using Lemma 3.10 we have (−z L ) x+N −x −N x <h L −N +(y−s)L 1/3 |K(t, N, x; t , N , x )| ≤ C 1 e −c 1 L 1/3 log L(s − y) + (s − y) −1 C 2 e −c 2 (s−y) 1 + O(L −1/3 ) for some C i , c i > 0. The first term decays rapidly for large L, and the second term can be made small for fixed y by choosing a sufficiently large s. Take the n-th term in (3.18) where some of the x i 's are summed from −∞ to y − N − sL 1/3 , and expand the n × n determinant along a column x j corresponding to x j < h L − N + (y − s)L 1/3 . The resulting n − 1 determinants are estimated via Hölder and Hadamard's inequalities. Thus, the remaining terms in the Fredholm expansion are negligible and can be included in the error in the right-hand side of (3.17). This completes the proof.
The final step of the proof of Theorem 1.13 (which would also imply Theorem 1.12) is to show that the approximation of the probability in Proposition 3.11 implies the convergence of the probabilities P(h(t, N ) > Lh(τ, η) + yL 1/3 ) to the GUE Tracy-Widom distribution function. This convergence would clearly follow if Observe that all sums in the definitions of z L , h L , d L in Section 3.2 are Riemann sums of the integrals from 0 to η from Section 1.5. The mesh of these integrals is of order L −1 , and due to the piecewise C 1 assumption on ξ(·), integrals are approximated by Riemann sums within O(L −1 ). This implies (3.19), which is the last step of in the proof of Theorem 1.13.
Combining the above formulas yields the hydrodynamic equation (1.9) for the limiting density.