Determinantal structures in space inhomogeneous dynamics on interlacing arrays

We introduce a space inhomogeneous generalization of the dynamics on interlacing arrays considered by Borodin and Ferrari. We show that for a certain class of initial conditions the point process associated to the dynamics has determinantal correlation functions and we calculate explicitly, in the form of a double contour integral, the correlation kernel for one of the most classical initial conditions, the densely packed. En route to proving this we obtain some results of independent interest on non-intersecting general pure-birth chains, that generalize the Charlier process, the discrete analogue of Dyson's Brownian motion. Finally, these dynamics provide a coupling between the inhomogeneous versions of the TAZRP and PushTASEP particle systems which appear as projections on the left and right edges of the array respectively.

The one that we will be concerned with in this contribution is due to Borodin and Ferrari 2 [7] (see also the independent related work of Warren and Windridge [53] and Warren's Brownian analogue of the dynamics [52]) based on some ideas from [23]. This is the simpler out of the three approaches to describe (see Definition 1.2 for a precise description) and in some sense, see [20], the one with "maximal noise". Many of the ideas and results from the important paper [7] have either directly generated or have been made use of in a very large body of work, see for example [12], [4], [17], [3], [15], [24], [16], [5], [49], [22], [10], [11], [8] for further developments and closely related problems. The other two approaches can be concisely described as follows: one of them (which is historically the first out of the three) is based on the combinatorial algorithm of the RSK correspondence, see [29], [40], [41], [42] and the other, which has begun to develop very recently, is based on the Yang-Baxter equation, see [20], [21]. Now, in the past few years, there has been considerable interest in constructing new integrable models in inhomogeneous space or adding spatial inhomogeneities, in a natural way, to existing models while preserving the integrability, see [18], [32], [26], [27], [2]. In this paper we do exactly that for the original (continuous time) dynamics of Borodin and Ferrari (see Definition 1.2). We show that for a certain class of initial conditions the point process associated to the dynamics has determinantal correlation functions. We then calculate explicitly, in the form of a double contour integral, the correlation kernel for one of the most classical initial conditions, the densely packed. This allows one to address questions regarding asymptotics and it would be interesting to return to this in future work. Here, our focus is on developing the stochastic integrability aspects of the model.
Finally, the projections on the edges of the interlacing array give two Markovian interacting particle systems of independent interest, see Remark 1.3 for more details. On the left edge we get the inhomogeneous TAZRP (totally asymmetric zero range process) or Boson particle system, see [4], [17], [51] and on the right edge we get the inhomogeneous PushTASEP, see [7], [6], [16], which is also studied in detail in the independent work of Leonid Petrov [46] which uses different methods 3 .
In the next subsection we give the necessary background in order to introduce the model and state our main results precisely.

Background and main result
We define the discrete Weyl chamber with non-negative coordinates: where Z + = {0, 1, 2, . . . }. We think of the coordinates x i as positions of particles and will use this terminology throughout, see Figure 1 for an illustration. We say that y ∈ W N and x ∈ W N+1 interlace 1 By this rather vague term we mean that there exist useful explicit formulae for the expectations of at least some observables. See also the introductions of [20] and [21] for a lengthier historical overview and comparison between the three approaches. 2 Thus the name "Borodin-Ferrari push-block dynamics" that we use throughout the paper. 3 Petrov finds an integrable structure underlying PushTASEP by making a connection to the theory of Schur processes, see [43], [44]. He also performs some asymptotics. The overlap in terms of results and techniques between the two papers is minuscule. and write y ≺ x if: We define the set of Gelfand-Tsetlin patterns (interlacing arrays) of length N by 4 : The basic data in this paper is a rate function which we think of as the spatial inhomogeneity of the environment. It governs how fast or slow particles jump when at a certain position. We enforce the following assumption throughout the paper. Definition 1.1. Assumption (UB). We assume that the rate function λ : Z + → (0, ∞) is uniformly bounded away from 0 and ∞: We now introduce the inhomogeneous space push-block dynamics in GT N . This is the continuous time Markov jump process in GT N described as follows: Definition 1.2. Borodin-Ferrari inhomogeneous space push-block dynamics. Let λ(·) satisfy (UB). Let M N be the initial distribution (possibly deterministic) of particles on GT N . We now describe Markov dynamics in GT N denoted by (X N (t; M N ); t ≥ 0) = X 1 (t), . . . , X N (t) ; t ≥ 0 where the projection on the k th level is given by X k (t); t ≥ 0 = X k 1 (t), X k 2 (t), . . . , X k k (t) ; t ≥ 0 . Each particle has an independent exponential clock of rate λ(⋆) depending on its current position ⋆ ∈ Z + for jumping to the right by one to site ⋆ + 1. The particles interact as follows, see Figure 1 for an illustration: if the clock of particle X n k rings first then it will attempt to jump to the right by one.
• (Blocking) In case X n−1 k = X n k the jump is blocked (since a move would break the interlacing; lower level particles can be thought of as heavier).
• (Pushing) Otherwise it moves by one to the right, possibly triggering instantaneously some pushing moves to maintain the interlacing. Namely, if the interlacing is no longer preserved with the particle labelled X n+1 k+1 then, X n+1 k+1 also moves instantaneously to the right by one and this pushing is propagated (instantaneously) to higher levels, if needed. Remark 1.3. Inhomogeneous Boson and PushTASEP. It is easy to see that the particle systems at the left X 1 1 (t), X 2 1 (t), . . . , X N 1 (t) ; t ≥ 0 and right X 1 1 (t), X 2 2 (t), . . . , X N N (t) ; t ≥ 0 edge respectively in the GT N -valued dynamics of Definition 1.2 enjoy an autonomous Markovian evolution.
The left edge process is called the inhomogeneous TAZRP (totally asymmetric zero range process) or Boson particle system, see [4], [17], [51]. In particular, in [51] a contour integral Figure 1: A configuration of particles in GT 4 . If the clock of the particle labelled X 3 1 rings, which happens at rate λ(⋆ + 1), then the move is blocked since interlacing with X 2 1 would be violated. On the other hand, if the clock of the particle X 2 2 rings, which happens at rate λ(⋆ + 3), then it jumps to the right by one and instantaneously pushes both X 3 3 and X 4 4 to the right by one as well, for otherwise the interlacing would break.
expression (for a q-deformation of the model) is obtained for its transition probabilities. The fact that, as we shall also see in the sequel, the distribution of particles at fixed time t ≥ 0 is a marginal of a determinantal point process (with explicit kernel) is essentially contained in the results of [32] (which makes use of different methods).
The right edge particle system is called inhomogeneous PushTASEP, see [7], [6] and also [16] for a q-deformation of the homogeneous case. The fact that this particle system has an underlying determinantal structure is new, but is also obtained in the independent work of Petrov [46] that uses different methods.

Remark 1.4.
A generalization of the model. Borodin and Ferrari in fact considered a more general model where all particles at level k jump at rate β k (independent of their position). Using the techniques developed in this paper we can study the following model that allows for level inhomogeneities as well 5 . Let {α i } i≥1 be a sequence of numbers such that: The dynamics are as in Definition 1.2 with the modification that each of the particles at level k jumps at rate α k + λ(⋆) depending on its position ⋆. Since our main motivation in this work is the introduction of spatial inhomogeneities we will only consider the level homogeneous case of Definition 1.2 in detail. However, in the sequence of Remarks 2.21,2.24,3.3,3.10 we indicate the essential modifications required at each stage of the argument to study the more general model.
Observe that, for any n ≤ N the process described in Definition 1.2 restricted to GT n is an autonomous Markov process. We consider the natural projections: 5 We believe (but do not have a rigorous argument) that this is the most general model of continuous time Borodin-Ferrari dynamics (namely the particle interactions being as in Definition 1.2) with determinantal correlations and which can be treated with the methods developed here. It is plausible however that if one considers discrete time dynamics that there exists an even more general model in inhomogeneous space involving geometric or Bernoulli jumps, as is the case in the homogeneous setting, see [7]. forgetting the top row x (N+1) and we write: for the corresponding projective limit, consisting of infinite Gelfand-Tsetlin patterns. We Suppose we are given such a consistent sequence of distributions {M N } N≥1 . Then, by construction since the projections on any sub-pattern are autonomous, the processes (X N (t; M N ) ; t ≥ 0) N≥1 are consistent as well: and we can correctly define (X ∞ (t; {M N } N≥1 ) ; t ≥ 0), the corresponding process on GT ∞ . Now, we will be mainly concerned with the so-called densely packed initial conditions {M dp N } N≥1 defined as follows: M dp N = δ (0≺(0,1)≺(0,1,2)≺···≺(0,1,...,N−1)) . Clearly {M dp N } N≥1 is a consistent sequence of distributions. We simply write (X ∞ (t); t ≥ 0) for the corresponding process on GT ∞ .
Observe that, X ∞ (t) for any t ≥ 0 gives rise to a random point process on N × Z + which we denote by P t ∞ . We will use the notation z = (n, x) to denote the location of a particle, with n being the level/height/vertical position while x is the horizontal position. Finally, it will be convenient to introduce the following functions, which will play a key role in the subsequent analysis. Definition 1.5. For x ∈ Z + , we define: Remark 1.6. Clearly ψ x (w; λ) = 1/p x+1 (w; λ) but it is preferable to think of them as two distinct families of functions. Observe that p x (w) is a polynomial of degree x and p x (0) = 1.
We have then arrived at our main result. Theorem 1.7. Let λ : Z + → (0, ∞) satisfy (UB). Consider the point process P t ∞ on N × Z + obtained from running the dynamics of Definition 1.2 for time t ≥ 0 starting from the densely packed initial condition. Then for all t ≥ 0, P t ∞ has determinantal correlation functions. Namely, for any k ≥ 1 and distinct points z 1 = (n 1 , x 1 ), . . . , z k = (n k , x k ) ∈ N × Z + : where the correlation kernel K t (·, ·; ·, ·) is explicitly given by: where C λ is a counter clockwise contour encircling 0 and the points {λ(x)} x≥0 while C 0 is a small counter clockwise contour around 0 as in Figure 2. Remark 1.8. In the homogeneous case, λ(·) ≡ 1, the correlation kernel in Theorem 1.7 above reduces to a kernel from [7] (there is a number of different kernels in [7] which give rise to equivalent determinantal point processes). To see this, it is most convenient to compare the expression for the kernel in Proposition 4.2 of [7] with the formulae from Section 3 herein (in particular see the expressions in display (4.4) in [7] and Lemmas 3.4, 3.5 and display (35) in this paper).

Intermediate results and strategy of proof
The proof of Theorem 1.7 essentially splits into two parts that we now elaborate on.
Firstly showing that under certain initial conditions M N on GT N that we call Gibbs, of which the densely packed is a special case, the distribution of X N (t; M N ) for any t ≥ 0 is explicit (and again of Gibbs form), see Proposition 2.22. By applying an extension of the famous Eynard-Mehta theorem [28] to interlacing particle systems, see [19], [9], it is then fairly standard, see Section 3.1, that under some (rather general) Gibbs initial conditions the point process associated to X N (t; M N ) has determinantal correlation functions (with a not yet explicit kernel).
A key ingredient in the argument for this first part is played by a remarkable block determinant kernel given in Definition 2.5 that we call the two-level Karlin-McGregor formula (for the original, single level, Karlin-McGregor formula see [31]). This terminology is due to the fact that, as we see in Section 2 and in particular Theorem 2.20, this provides a coupling between two Karlin-McGregor semigroups so that the corresponding processes interlace.
An instance of such a formula, in the context of Brownian motions, was first discovered by Warren in [52]. Later it was understood that it can be further developed to include general one-dimensional diffusions and birth and death chains and this was achieved in [1] and [2] respectively. Part of our motivation for this paper was to enlarge the class of models for which such an underlying structure has been shown to exist 6 .
We note that, a crucial role in the previous works [1], [2] was played by the notion of Siegmund duality, see [48] (also [34], [35] for other uses of duality in integrable probability) alongside with reversibility and their absence in the present setting is an important conceptual difference 7 . However, as it turns out an analogue of this duality suitable for 6 Currently this class has been shown to include essentially all examples of Borodin-Ferrari type dynamics in continuous time with determinantal correlations studied in the literature (both in continuous and discrete space, see [1] and [2] respectively). We expect such a formula to exist in the fully discrete setting (both time and space) as well and we plan to investigate this in future work. In this case, the Borodin-Ferrari dynamics give rise to shuffling algorithms for sampling random boxed plane partitions [10], [11] or domino tilings [38], [8]. 7 In particular the statements of the results look different to the ones in [1] and [2]. More precisely, in [1] and [2] we obtain couplings between a Karlin-McGregor semigroup associated to a diffusion (or birth and death chain) and the one associated to (a Doob transform of) its Siegmund dual (which is in general different to the original process, except for a smaller sub-class of self-dual ones). In the present paper the couplings are between two Karlin-McGregor semigroup associated to a general pure-birth chain (and Doob transforms thereof), see Section 2 for more details. our purposes does exist and is proven in Lemma 2.2.
Moreover in Section 2 we prove some results on general non-intersecting pure-birth chains that generalize the Charlier process [33], [7], the discrete analogue of Dyson's Brownian motion [25]. In particular, we construct harmonic functions (and also more general eigenfunctions) that are given as determinants with explicit entries.
Finally, we should mention that all of these formulae have in some sense their origin in the study of coalescing stochastic flows of diffusions and Markov chains, see [37]. Recently a general abstract theory has been developed for constructing couplings between intertwined Markov processes based on random maps and coalescing flows, see [36]. It would be interesting to understand to what extent this is related to the present constructions and more generally to analogous couplings in integrable probability.
The second part of the proof involves the solution of a certain biorthogonalization problem which gives the explicit form of the correlation kernel K t and is performed in Section 3. An important difference to the works [7], [12], [9], [13] where corresponding biorthogonalization problems were analysed is that the functions involved in the current problem are not translationally invariant in the spatial variable. This is where we make use of the functions {ψ x } x∈Z + and {p x } x∈Z + which arise in the spectral theory of a general pure-birth chain (see display (8) for the spectral expansion of the transition kernel of the chain in terms of them). These provide both intuition and also make most of the (otherwise quite tedious) computations rather neat, see in particular the sequence of Lemmas 3.4, 3.5, 3.6, 3.7, 3.8 and their proofs.
Acknowledgements. I would like to thank Maurice Duits, Patrik Ferrari and Jon Warren for some early discussions which motivated the problems studied in this paper. Moreover, I am very grateful to Leonid Petrov for sending me his preprint and for some interesting questions and remarks. I would also like to thank Alexei Borodin for some very interesting suggestions and pointers to the literature. Finally, I am very grateful to an anonymous referee for a careful reading of the paper and some very useful suggestions and remarks which have improved the presentation. The research described here was supported by ERC Advanced Grant 740900 (LogCorRM).

The one dimensional chain
We define the forward and backwards discrete derivatives We define the following pure-birth chain which is the basic building block of our construction: this is a continuous time Markov chain on Z + which when at site x jumps to site x + 1 with rate λ(x). Its generator is then given by (the subscript indicates the variable on which it is acting): Due to (UB) non-explosiveness and thus existence and uniqueness of the pure-birth process is immediate (simply compare with a Poisson process with constant rate M).
We write e tL ; t ≥ 0 for the corresponding Markov semigroup and abusing notation we also write e tL (x, y) for its transition density, namely the probability a Markov chain with generator L goes from site x to y in time t. This is the unique solution to both the Kolmogorov backwards and forwards equations, see Section 2.6 of [39]. The backwards equation, which we make use of here, reads as follows, for t > 0, x, y ∈ Z + : It is an elementary probabilistic argument that e tL (x, y) is explicit, see Section 2 of [51] (or alternatively simply check that the expression below solves the Kolmogorov equation).
We have the following spectral expansion for it, see display (2.1a) of [51]:

Remark 2.1.
Here we can simply pick any counter clockwise contour which encircles the points {λ(x)} x≥0 and not necessarily 0 as well. The fact that C λ encircles 0 will be useful in the computation of the correlation kernel later on.
Throughout the paper we use the notation 1 [[0,y]] (·) for the indicator function of the set {0, 1, 2, . . . , y}. We have the following key relation for the transition density e tL (x, y).
Proof. We let and show that this solves the Kolmogorov backwards equation. By uniqueness the statement follows. The t = 0 initial condition follows from: Finally, as required.

Remark 2.3.
It is also possible to prove Lemma 2.2 using the explicit spectral expansion (8), the argument is similar to the one in the proof of Lemma 3.4.

Two level dynamics
We begin with a classical definition due to Karlin . This semigroup has the following probabilistic interpretation: it corresponds to N independent copies of a chain with generator L killed when they intersect, see [31], [30]. A conditioned upon non-intersection version of this semigroup will govern the dynamics on projections on single levels of our interlacing arrays, see Theorem 2.20 and Proposition 2.22 below.
We now move on to study two-level dynamics. We first consider the space of pairs (y, x) which interlace: We have the following key definition that we call the two-level Karlin-McGregor formula, since as we shall see in the sequel this provides a coupling for two Karlin-McGregor semigroups, so that the corresponding processes interlace.
where the matrices and N × N respectively are given by: Observe that, the equivalence of the two representations for A t is by definition while for D t is due to Lemma 2.2.
We begin our study of U t by proving some of its basic properties: Proof. We prove positivity first. A direct verification from Definition 2.5 appears to be hard (although it would be interesting to have one). Instead, we give a simple probabilistic argument.
Let ((S t (x 1 ), . . . , S t (x m )); t ≥ 0) be a system of m independent chains with generator L, starting from (x 1 , . . . , x m ) ∈ W m , and which coalesce and move together once any two of them meet. We denote their law by P. Let z, z ′ ∈ W m and t ≥ 0. Then, we have: The claim is a consequence of the Karlin-McGregor formula, see Proposition 2.5 in [2] for a proof in a completely analogous setting. Observe that (10) allows to give the following probabilistic representation for U N,N+1 t (by applying discrete derivatives to the RHS of (10) we match with the expression from Definition 2.5): is then a consequence of the fact that the events are increasing both as the variables y i decrease and as the variables x ′ j increase. Finally, we need to prove that for any (y, x) ∈ W N,N+1 and t ≥ 0: We claim that for any (y, x) ∈ W N,N+1 and t ≥ 0: Since P N t ; t ≥ 0 is sub-Markov the statement of the proposition follows. Now in order to prove the claim we take the sum {x ′ :(y ′ ,x ′ )∈W N,N+1 } in the explicit form of the kernel from Definition 2.5 and use multilinearity of the determinant. Then, the claim follows from the relations below: and simple row-column operations.
We now introduce the inhomogeneous two-level dynamics: Definition 2.7. Two level inhomogeneous push-block dynamics. This is the continuous time Markov chain ((Y(t), X(t)) ; t ≥ 0) in W N,N+1 , with possibly finite lifetime, described as follows. Each of the 2N + 1 particles evolves as an independent chain with generator L subject to the following interactions. The Y-particles evolve autonomously. When a potential move by the X-particles would break the interlacing it is blocked, see Figure 3. While if a potential move by the Y-particles would break the interlacing then the corresponding X-particle is pushed to the right by one, see Figure 4. The Markov chain is killed when two Y-particles collide, at the stopping time: Figure 3: (Blocking) A jump of x i is blocked by y i so that the interlacing is maintained.
Here, the clock of x i rings with rate λ(⋆).
x i+1 t is a member of by Lemma 2.6) follows by a generic argument presented in a completely analogous setting in Section 3 in [2], see also [15].
First, observe that we have the t = 0 initial condition: This follows directly from the form of U N,N+1 t (y, x), (y ′ , x ′ ) , by noting that as t ↓ 0, the diagonal entries converge to 1 x i = x ′ i , 1 y i = y ′ i , while all other contributions to the determinant vanish.
Moreover, observe that we have the Dirichlet boundary conditions when two Ycoordinates coincide: Moving on, note that (here we are abusing notation slightly by using the same notation for both the matrices and their scalar entries) we have the following, for any x, x ′ ∈ Z + fixed and t > 0: To see the relation for C t observe that: We define int W N,N+1 to be the set of all pairs (y, x) ∈ W N,N+1 which when any of the x or y coordinates is increased by 1 they still belong to W N,N+1 ; namely the pairs (y, x) ∈ W N,N+1 so that (y, x) + (e i , 0), (y, x) + (0, e i ) ∈ W N,N+1 with e i being the unit vector in the i-th coordinate. Observe that in int W N,N+1 each of coordinates evolves as an independent chain with generator L which do not interact. Then, by the multilinearity of the determinant and relations (12) and (13) we obtain: It remains to deal with the interactions. We will only consider one blocking and one pushing case, as all others are entirely analogous. First, the blocking case with x 1 = y 1 = x. In order to ease notation and also make the gist of the simple argument transparent we further restrict our attention to the rows containing x 1 , y 1 . In fact, it is not hard to see that it suffices to consider the 2 × 2 matrix determinant given by, with x ′ , y ′ ∈ Z fixed: By taking the d dt -derivative of the determinant, we obtain using (12) and (13): On the other hand, what we would like to have according to the dynamics in Definition 2.7 is simply the following: We thus must show that: which corresponds to particle x 1 being blocked when x 1 = y 1 and x 1 tries to jump (see the configuration in Figure 3). In order to obtain (14) we shall work on the RHS. We multiply the second row by −λ(x) −1 and add it to the first row to obtain and analogously for the second column, which then gives us the LHS of (14) as desired.
Similarly, we consider a pushing move with y 1 = x, x 2 = x + 1 and x ′ , y ′ ∈ Z + fixed (see the configuration in Figure 4): We calculate using the relations (12) and (13): From the dynamics in Definition 2.7 we need to have the following: Hence, we need to show: which follows from (14) after relabelling x → x + 1.
We now need a couple of definitions whose purpose will be clear shortly.

Definition 2.9.
We define the positive kernel Λ N+1 N from W N+1 to W N by its density (with respect to counting measure): Abusing notation, we can also view Λ N+1 N as a kernel from W N+1 to W N,N+1 , in which case we write Λ N+1 N (x, (y, z)) = Λ N+1 N (x, y). Observe that, this is supported on elements (y, z) ∈ W N,N+1 such that z ≡ x.

Definition 2.10. For any N ≥ 1 and x ∈ W N define the functions h
The functions h N can in fact be written as determinants whose entries are defined recursively: Lemma 2.11. Let N ≥ 1 and x ∈ W N . Then, where the functions I i are defined by, for x ∈ Z + : Proof. Direct computation by induction using multinearity of the determinant.
Remark 2.13. Lemma 2.11 implies that the sequence of functions {I i (·; λ)} i≥1 forms a (discrete) extended complete Chebyshev system on Z + , see [30]. On the real line and under certain assumptions, such systems have been classified and are characterized through a recurrence like (15) (with integrals instead of sums), see [30].

Remark 2.14. It is possible to express the entries of the determinant representation for h N (x; λ) in terms of contour integrals as we shall see in Section 3. This is essential in order to perform the computation of the correlation kernel.
Proof. This computation is implicit in the proof of Lemma 2.6.
Similarly we have: Proposition 2.16. For t ≥ 0, we have the equalities of positive kernels, Proof. We take the sum {y:(y,x)∈W N,N+1 } in the explicit form of the kernels and use multilinearity of the determinant. Then, the statement follows from the relations below and simple row-column operations.  In order to proceed we require a general abstract definition. For a possibly sub-Markov semigroup (P(t); t ≥ 0) having a strictly positive eigenfunction h with eigenvalue e ct (i.e. P(t)h = e ct h) we define its Doob h-transform by e −ct h −1 • P(t) • h; t ≥ 0 . We note that this is an honest Markovian semigroup. Thus, Proposition 2.18 allows us to correctly define the Doob h-transformed versions of the semigroups and kernels above: Note that, by their very definition, all of these are now Markovian. Moreover, as we have done previously, we can also view L N+1 N as a Markov kernel from W N to W N,N+1 . We observe that for the distinguished special case λ(·) ≡ 1, P N t ; t ≥ 0 is the semigroup of the well-known Charlier process, see [33], [41], [40], [7] the discrete analogue of Dyson's Brownian motion [25]. With all these preliminaries in place we have: Proposition 2.19. For t ≥ 0, we have the intertwining relations between Markov semigroups: Proof. These relations are straightforward consequences of Propositions 2.15, 2.16 and 2.17 respectively.
Observe that, the h-transform by h (N,N+1) (y, x) = h N (y) conditions the Y-particles to never collide and the process with semigroup U N,N+1 t ; t ≥ 0 has infinite lifetime. Under this change of measure the evolution of the Y-particles is autonomous with semigroup P N t ; t ≥ 0 , while the X-particles evolve as N + 1 independent chains with generator L interacting with the Y-particles through the same push-block dynamics of Definition 2.7. We now arrive at the main result of this section. N+1 . Then, the projection on the X-particles is distributed as a Markov process with semigroup P N+1 t ; t ≥ 0 and initial condition M N+1 . Moreover, for any fixed time T ≥ 0, the conditional distribution of (X(T), Y(T)) given X(T) satisfies:

Theorem 2.20. Consider a Markov process ((Y(t), X(t))
Proof. Let S be the operator induced by the projection on the x-coordinates: (we do not indicate dependence on N). Observe that, Then, the first statement of the theorem, by virtue of the intertwining relation (23), is an application of the theory of Markov functions due to Rogers and Pitman, see Theorem 2 in [47] (applied to the function s above). Finally, for the conditional law statement (25) see Remark (ii) following Theorem 2 of [47].
We observe that the strictly positive eigenfunction h Doob h-transforms L 2 to L 1 : We define, for n ≥ 1: By an inductive argument, making use of Proposition 2.17, we can show that h Then, all of the results above have natural extensions involving these quantities (whose precise statements we omit) to the level inhomogeneous setting.

Consistent multilevel dynamics
We have the following multilevel extension of the results of the preceding subsection. Proposition 2.22. Let P k t ; t ≥ 0 and L k k−1 denote the semigroups and Markov kernels defined in (20) and (19) above and let M N (·) be a probability measure on W N . Define the following Gibbs probability measure M N on GT N with density: Consider the process (X N (t; M N ) ; t ≥ 0) = X 1 (t) , X 2 (t) , . . . , X N (t) ; t ≥ 0 in Definition 1.2.
Then, for 1 ≤ k ≤ N the projection on the k-th level X k (t); t ≥ 0 is distributed as a Markov process evolving according to P k t ; t ≥ 0 . Moreover, for any fixed T ≥ 0, the law of X 1 (T), . . . , X N (T) is given by the evolved Gibbs measure on GT N : Proof. The proof is by induction. For N = 2, this is Theorem 2.20. Assume the result is true for N − 1 and we prove it for N. We first observe that the induced measure on GT N−1 is again Gibbs. Then, from the induction hypothesis X N−1 (t); t ≥ 0 is a Markov process with semigroup P N−1 t ; t ≥ 0 . Moreover, the joint dynamics of X N−1 (t), X N (t); t ≥ 0 are those considered in Theorem 2.20 (with semigroup U N−1,N t ) and thus by the aforementioned result, we obtain that X N (t); t ≥ 0 is distributed as a Markov process with semigroup P N t ; t ≥ 0 . Furthermore by the same theorem we have that, for fixed T ≥ 0, the conditional law of X N−1 (T) given X N (T) is L N N−1 X N (T), · . Hence, since the distribution of X N (T) has density M N P N T (·), we get by the induction hypothesis, that the fixed time T ≥ 0, distribution of X 1 (T), . . . , X N (T) is given by (27) as desired.
We observe that the densely packed initial condition M dp N is clearly Gibbs. We close this subsection with a couple of remarks on generalizations of this result.

Inhomogeneous Gelfand-Tsetlin graph and Plancherel measure
This subsection is independent to the rest of the paper and can be skipped. However, it provides some further insight into the constructions of the present work and how they fit into a wider framework. We begin with some notation. Let denote the discrete chamber without the non-negativity restriction. The definitions of interlacing in this setting and ofW N,N+1 are also completely analogous (we simply drop non-negativity).
Definition 2.25. We consider a graded graph Γ = Γ λ with vertex set ⊎ N≥1W N . Two vertices x ∈W N+1 and y ∈W N are connected by an edge if and only if they interlace. For all N ≥ 1 we assign a weight/multiplicity, denoted by mult λ (y, x), to each edge (y, x) ∈W N,N+1 , and more generally to all pairs (y, x) ∈W N ×W N+1 : The distinguished case Γ 1 with λ(·) ≡ 1 is the Gelfand-Tsetlin graph 8 , see [50], [14]. This describes the branching of irreducible representations of the chain of unitary groups, see [50], [14]. We propose to call the more general case Γ λ defined above the inhomogeneous Gelfand-Tsetlin graph (in fact it is a family of graphs, one for each function λ).
We now define the dimension dim λ N (x) of a vertex x ∈W N , inductively by: We can associate a family of Markov kernels {Λ N+1→N } N≥1 fromW N+1 toW N to the graph Γ λ given by: Observe that, by the very definitions, when restricting to the positive chambers W N (namely considering the subgraph Γ + λ = ⊎ N≥1 W N ) we have: We say that a sequence of probability The extremal points of the convex set of consistent probability measures form the boundary of the graph Γ λ . In the homogeneous case λ(·) ≡ 1, the boundary of the Gelfand-Tsetlin graph Γ 1 has been determined explicitly and is in bijection (see [50], [14], [45] for more details and precise statements) with the infinite dimensional space Ω: and we also write: The extremal consistent sequence of probability measures M N γ + N≥1 corresponding to γ + ≥ 0 with all the other parameters on Ω identically equal to zero is called the Plancherel measure 9 for the infinite dimensional unitary group, see [13]. The connection to the present paper is through the following, see [7] where the right hand side is defined for λ(·) ≡ 1. Now, due to observation (28) and the fact that P N γ + ((0, 1, . . . , N − 1), ·) is supported on W N the following is an immediate consequence of the intertwining relation (24)  Thus, the sequence P N γ + ((0, 1, . . . , N − 1), ·)

N≥1
can be viewed as the analogue of the Plancherel measure for the more general graphs Γ λ . It would be interesting to understand whether this sequence is actually extremal for Γ λ for general λ. A more ambitious question would be whether there exists a complete classification of extremal consistent measures for Γ λ , in analogy to the case of the Gelfand-Tsetlin graph Γ 1 .

Remark 2.27. Analogous constructions exist for the level-inhomogeneous generalization of the
Gelfand-Tsetlin graph c.f. Remarks 1.4,2.21,2.24. 3 Determinantal structure and computation of the kernel

Eynard-Mehta Theorem and determinantal correlations
We will make use of one of the many variants of the famous Eynard-Mehta Theorem [28], and in particular a generalization to measures on interlacing particle systems, see [19], [9]. More precisely, we will use Lemma 3.4 of [9]. For the convenience of the reader and to set up some notation we reproduce it here: where x n n+1 are some "virtual" variables, which we also denote by virt, and Z N is a non-zero normalization constant. Then, the correlation functions are determinantal. To write down the kernel we need some notation. Define, where (a * b)(x, y) = z∈Z a(x, z)b(z, y). Also, define for 1 ≤ n < N: are linearly independent and generate the n-dimensional space V n . For each 1 ≤ n ≤ N we define a set of functions {Φ n j (x), j = 0, . . . , n − 1} determined by the following two properties: • For 1 ≤ i, j ≤ n − 1 we have: Finally, assume that φ n (x n n+1 , x) = c n Φ (n+1) 0 (x) for some c n 0, n = 1, . . . , N − 1. Then, the kernel takes the simple form: From Proposition 2.22 we get that Law X N (t; M dp N ) for fixed time t ≥ 0 is given by, where we use the notation △ N = (0, 1, . . . , N − 1): Using the spectral expansion (8) of e tL (x, y) and row operations (recall that p x (w) is a polynomial of degree x in w) we can rewrite the display above as follows, for a (different) non-zero constant Z N : where the functions {Ψ N N−i (·)} N i=1 (we suppress dependence on the time variable t since it is fixed) are given by: Moreover, we note that it is possible (see for example [52]) to write the indicator function for interlacing as a determinant, for y ∈ W n−1 , x ∈ W n : Thus, we can write the measure above in a form that is within the scope of Proposition 3.1: In particular this implies determinantal correlations. The explicit computation of the kernel is performed in the next section.

Computation of the correlation kernel
Our aim now is to solve the biorthogonalization problem given in Proposition 3.1 and obtain concise contour integral expressions for the families of functions appearing therein. This is achieved in the following sequence of lemmas. Firstly in order to ease notation, since we are in the level homogeneous case, we define: φ (n 2 −n 1 ) (z 1 , z 2 ) = φ (n 1 ,n 2 ) (z 1 , z 2 ).
Lemma 3.5. For 1 ≤ k ≤ N, we have: Lemma 3.6. For 1 ≤ k ≤ N, we have: We define a family of functions Φ · · (·) on Z + , for 1 ≤ n ≤ N (again we suppress dependence on the variable t since it is fixed): Lemma 3.7. The functions Φ are biorthogonal to the Ψ's. More precisely, for any Finally, it is clear that we have the following. Assuming the auxiliary results above we first prove Theorem 1.7. The proofs of these results are given afterwards.
Proof of Theorem 1.7. Making use of Proposition 3.1 (by virtue of the preceding auxiliary lemmas) we get that: The term φ (n 2 −n 1 ) (x 1 , x 2 ) is given by, from Lemma 3.5: It then remains to simplify the sum: Proof of Lemma 3.4. It suffices to show that: It is in fact equivalent to prove that: For any R > 1, we consider a counter clockwise contour C ≥R λ that contains 0 and {λ(x)} x≥0 and for which the following uniform bound holds: Such a contour exists because of assumption (UB); we can simply take a very large circle. Clearly, the left hand side of display (37) is equal to (since we can deform the contour C λ to C ≥R λ without crossing any poles): We now claim that uniformly for w ∈ C ≥R λ : Assuming this, display (37) immediately follows and thus also the statement of the lemma. Now, after a simple relabelling (more precisely by writing a i = λ(y + i + 1)) in order to establish the claim it suffices to prove the following result. Let {a i } i≥0 be a sequence of numbers in [s, M] and let C ≥R a be the contour defined above. Then, uniformly for w ∈ C ≥R a we have: For the inductive step, we first compute: We can then deform the contour C λ to C ≥R λ as in the proof of Lemma 3.4 and use (39) to conclude.
Proof of Lemma 3.6. We make use of Lemma 3.5 and apply the same arguments as in the proof of Lemma 3.4 to compute: Proof of Lemma 3.7. We first write using the explicit expression: We claim that for any l ∈ Z + , uniformly for u ∈ K, where K is an arbitrary compact neighbourhood of the origin, we have: Then, (40) becomes 1 2πı C 0 e tu 1 u j+1 e −tu u i du = 1(i = j).
In order to establish the claim, it is equivalent (by taking finite linear combinations, since p x (·) is a polynomial of degree x) to prove that for any l ∈ Z + , uniformly for u ∈ K (an arbitrary compact neighbourhood of the origin): Now, observe that for fixed u, the function x → p x (u) is an eigenfunction, with eigenvalue −u, of the generator L: Thus, for u fixed we have: e tL p · (u) (l) = e −tu p l (u), l ∈ Z + . Now, in order to show that the convergence is uniform for u ∈ K we proceed as follows. We first estimate: On the other hand, making use of the spectral expansion (8), for any l ∈ Z + , we have the following bound, where R can be picked arbitrarily large: Here C denotes a generic constant independent of x (we suppress dependence of C on l). We pick R large enough so that: Then using the Weirstrass M-test, since for all x ∈ Z + sup u∈K |e tL (l, x)p x (u)| ≤ Cη x , where η < 1, we get that, for any l ∈ Z + : x≥0 e tL (l, x)p x (u) = e −tu p l (u), uniformly for u on compacts K, as required.