Hidden diagonal integrability of $q$-Hahn vertex model and Beta polymer model

We study a new integrable probabilistic system, defined in terms of a stochastic colored vertex model on a square lattice. The main distinctive feature of our model is a new family of parameters attached to diagonals rather than to rows or columns, like in other similar models. Because of these new parameters the previously known results about vertex models cannot be directly applied, but nevertheless the integrability remains, and we prove explicit integral expressions for $q$-deformed moments of the (colored) height functions of the model. Following known techniques our model can be interpreted as a $q$-discretization of the Beta polymer model from arXiv:1503.04117 with a new family of parameters, also attached to diagonals. To demonstrate how integrability with respect to the new diagonal parameters works, we extend the known results about Tracy-Widom large-scale fluctuations of the Beta polymer model.


Overview
During the last decades a new conjectural universality law has been actively researched under the name KPZ universality, originating from [31]. Unfortunately, despite the significant breakthroughs made in recent years we are still unable to prove or sometimes even formulate the desired universality statements. However, for some exceptional models one can directly verify the large-scale properties expected from the KPZ class models by finding explicit expressions describing the precise behavior of natural observables. Such expressions come from intricate algebraic or combinatorial properties, usually referred to as integrability, and various recent examples of large-scale analysis of integrable models can be found in [1, 3-7, 17, 25, 26, 36, 44-46] and references therein.
We focus on two integrable models, namely, the q-Hahn vertex model and the directed Beta polymer model. The former model was introduced in [38] (see also [22,24] for alternative descriptions) and it belongs to the family of stochastic solvable vertex models, which are of a particular interest for the following twofold reason. On one hand, vertex models have a rich algebraic integrability structure coming from the celebrated Yang-Baxter equation for R-matrix of the underlying quantum groups. In particular, this allows to find explicit integral expressions for q-deformed moments of naturally defined height functions, cf [7,12,16,21]. On the other hand, the stochastic vertex models can be degenerated to numerous other integrable models, including ASEP, q-Hahn and q-boson particle systems, and the mentioned directed Beta polymer model, cf. [11,16,24]. During this transition the explicit integral expressions for the vertex models reduce to expressions describing natural observables of the other models, which allow to explicitly study various properties. In this way integrability of a number of probabilistic systems can be inferred from the Yang-Baxter integrability of the vertex models.
The directed Beta polymer was originally introduced in [4] using two equivalent descriptions: it can either be viewed as a polymer model with Beta distributed weights or as a random walk in a random time-dependent Beta environment. This dual description already distinguished the Beta polymer model among the other integrable discrete polymers, namely, log-Gamma [39], strict-weak [25,36] and inverse-Beta [41] polymers, which cannot be readily interpreted as random walks in random environments (RWRE). Moreover, there is a whole cluster of integrable models obtainable as degenerations of the Beta polymer, this cluster includes strict-weak polymer, uniform sticky Brownian motions [17,28,33], Bernoulli-Exponential directed first passage percolation [4] as well as other percolation models. Going in the other direction, the Beta polymer itself can be obtained as a q → 1 limit of the q-Hahn vertex model, and this connection is closely related to the approach used originally in [4] to obtain integrability and asymptotics. 1 Overall, these numerous connections to other models give rise to a lot of interesting properties and conjectures, making Beta polymer a popular object in recent research, cf. [9,17,18,32].
The present work is devoted to a new layer of integrability of the q-Hahn vertex model, which has remained unnoticed until now. Namely, we find a new family of parameters, which are attached to the diagonals of the square lattice defining the model and which can be added to the model without breaking the integrability. To capture the most general situation containing the new parameters we introduce diagonally inhomogeneous colored q-Hahn vertex model, and the first main result of this work is the proof of integral expressions for the q-deformed moments of the colored height functions of this model. Such explicit expressions justify our claims about integrability and encapsulate the behavior of the model in a form suitable for subsequent analysis. As a direct application, we demonstrate that the shift-invariance property from [11,27] holds for our model, indicating that a broad class of symmetries is compatible with the additional parameters.
Existence of such integrable parameters attached to diagonals was surprising to us, and it is so far unique for the solvable vertex models: integrability is usually justified by the Yang-Baxter equation, which is well-suited only for the parameters attached to rows and columns of the model, but other parameters are typically incompatible with the Yang-Baxter equation. Accordingly, all other known solvable vertex models only have parameters attached to rows or columns. The q-Hahn model itself was known to have families of natural integrable column and row parameters coming from fusion of the six-vertex model, however, that construction did not indicate existence of the diagonal parameters in any way. Moreover, to the best of our knowledge the new diagonally inhomogeneous q-Hahn model cannot be obtained from any other model via known procedures. We believe that this is one of the reasons why the integrability with respect to diagonals has been hidden.
As an immediate consequence of our first result, we are able to degenerate the diagonal integrability to the Beta polymer model, finding an integrable inhomogeneous extension containing three families of parameters, attached to columns, rows and diagonals. Previously, only Beta polymer models inhomogeneous along two directions at most were known to be integrable. Moreover, to the best of our knowledge, this is the first integrable example of a directed polymer model inhomogeneous in three directions (there exist several related models with two families of parameters, studied in [5,23,42]). We expect that a lot of currently known properties of the Beta polymer can be extended to our inhomogeneous version, moreover, some new properties might emerge based on the new extension. To reinforce this claim, as the second main result of this work we find the Tracy-Widom asymptotics for the inhomogeneous model, extending one of the original results about the Beta polymer from [4].
Below we briefly describe the models and our results.

Main results: integral expressions for the diagonally inhomogeneous q-Hahn model
Let q ∈ (0, 1) and (μ 0 , μ 1 , μ 2 , . . . ), (κ 1 , κ 2 , . . . ), (λ 1 , λ 2 , . . . ) be real parameters satisfying The colored diagonally inhomogeneous q-Hahn vertex model consists of a random ensemble of colored up-right paths in a positive quadrant Z 2 ≥0 , with colors labeled by positive integers. The paths are sampled according to the following procedure: -All paths start from the left boundary of the model {0}×Z ≥1 , with b j paths of color j starting at (0, j). The numbers b j are independent and distributed according to -The remaining configuration is sampled sequentially and independently around each integer point (i, j), starting from the bottom and going along the rows. At each integer point all paths coming from the left turn upwards, while the paths coming from below might turn right following the probabilities where A = (A 1 , A 2 , . . . ) is a finite integer sequence with A c equal to the number of paths of color c entering from below, B, C, D similarly encode the colors of paths along the left, top and the right edges respectively, and we set | A| = A 1 + A 2 + · · · , | D| = D 1 + D 2 + · · · .
The name for the model comes from the probabilities above: one can notice that, for a specific choice of parameters, the probabilities resemble the orthogonality weights of the q-Hahn polynomials. On the other hand, the probabilities depend on the parameters μ i , κ j , λ j−i attached to columns i = const, rows j = const or diagonals j − i = const respectively, hence we call the model inhomogeneous, see Fig. 1 for a schematic description of the parameters. Finally, the diagonal parameters make our model unique among the other vertex models, hence inhomogeneity in this direction is emphasized.
The configurations of our model can be described in terms of colored height functions h (x,y) ≥c , which are assigned to the facets (x, y) ∈ Z + 1 2 2 formed by rows and columns and count the number of paths of color ≥ c passing below the facet, as depicted on Fig. 1. We consider the following observable: for collections x 1 ≤ x 2 ≤ · · · ≤ x k , y 1 ≥ y 2 ≥ · · · ≥ y k , x a , y a ∈ Z ≥0 + 1 2 ; c 1 ≤ c 2 ≤ · · · ≤ c k , c a ∈ Z ≥0 ;  (1.1) where the positively oriented contours Γ b are chosen to enclose μ −1 i and qΓ a for a < b, with other singularities being outside of the contours, and T τ denotes the polynomial representation of the Hecke algebra generated by the Demazure-Lusztig operators The proof is based on the idea which has already appeared in [12]: we show that both sides of (1.1) satisfy the same discrete recurrence relations, called local relations. However, the required local relations and our way of proving them are novel: we use Yang-Baxter equations for this purpose, directly connecting the existence of local relations with the integrability of vertex models. Surprisingly, one of the needed equations is the recently discovered deformed Yang-Baxter equation from [13], which has appeared in the context of spin q-Whittaker functions. Theorem 1.1 and the recurrence relations based on the Yang-Baxter equations indicate that our model is integrable, however, they do not clearly explain where the integrability and the diagonal parameters come from. To partially remedy this, let us briefly outline a more constructive approach to Theorem 1.1, which however does not work in full generality. The main idea comes from [16,21]: By definition, the average of an observable O can be written as It turns out that for certain observables O similar to Q (x,y) ≥τ.c one can relate both O(λ) and P[λ] to spin deformations of Hall-Littlewood or q-Whittaker symmetric functions, which are constructed using higher spin six-vertex model (cf. [15,16,19]) or q-Hahn model (cf. [13,20,35]). In this situation to compute E [O] one can first find a suitable Cauchy-type summation identity, which identifies the sum above with a single evaluation of a higher-spin symmetric function, and then find an integral expression for the latter.
For this constructive approach to work, one needs Cauchy-type summation identities and integral expressions for higher-spin functions, which in the context of vertex models are usually found with a "zipper"-argument based on a suitable Yang-Baxter equation. Unfortunately, we do not know how to do it for our model in full generality. However we can find suitable Yang-Baxter equations in two situations: In the previously known case when all parameters λ d attached to diagonals are equal this was done in [16,21] using Yang-Baxter equations coming from the colored six vertexmodel. In the other case, when all column parameters μ i are equal and we ignore the colors of the model, the integral expression can be constructed using the deformed Yang-Baxter equation from [13], which has an additional degree of freedom giving rise to the diagonal parameters. The exact form of our model and Theorem 1.1 were initially heuristically conjectured by taking superposition of these two cases.
and then we set The delayed Beta polymer partition function Z (r ) x,y is defined as where the sum is over all directed lattice paths π from (0, r ) to (x, y), and f (π ) = min{i | π i+1 − π i = (0, 1)}, so the path starts accumulating its weight with a delay, waiting until the first vertical step, see Fig. 2.
The partition function Z (r ) x,y also has an alternative description, in terms of a random walk in a random environment (RWRE). Fix x, y and let η i, j ∼ Beta( σ i − ρ j , ρ j − ω j−i ) be the random variables as in the polymer description above. Then, conditionally on the random environment {η i, j } i, j , we consider a random walk X t on Z, starting at X 0 = 2x − y and doing independent ±1 steps with probabilities where we use P to denote the probability space of the random walk, conditioned on the environment {η i, j } i, j . With this setup, we have Z (r ) x,y = P(X y−r ≥ −r ), which can be checked as in [4]. 2 Note that in the RWRE description the parameters σ i , ρ j , ω d are attached to the lines t − X = const, t = const and t + X = const respectively. Under a q → 1 limit the diagonally inhomogeneous q-Hahn model converges to the inhomogeneous Beta polymer, and the expression from Theorem 1.1 produces an explicit integral expression for the moments: 2 One way to check it is to identify the polymer paths contributing to Z (r ) x,y with trajectories (X t , t) of the random walk reaching the half-line X t ≥ r , X t + t = y − 2r , using the change of coordinates (i, j) = ( X +y−t 2 , y − t). x,y . Due to the delay, the dashed edges do not contribute their weights where the contour S b surrounds σ i and S a + 1 for a < b. Such integral expressions are well-suited for asymptotic analysis and the standard workflow consists of first obtaining a Fredholm determinant expression for the Laplace transform of Z (r ) x,y (which is analytic since Z (r ) x,y ∈ (0, 1) almost surely) and then performing a steep descent analysis. Below we list our results.
Consider the large scale limit of Z xt , yt when the parameters σ i , ρ j , ω d have finite number of possible values σ i , ρ j , ω d repeated with frequencies α i , β j , γ d . For example, when γ 1 = γ 2 = 1 2 we assume that half of the diagonals have parameter ω d = ω 1 while the other half has parameter ω 2 ; see Sect. 8 for a precise description of the model. We conjecture (Conjecture 8.2 in the text) that the following limit relation holds: where the slope x/y ∈ (0, 3 F GU E (s) denotes the GUE Tracy-Widom distribution [43] and the constants I , c are defined in an implicit way involving an auxiliary parameter θ ∈ (max i σ i , ∞) and polygamma functions Ψ k : 3 The range for the slope looks more natural in the RWRE description: corresponds to the average drift of the walk, and we are looking at fluctuations at non-typical velocities. See [4] for more details about the RWRE interpretation of the model.
Under assumptions dictated by purely technical reasons we are also able to partially verify our conjecture: Theorem 1.2 (Theorem 8.4 in the text) When θ ∈ (0, 1 2 ) and with the notation specified above.
In other words, we are able to partially verify the conjectural limit statement for the model having only diagonal parameters ω d < −1, with Beta distributions Beta(1, −1 − ω d ). Note that even in the homogeneous case our result is more general than the analogous one from [4], where the full proof is given only for the model with Beta(1, 1) distributions. Shortly after this paper has been written, a more general limit result for the homogeneous Beta polymer has appeared in [37], and it covers the model with Beta(α, β) distributions for arbitrary α, β ∈ R >0 . It would be interesting to see in what generality the methods and results from [37] can be transferred to the inhomogeneous model.

Layout of the paper
Sections 2-5 are devoted to the integrability of the q-Hahn vertex model with diagonal parameters: In Sect. 2 we describe the needed background on vertex models and Yang-Baxter equations, which is immediately used in Sect. 3 to establish local recurrence relations for the q-Hahn model. In Sect. 4 we define the q-Hahn vertex model and prove one of the central results of this work, Theorem 4.2, using relations from the previous section. Section 5 is a brief overview of the shift-invariance property for the vertex models, which can be established using Theorem 4.2. Then we pass to the Beta polymer side of the work, starting with Sect. 6 where we describe the transition from the q-Hahn vertex model to the Beta polymer model, and continuing with Sect. 7 where we use integral expressions coming from Theorem 4.2 to write a Fredholm determinant expression for the Laplace transform of the polymer partition function. Finally, Sect. 8 is devoted to the large-scale analysis of the inhomogeneous Beta polymer model, with the results stated in Conjecture 8.2 and Theorem 8.4.

Notation
Throughout the work we follow the following standard notation. For k ∈ Z >0 the symmetric group of rank k consists of permutations of k elements and is denoted by S k . For a permutation π ∈ S k its length is denoted by l(π ) and is defined as the number of pairs (i, j) such that i < j and π(i) > π( j). A partition λ is a finite integer sequence λ 1 ≥ λ 2 ≥ · · · ≥ λ l > 0, with l = l(λ) being called the length of λ. If λ 1 + · · · + λ l(λ) = n we write λ n. The part multiplicities m k (λ) are defined by For any n ∈ Z the q-Pochhammer symbol is defined as In particular, (x; q) 0 := 1. We also set (x; q) ∞ : . For n, m ≥ 0 the q-binomial coefficient n m q is given by With the assumption |q| < 1, all expressions above make sense both numerically and formally in the space of power series in q. For a pair of parameters (α, β) ∈ R 2 ≥0 the Beta distribution is defined by the following density function with respect to the Lebesgue measure: where Γ (z) is the Gamma function. We use Beta(α, β) to denote a Beta-distributed random variable with the corresponding parameters. For a trace-class operator K over L 2 (X ) we use det(I + K ) L 2 (X ) to denote its Fredholm determinant. The only operators K considered in this work are integral operators of the form where K (x, y) is called the kernel of K . In this case the Fredholm determinant is given explicitly by We often use Fredholm determinants over complex integration contours, so for a contour C ⊂ C we let L 2 (C) denote the space of functions on C with the complex valued integration measure dz 2π i . The Tracy-Widom GUE cumulative distribution function F GU E (r ) is defined by where Ai(z) is the Airy function and the integration contours in the last expression do not intersect and go to infinity along the prescribed directions.

Vertex models and Yang-Baxter equation
The focus of this section is the Yang-Baxter equation, which can be seen as the fundamental reason for integrability of the models presented in this work. More precisely, here we introduce our notation for vertex models, describe vertex weights and provide known forms and consequences of the Yang-Baxter equation.

Vertex configurations and weights
Vertex models are described in terms of vertices, which we graphically represent as intersections of oriented lines. Segments of the lines separated by the vertices are called edges, so each vertex has a pair of incoming edges and a pair of outgoing edges. A configuration of a vertex is an assignment of four labels to the adjacent edges. Finally, to each vertex configuration we assign a vertex weight, depending on these four labels.
Configurations of our models are labeled by compositions assigned to the edges, that is, by n-tuples of nonnegative integers A = (A 1 , A 2 , . . . , A n ) ∈ Z n ≥0 . We often present configurations as ensembles of oriented colored lattice paths directed along the edges, with colors labeled by integers 1, 2, . . . , n. 4 The two descriptions are identified by encoding the colors of paths occupying a single edge by a composition A = (A 1 , . . . , A n ), where A i is the number of paths with color i.
Throughout the text we use the following notation regarding compositions. For a composition A = (A 1 , . . . , A n ) we set | A| := A 1 + A 2 + · · · + A n . More generally, for integers 1 ≤ i, j ≤ n let Finally, given a composition A we can define an ordered | A|-tuple of integers 1 A 1 2 A 2 . . . n A n which consists of 1 repeated A 1 times, followed by A 2 copies of 2, A 3 copies of 3 and so on. The main family of vertex weights used in this work is the family of q-Hahn 5 vertex weights, denoted in one of the following ways and given explicitly by 1 − szq I [1,n] 1 − sz where i, j are integers from [1, n] satisfying i < j, I = (I 1 , . . . , I n ) is a composition and all unlisted weights are assumed to be 0 since they violate the conservation law. The parameter q is the same quantization parameter as in the q-Hahn weights, while the parameters z and s are spectral and spin parameters respectively. Both vertex weights W t,s and w z;s are particular cases of more general weights and are given explicitly by where the sum is over compositions P = (P 1 , . . . , P n ) such that P i ≤ min(B i , C i ) and we set The weights W t,s and w u;s are recovered using the following specializations: Though it is not needed for our purposes, let us briefly elaborate on the origins of these weights. The weights W N ,M z ( A, B; C, D) originate from the quantum group U q ( sl n+1 ) and can be seen as the matrix coefficients of a renormalized R-matrix acting on symmetric tensor representations, see [30]. The explicit expression (2.8) was obtained in [14] using methods of three-dimensional solvability. See [19,Appendix C] for a condensed summary of the properties of W N ,M z ( A, B; C, D) relevant for the current work. The colored higher spin six-vertex weights are coefficients of the Rmatrix above when one of the representations is the standard representation, and these weights are related to a spin deformation of the non-symmetric Hall-Littlewood polynomials, see [19] and Remark 3.8. The q-Hahn weights are not as well-understood as the other weights; we limit ourselves to stating that in one-color case they are related to spin deformations of the q-Whittaker symmetric functions, see [13,20,35].

Yang-Baxter equations
The vertex weights described above are distinguished by an algebraic relation called the Yang-Baxter equation. In the most general situation relevant to this work the Yang-Baxter equation is given in [19, Appendix C]: where both sums are taken over arbitrary triples of compositions K 1 , K 2 , K 3 . Note that both sums are actually finite due to the non-negativity of compositions and the conservation law. Setting x = y = z and using specialization (2.9) we obtain the Yang-Baxter equation for the weights W t,s : Throughout the text we usually write such equations using graphical notation. For instance, (2.11) is represented by the following diagrammatic equation: More precisely, the diagrams like the above one denote partition functions of the corresponding models. That is, given such a diagram with a fixed labelling of the boundary edges and an assignment of vertex weights to the vertices, the corresponding partition function is defined by taking the sum of products of vertex weights over all possible configurations of the internal edges.
Another instance of the Yang-Baxter equation is obtained by setting N = 1, y = z = 1 in (2.10) and applying (2.9). The result is graphically represented by Note that we again use diagrams to write the equation, but this time the lines have different thickness. Throughout the text we follow the following convention: thick edges can be labelled by any compositions, while thin edges are labeled by compositions I with |I| ≤ 1.
In [13] it was shown that the Yang Baxter equations involving the weights W t,s can be sometimes deformed by adding an additional parameter to the equation in a nontrivial way. More precisely, the following deformed Yang-Baxter equations hold:  (2.15) where η in both cases is the added deformation parameter.
Since the context of [13] required only one-colored models, the deformed Yang-Baxter equations were only proved in that case, while for the colored situation the proofs can be repeated verbatim. However, for the completeness sake, we provide the proofs for the colored setting in "Appendix A" of the arXiv version of this paper.

Stochasticity
Apart from the Yang-Baxter equation the vertex weights W t,s and w z;s satisfy another property, namely, they are stochastic: the sum of vertex weights with the same valid incoming configuration is equal to 1. In other words, the following identities hold where for the first identity A, B are arbitrary compositions, and for the second one I is an arbitrary composition and j ∈ [0, n].
The importance of the stochasticity comes from the probabilistic interpretation of the vertex models, where the stochastic weights can be used to define a stochastic sampling rule for an outgoing configuration of a vertex given an incoming one. We postpone the detailed discussion of this construction until Sect. 4.
There are several ways to establish relations (2.16), (2.17). One option is to verify everything by a direct computation, which can be readily performed for the weights w z;s , but gets more involved for the weights W t,s . Another approach is to use the general weights W N ,M z , which turn out to be stochastic as well: where A, B are arbitrary compositions satisfying the constraints (2.7). The latter fact can be derived, for example, from the fusion construction of the weights W N ,M z , see [11,Appendix].
Finally, for the weights W t,s one can use the Yang-Baxter equation to prove (2.16). This approach is described in Sect. 3, where a more general local relation is considered.

Exchange relations and Hecke algebras
Given a vertex model, it is usual to consider row-to-row transfer matrices. For a fixed L let V L denote an infinite dimensional vector space with basis enumerated by Ltuples of compositions. We denote the elements of this basis by |I 1 , . . . , I L and we let I 1 , . . . , I L | ∈ V * L denote the elements of the basis dual to {|I 1 , . . . , I L } I 1 ,...,I L . For a fixed pair of families where 1 ≤ i ≤ n and the partition function on the right-hand side consists of L vertices with the weights w uξ 1 ;s 1 , w uξ 2 ;s 2 , . . . , w uξ L ,s L . For any 1 ≤ i < j ≤ n the operators C i , C j satisfy the following commutation relations: is an operator acting on rational functions in w 1 , w 2 by Operators T (w 1 ,w 2 ) are closely related to the polynomial representation of the Hecke algebra via Demazure-Lusztig operators. The Hecke algebra H k is a C-algebra with basis {T π } π ∈S k enumerated by permutations π ∈ S k and satisfying the following defining relations where T i = T σ i denote the basis elements corresponding to the simple transpositions σ i exchanging i and i + 1. Note that the elements T i generate H k , and thus we can define a representation of H k on the rational functions in w 1 , . . . , w k by setting where s i denotes the operator exchanging the variables w i and w i+1 . 8 Note that these operators preserve polynomials in w 1 , . . . , w k , hence for any rational function f in w 1 , . . . , w k the function T π f is well defined on the domain of f , that is, T π f does not have any additional singularities compared to f .

Local relations
In this section we use the Yang-Baxter equation to derive two recurrence relations: the first relation captures the local behavior of q-moments of the colored height functions, defined in the next section, while the other relation compares rational function of certain form. While the underlying objects behind these relations are different, surprisingly, the relations themselves turn out to be identical. This observation will play a central role in the proof of one of the main results of this work, namely Theorem 4.2.

Ordered permutations
Both local relations in this section can be described in terms of order-preserving length p subsequences of a generic sequence c = (c 1 , c 2 , . . . , c r ), which can be naturally encoded by subsets But for our purposes it is more convenient to encode such subsequences using permutations. Recall that the action of the symmetric group on sequences is defined by Then any subsequence of c can be constructed as for a permutation π satisfying This leads us to the following definitions.
-ordered. With this notation the following can be readily verified:

Lemma 3.1 There is a one-to-one correspondence between integer subsets
So, each order-preserving length p subsequence of c = (c 1 , . . . , c r ) can be uniquely described as π.c [1, p] for a [1, p] × [p + 1, r ]-ordered permutation π . 9 We denote the set of such permutations by S p|r − p .

Local relation: q-moments
The first local relation describes the local behavior of a certain observable of the weights W s,t ( A, B; C, D), which are treated as a stochastic sampling rule for compositions (C, D) given ( A, B).

Proposition 3.2 For any compositions
where the sum in the right-hand side is over compositions P = (P 1 , . . . , P n ) such that P i ≤ R i .
The proof of Proposition 3.2 will be given shortly, but first we want to explain the probabilisitic interpretation of (3.1). For a sequence c = (c 1 , c 2 , . . . , c r ) and a composition [cr ,n] . and let E t,s denote the expectation with respect to this probability measure.
and arbitrary compositions A, B the following relation holds where the sum is over subsets I ⊂ {1, 2, . . . , r }, l(I) denotes the number of elements in I and Equivalently, (3.4) where the sum is over all [1, p] × [p + 1, r ]-ordered permutations and

Reduction of Corollary 3.3 to Proposition 3.2
The second part follows at once from Lemma 3.1. So we focus on the first part.
Then for any composition D we have Similarly, if c[I] = 1 P 1 2 P 2 . . . n P n for a composition P, then Q ≥c [1,i] , so to reduce Corollary 3.3 to Proposition 3.2 it is enough to show for any compositions P ≤ R that I⊂{1,...,|R|}: The latter relation readily follows from the standard expression for the q-binomial coefficients, cf. [29, Theorem 6.1]: applied for 1 ≤ k ≤ n to the decomposition I = I 1 I 2 · · · I n , with I k := I ∩ R [1,k−1] + 1, . . . , R [1,k] .

Proof of Proposition 3.2
The main idea is to see the relation (3.1) as a particular instance of the Yang-Baxter equation. We start with the Yang-Baxter equation (2.12) with the following parameters and boundary conditions where t, s, A, B, R are the same as in (3.1), while ε is an arbitrary parameter and X, Y are arbitrary compositions satisfying The next step is to take ε → 0 and to analytically continue the labels of the vertical edges. Namely, consider the following vertex weights: where α = (α 1 , . . . , α n ) ∈ C n and β = (β 1 , . . . , β n ) ∈ C n are analytically continued compositions and we use the notation αq I = (α 1 q I 1 , . . . , α n q I n ). The weights W t are obtained from the weights W t,ε by taking ε → 0 and performing analytic continuation: using (2.3) one can readily verify that . . , α n , with the highest degree term 10 having an explicit expression Multiplying both sides of (3.5) by ε −2|R| , taking ε = 0 and applying (3.6) we obtain C, D: 10 The degree of a monomial α l 1 1 . . . α l n n is l 1 + · · · + l n . for any α = (q X 1 , q X 2 , . . . , q X n ), where X = (X 1 , . . . , X n ) is a composition satisfying X + R − A − B ≥ 0. Note that we have explicitly indicated the summations over internal edges in both sides of (3.8), specifying the internal labels. In particular, both sides are clearly finite sums of polynomial functions in α 1 , . . . , α n over C(q, s, t), hence the equality for α = q X can be continued to any α = (α 1 , . . . , α n ) ∈ C n .

Remark 3.4
We have not used (2.16) during the proof of Proposition 3.2, so one can use the Yang-Baxter equation (2.12) to show that the weights W t,s are stochastic.

Remark 3.5
The local relation above also holds for more general weights W N ,M z from (2.8), with the analogue of Proposition 3.2 having the following form: (3.9) Using the same methods as in Sect. 4 below, this identity can be written in a form similar to (4.15), describing the local behavior of the height functions around a vertex with weights W N ,M z . Since (3.9) is not used in the present work we omit its proof, but the identity from Proposition 3.2 constitutes the majority of that proof: plugging the explicit expression (2.8) into the left-hand side of (3.9) one obtains a double summation which can be taken in the order allowing two subsequent applications of Proposition 3.2, readily leading to the right-hand side of (3.9) after algebraic manipulations.

Remark 3.6
For the six-vertex model the recurrence relation analogous to Corollary 3.3 can be found, in increasing generality, in [10][11][12] (in the former two the relation is called the four point relation). We note that the local relation for the q-Hahn model from the current work does not directly follows from those earlier results, and vice versa, the local relations for the six-vertex model cannot be obtained as degenerations of Corollary 3.3.

Local relation: rational functions
Now we turn to the second relation, which considers certain rational functions in w 1 , . . . , w r . Recall that the action of the Hecke algebra H r on the rational functions in w 1 , . . . , w r is generated by

Proposition 3.7
For any integer r ∈ Z ≥0 and any complex parameters t, s, λ such that s 2 = q −n for any n ∈ Z ≥0 the following equality of rational functions in w 1 , . . . , w r holds where the sum is taken over [1, p] × [p + 1, r ]-ordered permutations τ , that is, permutations satisfying Proof The main idea of the proof is to take two columns of vertices with weights w u;s , find an explicit expression for their partition function using exchange relations and then apply the deformed Yang-Baxter equation (2.15) to exchange the columns. The resulting identity turns out to be equivalent to the claim. Throughout the proof we use additional complex parameters x, y, ε and we fix a composition 1 r = (1, 1, . . . , 1) ∈ Z r . We start with defining the two-column functions mentioned above. For any composition P ≤ 1 r define and C i are the row operators (2.18). Alternatively, we can depict these functions as . . ; y . Now we want to use the exchange relations to obtain an explicit expression for f 1 r − P| P . Since x, y are fixed in this part of an argument, for now we write Ξ, S instead of Ξ x,y , S x,y . The exchange relation (2.22) implies that for any permutation π such that l(π ) < l(π σ i ), where σ i is the transposition exchanging i and i + 1. The index r − i of the inverse Hecke algebra operator comes from the fact that w are ordered in the reverse order compared to C. Iterating the identity above for a reduced expression τ = σ i 1 . . . σ i l we obtain where τ is the image of τ under the automorphism of S r sending σ i to σ r −i .
For now we focus on the partition function in the right-hand side of (3.10), namely It turns out that the partition function above can be explicitly computed for a certain choice of τ . Namely, for a composition P ≤ 1 r let τ P ∈ S r be the unique permutation satisfying In other words, τ P is the minimal permutation such that τ −1 P (r − | P| + 1), . . . , τ −1 P (r ) are the colors present in P. Then for τ = τ P the partition function (3.11) has only one configuration with non-vanishing contribution: the path of color τ −1 P (i) for 1 ≤ i ≤ r − |P| enters at row i and immediately goes upwards, exiting through the first column, while the path of color τ −1 P (i) for r − |P| < i ≤ r enters at row i, crosses the first column and exits through the second one. Since each horizontal edge can be occupied by at most one path, this is the unique possible configuration. Computing the vertex weights for this configuration gives where we set p := | P| and we use the definition of τ P to deduce Plugging the expression above for the "frozen" partition function into (3.10) we get The next step is based on the deformed Yang Baxter equation (2.15), used in the following form: Here the tilted vertices have the weights W /x, /y , while the other vertices have wweights with indicated parameters. Introducing a new operator R with coefficients we can rewrite (3.12) as Using this exchange relation, we obtain where in the last equation we have explicitly evaluated the action of R on 1 r , 0|.
Plugging the explicit expressions for the functions f 1 r − P| P gives The remainder of the proof consists of algebraic manipulations bringing the identity above to the desired form. First, we set x = t and y = 1. Then, recalling the discussion from Sect. 3.1, the summation over P ≤ 1 r is equivalent to the summation over subsets I ⊂ {1, 2, . . . , r }, which is in turn equivalent to the summation over pairs ( p, π), where p ≤ r and π ∈ S r − p| p . Moreover, tracing back these equivalences, one can readily see that P corresponds to the pair (| P|, τ P ), so we get The claim follows since the action of the Hecke algebra commutes with the multipli-

q-moments of height function
In this section we describe and prove the first main result of this work, namely, the integral expression for q-moments of the colored height functions of the colored diagonally inhomogeneous q-Hahn vertex model.

The model
As noted before, the weights W t,s ( A, B; C, D) are stochastic, so they can be used to interpret a vertex as a random sampling of an outgoing configuration given an incoming one. Now we extend this sampling construction by considering a grid of vertices with a certain choice of the spin parameters of the weights W t,s ( A, B; C, D). Then, given an incoming boundary conditions, we can sample a random configuration of the grid, that is, an assignment of a random composition to each lattice edge. Let us describe the construction in more detail.
A configuration Σ of the model is an assignment of compositions to the edges listed above. We use A (i, j) (resp. B (i, j) ) to denote the composition of the vertical (resp. horizontal) edge starting at (i, j).
Assume that we are given q ∈ (0, 1), a composition I such that |I| = N and three families of real parameters (μ 0 , μ 1 We construct a random configuration Σ of the model according to the following sampling procedure: are set to 0. -To construct compositions of the incoming left horizontal edges we first sample a collection of independent random nonnegative integers  [1,c] . In other words, the left incoming edges of the first I 1 rows have color 1, the next I 2 rows have color 2 and so on. -The compositions of the remaining edges are sampled sequentially in N 2 steps, which correspond to the vertices (i, j) of the model and are performed in lexicographical order with respect to (i, j), starting from the vertex (1, 1). -At each step we have a vertex (i, j) with already sampled compositions A (i, j−1) and due to the order of the steps. Then we treat the vertex weights W √ as probabilities for a stochastic sampling algorithm, transforming the incoming compositions into the outgoing ones A (i, j) , B (i, j) : Note that all the probabilities are nonnegative for q ∈ (0, 1) and μ i , κ j , λ d satisfying (4.1). Since the weights at vertex (i, j) depend on σ i , κ j and λ j−i , we treat these parameters as attached to column i = const, row j = const and diagonal j − i = const respectively, see Fig. 3 for a depiction of the parameters of the weights in the model.
In the description of the model we have also used that (4.2) gives a well-defined probability distribution. This follows from the identity which can be proved by setting y = q −n , n ∈ Z ≥0 and reducing the claim to the q-binomial theorem. Finally, note that the vertex sampling is trivial outside of the region i < j: since the weight W t,s ( A, B; C, D) vanishes unless A ≥ D, a vertex (i, j) with the incoming bottom configuration A (i, j−1) = 0 is forced to have a deterministic sampling with −1, j) , that is, the top outgoing configuration is equal to the left incoming one. These rules force all vertices (i, j) with i > j to have edge labels 0 around them, while the sampling at the vertices (i, i) is deterministic. In particular, the model does not depend on the parameters {λ d } d≤0 , thus we can freely omit them from now on.

Remark 4.1
The q-Hahn model described in this section originates from the work [13], where the one-colored case of this model was used to construct inhomogeneous spin q-Whittaker functions. Existence of an integrable vertex model with diagonal parameters not attached to its lines is surprising: usually integrability is caused by the Yang-Baxter equation, where parameters are attached to the lines of the model. Actually, the q-Hahn model without diagonal parameters was already known and it was previously obtained from the six-vertex model directly via fusion, see [16,20].
However, the ordinary q-Hahn model without diagonal parameters was in a certain sense "incomplete", which was indicated in [20]. Namely, the q-Hahn model can be used to construct spin q-Whittaker functions, which are dual in the sense of the dual Cauchy identity to the spin Hall-Littlewood functions, constructed using higher spin six-vertex model. The latter functions have a natural inhomogeneous version, coming from a model with three families of parameters (one family is attached to rows and two families are attached to columns of the higher spin six-vertex model), but the corresponding inhomogeneous version of the spin q-Whittaker functions turned out to be unreachable using the ordinary q-Hahn model. This was remedied with the discovery of the deformed Yang-Baxter equations in [13], which have naturally led to the vertex model construction for the inhomogeneous spin q-Whittaker functions featuring the additional diagonal parameters.

The colored height functions and the formula for q-moments
The vertex model just introduced has a collection of natural observables, called colored height functions. We denote them by h ≥c indicates the number of paths of color ≥ c passing below (x, y). More precisely, for a given fixed configuration Σ we define the corresponding height functions following two local rules: (4.5) (4.6) Due to the conservation law around each vertex, these local rules give a well-defined height function h To simplify expressions in the following sections we also introduce the following notation for the multi-point height function: for a sequence of colors c = (c 1 , . . . , c k ), a pair of sequences of half integers x = (x 1 , . . . , x k ), y = (y 1 , . . . , y k ) and a configuration Σ we set

The integral expression for q-moments
The main result of this section expresses a certain observable of the q-Hahn model as a nested contour integral. The observable in question is denoted by Q   (x,y) ≥τ.c (Σ) and for a configuration of the model Σ, an ordered sequence of colors c = (c 1 , . . . , c k ), a permutation τ ∈ S k and a pair of sequences of half integers x = (x 1 , . . . , x k ), y = (y 1 , . . . , y k ) it is defined by , where we follow the notation To formulate the result we also need to specify integration contours. They are denoted by Γ 1 , Γ 2 , . . . , Γ k and are assumed to satisfy the following conditions: and for any permutation τ ∈ S k the following holds: the following holds: Note that, contrary to the colored case, in the one-colored situation the integrand has no singularities at κ −1 j , so we need no restrictions on the positions of κ −1 j with respect to the contours Γ a .

Remark 4.4
When all parameters {λ d } d≥1 are equal and there is no inhomogeneity assigned to diagonals, the q-Hahn model can be obtained from the colored higherspin six-vertex model using stochastic fusion. In this case Theorem 4.2 can be obtained from analogous integral expressions for the q-moments of the higher-spin model: for one-color case this was done in [16], while the colored situation is covered by [12,21]. But it seems that the general case with arbitrary parameters {λ d } d≥1 cannot be reached with fusion.

Explicit integral computations
The proof of Theorem 4.2 consists of two ingredients: explicit computations of the contour integrals in certain degenerate cases, and recurrence relations given in Sect. 3. Here we focus on the former.
We fix the contours Γ a as in the theorem, and let w = (w 1 , . . . , w k ) denote the collection of the integration variables. We need the following fact describing the operators adjoint to T τ with respect to the scalar product

Proposition 4.5 ([19])
Assume that F(w 1 , . . . , w k ) = F(w) and G(w 1 , . . . , w k ) = G(w) are rational functions having only singularities of the form w a = μ −1 i , w a = κ −1 j or w a = λ −1 d . Then for any permutation τ the following relation holds: 9) or, equivalently, Proof The complete proof is given in [19,Proposition 8.1.3], see also [12, Proposition 4.1]. The main idea is to prove that the generators T i are self-adjoint, which is enough since T i generate the Hecke algebra H k . Each T i is a sum of a multiplication by a constant and w i+1 −qw i w i+1 −w i s i . The former is clearly self-adjoint, while for the latter note that the multiplication by w i+1 −qw i w i+1 −w i removes the singularity of the integrand at w i+1 = qw i , allowing to exchange the contours Γ i and Γ i+1 without changing the integral. Performing a change of variables swapping w i and w i+1 , one can see that The next statement describes the situation when one of the integrals can be easily computed, reducing the number of integration variables by one. In the context of Theorem 4.2 we will use this computation when x k = y k and the integral has no residue at w k = λ −1 d for any d, see Step 2 in Sect. 4.5.

Proof
The proof is similar to [12, Lemma 4.2.1]. Note that T τ = T ρ T τ , so by Proposition 4.5 the left-hand side can be rewritten as (4.10) Now we are taking the integral with respect to w k by computing its residues outside of Γ k . Note that τ (k) = k, hence By the conditions of the claim the function is holomorphic with respect to w k outside of Γ k , so the only non-vanishing residue of the integrand in (4.10) outside of Γ k is the residue at w k = 0. Computing it, we obtain To finish the proof, note that for any function G(w) nonsingular at w a = 0 for every a we have where σ i is a simple transposition acting on the variables w by permutation. A straightforward induction yields Another integral we need to explicitly compute by taking residues is the right-hand side of (4.7) when x 1 = · · · = x k = 1 2 . For convenience, set p i := y i − 1 2 and recall that l a corresponds to the position of the color c a on the left boundary.
and any permutation τ we have (4.11) where for every integer j ∈ [1, N ] we set that is, r j is the number of intervals [l τ −1 (a) + 1, p a ] containing j.

Proof
The proof is by induction on k and N , where at each step we either reduce k or N by 1, keeping N , k ≥ 0. The base case k = 0 is trivial. Assume that k > 0 and let l = (l 1 , . . . , l k ) and p = ( p 1 , . . . , p k ). Set We have three cases: Case 1: l k < N and p 1 < N : Then for any a we have l a , p a < N , so r N = 0, N > 0 and we can freely reduce N by 1 without changing the claim.
Case 2: l k = N : Here we can compute the integral with respect to w k by applying Proposition 4.6 for the functions Since l k = N ≥ p a for any a, after the multiplication by f k (w) the poles of g a (w) at w = λ −1 j vanish, so the function f k (w)g a (w) has no nonzero residues outside of Γ k and Proposition 4.6 implies that where the permutation τ as defined by τ = ρτ for the same cycle ρ as in Proposi- Hence, the right-hand sides of (4.11) for the triples (l, p, τ ) and (l , p , τ ) coincide.
Case 3: l k < N and N = p 1 = p 2 = · · · = p r > p r −1 : To tackle this case we start computing the integrals with respect to w 1 , . . . , w r , taking residues inside Γ 1 , . . . , Γ r . Note that possible simple poles at w 1 = κ −1 j , . . . , w r = κ −1 j coming from F l (w) are cancelled by r a=1 N j=1 (1−κ j w a ) from G p (w), so all non-vanishing residues inside Γ 1 , . . . , Γ r come from the singularities of Taking into account the geometry of the contours, the integral with respect to w 1 has only one residue inside Γ 1 , which is at w 1 = μ −1 0 . Computing this residue for the expression above, we obtain Hence, after taking the residue at w 1 = μ −1 0 , the integral with respect to w 2 again has only one non-vanishing residue inside Γ 2 , this time at w 2 = qμ −1 0 . Continuing this argument, we see that the only non-vanishing contribution in the computation of the right-hand side of (4.11) comes from the residue at w 1 = μ −1 0 , w 2 = qμ −1 0 , . . . , w r = q r −1 μ −1 0 . Note that At the same time, note that r N = r for the triple of data (l, p, τ ), while the numbers r i for i < N are the same for both triples (l, p, τ ) and (l, p , τ ). So, factoring out the term (κ N /μ 0 ;q) r (λ N /μ 0 ;q) r , we reduce the claim for (l, p, τ ) to the claim for (l, p , τ ), and the latter one allows to freely reduce N by 1, finishing the step of induction.

Proof of Theorem 4.2
The main idea of the proof is to use induction, applying the local relations from Sect. 3 to reduce the coordinates x and y. This process has to eventually reach either the case x k = y k , when we can use Proposition 4.6 to reduce k, or the case x 1 = x 2 = · · · = x k = 1 2 , in which the claim holds by Proposition 4.7. Throughout the whole argument the parameters of the model, the contours Γ a and the sequence c are fixed.
Step 1. Fix x, y, τ as in the claim and assume that 1 2 < x k < y k . In this main step we use the local relations to show that the claim for (x, y, τ ) follows from the claims for (x , y , τ ) with smaller coordinates x a ≤ x a , y b ≤ y b , y = y.
More precisely, we start by finding integers m ∈ Z ≥0 , r ∈ Z >0 such that x m < x m+1 = · · · = x m+r ≤ x m+r +1 , y m ≥ y m+1 = · · · = y m+r > y m+r +1 , where m + r ≤ k and if m = 0 or m + r = k we omit inequalities involving x 0 , y 0 or It turns out that the claim for x, y and arbitrary τ follows from the claims for x − e [m+1,m+ p] , y − e [m+1,m+r ] with p = 0, . . . , r . To show it we proceed in three substeps: first we simplify the permutation τ and then apply the local relation to both sides of (4.7).
Step 1.1. For now we want to show that it is enough to consider permutations τ which are [m + 1, m + r ]-ordered: For the left-hand side of (4.7) note that for any permutation ρ ∈ S (4.12) Note that F x,y (w) is symmetric in w m+1 , w m+2 , . . . , w m+r , hence for any ρ ∈ S [m+1,m+r ] we have Multiplying (4.13) on the left by T −1 τ min for a minimal permutation τ min of a coset in S [m+1,m+r ] \S k , one can readily see that the right-hand side of (4.7) also depends only on the left coset of τ , rather than τ itself. So, we can assume that τ is a minimal representative of a left coset in S [m+1,m+r ] \S k , or, equivalently, τ is [m + 1, m + r ]-ordered.
Step 1.2. Now, assuming that τ is [m + 1, m + r ]-ordered, we want to apply an appropriate local relation to the left-hand side of (4.7). Let v = (x m+1 − 1 2 , y m+1 − 1 2 ) denote the vertex immediately to the south-west of the facet (x, y) := (x m+1 , y m+1 ) = · · · = (x m+r , y m+r ). It has vertex weights W √ λ/κ, Now we want to change the order in which we perform the sampling of the model: first we sample the outgoing configurations for the vertices (i, j) either below v or to the left of v, obtaining a configuration Σ v , and then we sample the remaining configuration Σ. Since all our samplings are independent, their order does not matter as long as each vertex is sampled after the vertices directly below and to the left of it, so the resulting configuration is distributed in the same way as in the initial model. But note that Σ v determines the value of the height functions at the facets (x a , y a ) for a ≤ m, because they are to the left of v, and at the facets (x a , y a ) for m + r < a because they are below v. At the same time, the height functions at (x, y) are determined only by Σ v and the labels of the edges around v, which are determined by the sampling at v . Let ( A, B) denote the compositions of the incoming bottom and left edges of v, determined by Σ v , and let (C, D) denote the compositions of the outgoing top and right edges of v. Then the local rule (4.6) from the definition of the height functions implies that (4.14) and hence Q x,y ≥τ.c (Σ) = Q x,y−e [m+1,m+r ] ≥τ.c where we use the notation Q ≥c ( D) from (3.2) and we set Taking expectations we obtain

Now we can apply Corollary 3.3 to the conditional expectation inside to get
So the local relation from Corollary 3.3 gives the following relation on the left-hand side of (4.7): Step 1.3. Continuing with the same setup, we now want to apply the other local relation to the right-hand side of (4.7). Using the notation from Step 1.1, namely (4.12), we rewrite the right-hand side as Note that so we can apply Proposition 3.
Since F x,y−e [m+1,m+r ] (w) is symmetric in w m+1 , . . . , w m+r and we get Let I x,y τ denote the right-hand side of (4.7). Then, the application of the local relation above gives Step 2. Applying Step 1 we can reduce the coordinates x, y until either x 1 = · · · = x k = 1 2 or x k = y k . Here we deal with the latter case, assuming that x k = y k and showing that in this situation we can reduce k by 1 and replace the data (x, y, c, τ ) by where ρ is the cycle sending τ −1 (k), τ −1 (k) + 1, . . . , k to τ −1 (k) + 1, . . . , k, τ −1 (k) respectively.
For the left-hand side of (4.7) note that for any color c we have  x, y, c, τ ) and (x , y , c , τ ) are equal. For the right-hand side of (4.7) we are going to apply Proposition 4.6. Set Since x k = y k , the function f k (w) has no singularities outside of Γ k , as well as the functions g a (w) for any a ∈ [1, k]. Thus, the conditions of Proposition 4.6 hold and we have where ρ and τ are exactly the permutations defined above. Plugging back the expressions for f a , g a , we get the right-hand side of (4.7) for (x , y , τ , c ), as desired.
Step 3. Repeating Steps 1 and 2, at some point we either reach k = 0, when the claim is trivial, or the situation when x 1 = x 2 = · · · = x k = 1 2 . In the latter case Proposition 4.7 implies that the right-hand side of (4.7) is equal to At the same time, recall that the compositions of the incoming edges are mono-colored, and the color of the incoming edge on row j is ≥ c τ −1 (a) if and only if j > l τ −1 (a) .
Recall that b j = |B (0, j) | denote the random integers used to define the compositions at the left boundary and distributed according to ( which is equal to the value of the right-hand side, concluding the base case. Note that for the last identity we have used (4.4).

Shift invariance
In this section we show how the integral expressions from Theorem 4.2 lead to certain symmetries of the diagonally inhomogeneous q-Hahn model, which are similar to the shift-invariance symmetries of the stochastic colored six-vertex model from [11] and [27]. Our argument is identical to [12,Section 7].
To state the result we need a finer description of the height functions. For simplicity, consider a q-Hahn vertex model with incoming colors along the left boundary defined by I = (1, 1, 1, . . . ), that is, in the incoming boundary composition at row j only jth entry is nonzero. Then, for each color c we can define a cutoff point ( 1 2 , c − 1 2 ), which separates the rows with the incoming colors on the left ≥ c and < c. Since such cutoff points uniquely determine the corresponding colors, instead of using a color c and a point (x, y) to identify a colored height function h We start with the informal statement of the shift-invariance symmetry, which is rigorously described in Proposition 5.1. Given a vector of height functions h (x τ (1) ,y τ (1) ) ≥c 1 , . . . , h (x τ (k) ,y τ (k) ) ≥c k with x a , y a , c a ordered as before in Theorem 4.2, the joint k-dimensional distribution is determined solely by 11 With the exception of the artificially added 0th column. ) on the left is equal to the joint distribution of ( h (3/2,9/ ) on the right. Here the two models have the same row, column and diagonal parameters, but with different orderings corresponding to ψ Row , ψ Col , ψ Diag . To illustrate assumptions of Proposition 5.1, note that for h (3/2,9/ on the left we have rows Row The data in the first part corresponds to the one-dimensional distributions of h (x,y) ≥c , which are determined by the listed parameters in view of Theorem 4.2. Since the onedimensional distributions are determined by the joint k-dimensional distributions, one can expect this data to be necessary. On the contrary, the second part corresponds to the relative position of the height functions with respect to each other on the plane, and this data is not a priori contained in the joint distribution. Moreover, we actually expect it to be redundant, see Remark 5.2. parameters of a q-Hahn model, and  let c = (c 1 , . . . , c k ), x = (x 1 , . . . , x k ) and y = (y 1 , . . . , y k ) be sequences satisfying Additionally, let τ ∈ S k be a permutation such that y τ (a) − x τ (a) ≥ c a . Similarly, let ( c 1 , . . . , c k ), x = ( x 1 , . . . , x k ), y = ( y 1 , . . . , y k ) and τ = τ be another collection of data satisfying identical constraints, which we treat as When ρ = τ both sides of (5.1) are equal by the assumptions of the claim: for the left-hand side the factors in the product over a can be rewritten as where ρ = ρ −1 , and the statement of Proposition 5.1 essentially just tells that the expression above is equal to the analogous expression for the right-hand side of (5.1) when ρ = τ .  y ρ(a) ) ≥ c a are degenerate and equal to 0. Then both sides of (5.1) have no nonzero residues with respect to w a outside of the integration contours, and hence vanish.
-For every a we have y ρ(a) − x ρ(a) ≥ c ρ(a) and y ρ(a) − x ρ(a) ≥ c ρ(a) . Then one can verify that the conditions of Proposition 5.1 still hold when τ is replaced by ρ and, in particular, the expressions of the form (5.2) match for ρ as well.

Remark 5.2
This section is based on [12,Section 7], where the shift-invariance for the stochastic colored six-vertex model is derived directly from the integral expressions for the q-moments of the height functions. The shift-invariance for the six-vertex model was initially introduced in [11] and in full generality it was proved in [27], where flip symmetries were used as a basic block for all other symmetries. Moreover, it was shown in [27] that the relative position of the height functions, encoded by τ in our notation, is irrelevant and joint distributions of height functions are determined solely by their one-dimensional distributions. Since we don't know how the q-Hahn model with diagonal parameters can be constructed from the six-vertex model, the more general invariance result from [27] cannot be applied to our model, but we suspect it to be true.

Reduction to Beta polymer
In this section we follow [11] in order to reduce the results from Sect. 4 to the Beta polymer model. See also [12,21] for similar arguments.

Continuous model
As before, we treat q-Hahn vertex model from Sect. 4 as a sequence of independent stochastic samplings corresponding to the vertices, depending on the parameters μ i , κ j , λ d . We slightly simplify the model by assuming that the incoming colors from the left are different and are encoded by the composition I = (1, 1, 1, . . . ). We consider the following limit regime: where σ i , ρ j , ω d are new real parameters satisfying It turns out that the limiting model in this regime is described in terms of the Beta distribution. Below we present without proofs the exact description, referring to [11] for the general case, and to [4] for the one-colored case.
to the following sampling procedure for the outgoing masses γ = (γ 1 , . . . , γ n ) and δ = (δ 1 , . . . , δ n ) given the incoming masses α = (α 1 , . . . , α n ) and β = (β 1 , . . . , β n ): where we set δ ≥k = l≥k δ l and α ≥k = l≥k α l . The remaining masses γ k are defined by The two statements above result in the following model, serving as the limit under (6.1) of the inhomogeneous q-Hahn vertex model and depending on the parameters σ i , ρ j , ω d : -Each lattice edge of the positive quadrant Z 2 ≥1 has a random real-valued vector u = (u 1 , u 2 , u 3 , . . . ) attached to it. Coordinate u i ≥ 0 is interpreted as the mass of color i; -To sample a configuration of the model we first sample independent Betadistributed variables η i, j ∼ Beta(σ i − ρ j , ρ j − ω j−i ) for all i ≥ 0, j ≥ 1; -For the edges entering the quadrant from below all masses are set to 0; -For the edge entering the quadrant at row j only the mass u j is non-zero and it is equal to − ln(η 0, j ); -Given the incoming masses of colors α, β entering into a vertex (i, j) from the bottom and from the left, the outgoing masses γ , δ exiting the vertex to the top and to the right are defined through where δ ≥k = l≥k δ l and α ≥k = l≥k α l . -This definition implies that only colors ≤ j pass through vertex (i, j), guaranteeing that at any vertex we are working only with finite sequences of masses.
We can define the colored height functions h (x,y) ≥c of the continuous model in the same manner as for the q-Hahn model: we set h (x, 1 2 ) ≥c = 0 and The following generalization of [11, Corollary 6.22] straightforwardly follows from the statements above and clarifies in what sense the continuous model is the limit of the q-Hahn model.
The model just described is an inhomogeneous generalization of the Beta polymer model from [4], with three families of the parameters. See Sect. 1 for a graphical description of Z (r ) x,y . It turns out that the continuous vertex model described earlier can be identified with the delayed inhomogeneous Beta polymer model, in the way identical to the homogeneous situation [11]: repeating the proof of [11, Proposition 7.1] verbatim one can verify that for (x, y) ∈ Z 2 ≥0 the height functions h

Proposition 6.5
For any sequences (r 1 , . . . , r k ), (x 1 , . . . , x k ) and (y 1 , . . . , y k ) satisfying and for any permutation τ ∈ S k the following holds: ≥r τ −1 (a) +1 and assuming I = (1, 1, . . . ). Set For the left-hand side of (4.7) we can apply Corollaries 6.3 and 6.4 to obtain x,y as ε → 0. For the right-hand side of (4.7) we perform a change of coordinates w a = 1 − εv a . One can readily see that for fixed contours S 1 , . . . , S k satisfying the conditions of the claim and a sufficiently small ε the contours 1 − εS a satisfy the conditions for the contours Γ a , and so in the right-hand side after a change of contours we have an integral over fixed compact contours, with integrand depending on ε. Point-wise we have , One can readily check that the integrand is uniformly bounded for fixed contours S a and sufficiently small ε, so by dominated convergence the claim follows.
As a particular case of Proposition 6.5, we have an integral expression for the moments of a single partition function Z (r ) x,y of the inhomogeneous Beta polymer. Note that, up to a shift of parameters, one dimensional distributions of Z (r ) x,y depend only on x and y − r (this is the vector between end points of paths inside the polymer), so without loss of generality we can set r = 0 and consider the moments of Z x,y := Z (0) x,y : Corollary 6.6 For any (x, y) ∈ Z 2 ≥0 such that y ≥ x and any k ∈ Z ≥0 we have

4)
where the integration contours S a encircle σ i , leave ω d outside and S b encircles both S a and S a + 1 for a < b.

Laplace transform for inhomogeneous Beta polymer partition function
From now on we focus on the partition function Z x,y = Z x,y of the inhomogeneous Beta polymer, with parameters σ i , ρ j , ω d . Our goal is Theorem 8.4, which partially establishes the Tracy-Widom large-scale behavior and which is proved in next section. Our result generalizes the corresponding result from [4], moreover, our proof follows the same outline: in this section we find a suitable expression for the Laplace transform of Z x,y in terms of a Fredholm determinant, which is asymptotically analyzed in next section to establish the limit theorem.
Our starting point is the integral expression (6.4) for the moments of Z x,y . Fix x, y and the parameters of the model σ i , ρ j , ω d , and assume that Since our aim is an integral expression for the Laplace transform of Z x,y which involves all the moments, it is convenient to replace the integration over shifted contours S 1 , . . . , S k in (6.4) by an integration over a single contour C. This is done with the help of the following statement: Then where the sum in the right-hand side is over partitions λ of k, l(λ) denotes the number of nonzero parts of λ and m i denotes the number of parts of λ equal to i.

Rough idea of the proof
The statement above follows a type of deduction called contour shift argument. The main idea is to shrink the contours S 2 , . . . , S k in the left-hand side to the smallest contour S 1 , picking up the residues at v b = v a + 1 for b > a.
Due to the conditions on the function f and the contours, all non-zero residues along this deformation come from the poles of a<b , and a careful bookkeeping produces the right-hand side.
A detailed explanation of the full argument can be found in [2, Proposition 3.2.1] and [8,Proposition 7.4], where more general q-deformed versions of the argument are presented. To adopt those claims to our setting one have to repeat the argument similar to the proof of Proposition 6.5, see [2,Proposition 6.2.7] and [4,Proposition 3.6]. Let and choose C = S 1 to be a sufficiently small contour around the points σ i such that the shifted contour C + 1 and the points ω d are completely outside of C. Such contour exists due to the conditions (6.2), (7.1). Corollary 7.2 With the notation above, we have Proof Direct application of (7.2).
Now we want to take a sum over k in (7.3) to obtain the Laplace transform. But to resolve possible convergence issues we need to estimate the integrand for large λ. Set where Γ (z) is the gamma function. Then

Proposition 7.4 For any u ∈ C we have
where det(I + K ) denotes the Fredholm determinant of K , and K u is the integral operator with the kernel Proof The proof essentially repeats [2, Proposition 3.2.8]. First we rewrite (7.3) as and then change the summation over partitions to a summation over arbitrary sequences, absorbing the multinomial coefficient: (7.6) Now we want to take the sum over k, but we need to be careful about possible convergence issues. For the left-hand side we have |Z x,y | ≤ 1 almost surely, so the sum converges absolutely for any u and by dominated convergence gives the Laplace transform: For the right-hand side we can use the conditions on the parameters σ i , ρ j , ω d and on the compact contour C to give a uniform in At the same time, by Lemma 7.3 we have d=1 ω d and the constants of O n −1 can be chosen independently of v ∈ C, since C is bounded. Hence for some constants Λ 2 , ε > 0 and any n ∈ Z ≥1 we have Gathering all pieces, we obtain | K u (v, n; v , n )| < Λ 1 Λ 2 n −εn |u| n and hence, by Hadamard's inequality, where ξ > 0 is the length of the contour C. Since for any u ∈ C, the right-hand side of (7.6) absolutely converges.
So, we have obtained an expression for the Laplace transform in terms of a Fredholm determinant over L 2 (Z ≥1 × C). But the summation over Z ≥1 is inconvenient for the asymptotical analysis, so we replace it by another integration using Mellin-Barnes integral formula. To proceed we add another assumption, namely, we assume that there exists a vertical line L = {z | Re[z] = h} for some h such that the contours C and C + 1 are separated by L. Note that this is still possible because of (7.1).

Proposition 7.5 For any u
we use the principal branch of the natural logarithm.
Proof The statement and the proof are similar to [4,Lemma 3.7], but our restrictions on u are different.
Since the contour C is compact, there exists N 0 ∈ Z >0 such that the contour C is entirely inside the rectangle formed by the lines L, has nonzero residues only at the poles of π sin(π(z−v)) , since all the residues of g (v) g(z) are to the left of the line L. Computing the residues at z = v + n for n = 1, . . . , N we get Now we want to take N → ∞ in both sides. For the left-hand side, using the same upper-bound on K (v, n; v n ) as in Proposition 7.4, we have so the absolute convergence holds for any u ∈ C. For the right-hand side we need to show that the integral along the symmetric difference of L and H N converges to 0. This symmetric difference D N consists of the line segments sequentially connecting the points Clearly, 1 z−v for z ∈ D N can be bounded uniformly in N , while π sin(π(z−v)) decays exponentially in Im[z]: where Λ 3 does not depend on N . To bound 1 g(z) we use Lemma 7.3 to get which is enough to conclude that the integral over D N tends to 0.

Proposition 7.6
For any u ∈ C such that arg(−u) ∈ (− π 2 ; π 2 ) we have where the integral operator K u is defined by the kernel Here C is a sufficiently small contour containing σ i for all i, and L is a vertical line separating C and C + 1.
Proof From Proposition 7.4 we have Due to the dominated convergence, which is valid since the contour C is compact and we have an integrable converging upper-bound (7.7), we can exchange the second summation and the integration. Since K u (v, n; v , n ) does not depend on n , we can use use the same absolute upper-bound (7.7) and multi-linearity of the determinant with respect to rows to get .
Then the claim follows from Proposition 7.5.

Limit theorem
This section is devoted to the partial result about the Tracy-Widom large-scale limit of the partition function Z x,y of the inhomogeneous Beta polymer model, when x, y → ∞ at a constant ratio and the parameters of the model have finite number of values repeated with certain frequencies. Our approach follows [4]: we perform a steepest descent analysis of the Fredholm determinant from Proposition 7.6.

Overview
Fix k and consider three families of parameters {σ 1 , . . . , σ k }, {ρ 1 , . . . , ρ k } and {ω 1 , . . . , ω k }. Let {α 1 , . . . , α k }, {β 1 , . . . , β k }, {γ 1 , . . . , γ k } be collections of corresponding "frequencies" satisfying For a fixed slope x/y we consider the partition function Z xt , yt as t → ∞, where parameters of the model { σ 0 , . . . , σ xt }, { ρ 1 , . . . , ρ yt }, { ω 1 , . . . , ω yt − xt } 13 are defined in the following way: For the column parameters, exactly α [t] i xt out of xt +1 parameters σ i are equal to σ i , where α [t] i is a sequence satisfying Similarly, the row parameters ρ j and the diagonal parameters ω d are repeated for β [t] j yt and γ [t] d (y − x)t times respectively, where Motivated by the homogeneous setting [4] we are looking for a result of the form where F GU E ( p) is the GUE Tracy-Widom distribution [43]. A probabilistic interpretation, which one can find more natural, is given by identification of Z xt , yt with the distribution of a random walk in a random Beta environment. From this point of view, we are looking at large deviations of the random walk, studying fluctuations of a random rate function (depending on the random environment) from the expected limit I (x/y) depending on the slope. See [4] for more details.
To prove a limit relation ( given that the limit is a continuous distribution function. So, the problem is reduced to studying the Laplace transform, which by Proposition 7.6 has an expression in terms of a Fredholm determinant. The latter can be analyzed using the steepest descent method: As t → ∞ the function h [t] (z) converges to the function which, for the right choice of I , has a second order critical point θ depending on the slope y/x in a nontrivial way. However, the dependence of the slope x/y on θ can be easily described, suggesting the following parametrization by θ ∈ (max i σ i , ∞): where Ψ (z) denotes the digamma function and Ψ k (z) denotes the kth polygamma function: (2) c(θ ) can be chosen to be positive.
Proposition 8.1 is proved in "Appendix B" of the arXiv version of this paper. 4 ) for z close to θ , and the general philosophy of the steepest descent method suggests that asymptotical behavior of the Fredholm determinant is dominated by its behavior in a small neighborhood of θ , which, after certain manipulation, can be shown to produce F GU E (r ). Unfortunately, the general case turns out to be technically challenging, and we do not pursue it here. In this section we prove Conjecture 8.2 for parameters σ i , ρ i , ω i , θ satisfying the following restrictions: Assumption 8. 3 We assume that -All σ i are equal to 0, while all ρ i are equal to −1:

Conjecture 8.2 For any
θ ∈ 0, 1 2 . The restrictions from Assumption 8.3 come from two technical difficulties. Namely, to justify the steepest descent argument one needs to -Find contours C θ , L θ passing through θ and satisfying Re[h(z) − h(v)] ≤ 0 for z ∈ L θ , v ∈ C θ , with equality only if v = z = θ ; -Prove that one can deform the initial contours C, L to C θ , L θ without changing the limit of the Fredholm determinant.
For the first step a suitable descent contour L θ for Re[h(z)] valid for any parameters σ i , ρ i , ω i was found in [4] and it is given by a vertical line passing through θ . At the same time, finding a descent contour for − Re[h(v)] is much more involved for general parameters: we know how to find it only if σ i −ρ j = 1 for any i, j and θ −max i σ i < 1, forcing the first part of Assumption 8.3. 15 Under these restrictions another contour from [4] turns out to be sufficient, we denote it by C θ and it is equal to a circle around 0 of radius θ . However, numerical computations and related results suggest that a descent contour for − Re[h(v)] should exist in general: see [17] for an example of a significantly more involved argument used to lift the analogous assumption in a different but related setting. The second step is complicated by potential singularities of the Fredholm determinant along the deformation of C, L to C θ , L θ . One can hope to estimate the residues at these singularities and show that their contribution asymptotically vanishes, but the required argument is technically challenging. The second part of Assumption 8.3 allows us to avoid this issue, since for θ < 1 2 the contours can be deformed without crossing any singularities of the relevant integrands.
The remainder of Sect. 8 is devoted to the proof of Theorem 8.4. As clarified above, it is enough to prove where u [t] = − exp(t I − t 1/3 cr ). The proof is structured as follows: first we describe the needed properties of the function h(z), as well as specify estimates on convergence of h (t) (z) to h(z). Then we prove various bounds on the kernel K u from Proposition 7.6. Finally, to prove Theorem 8.4 we reduce the computation of the limit to a small neighborhood of θ , perform a change of variables to get rid of lower degree terms and use standard algebraic manipulations to establish (8.3).

h-functions and descent contours
Here we list the needed properties of the functions h(z) and h [t] (z), which are defined by [t] θ L θ Fig. 6 Contours C θ and L θ , as well as the t-dependent deformation L [t] θ needed to avoid singularity at with x, y, I defined in terms of θ as before.
We start with the properties of h(z). Define the contour L θ as a vertical line passing through θ , oriented upwards: Under Assumption 8.3 we also define the contour C θ = {θ e iφ | φ ∈ [−π, π]} oriented counter-clockwise. See Fig. 6 for a depiction of the contours.
The proofs of the following three statements are technical and are omitted from this text, they can be found in "Appendix B" of the arXiv version of this paper. The first statement does not require Assumption 8.3.  Let us briefly sketch the main idea of the omitted proofs of Propositions 8.6-8.8. It turns out that the function h(z) is a convex combination of the analogous functions for the homogeneous model, namely This allows one to reduce the problem to the homogeneous case and analyze the functions h(z) in a way similar to [4,Section 5]. See "Appendix B" of the arXiv version for more details. The next statement reiterates the fact that θ is a critical point of h(z).

Lemma 8.9
For z in a neighborhood of θ we have where c is given by (8.2). Moreover, c > 0.
Proof Readily follows by considering the Taylor expansion of h(z), taking derivatives and using the definitions of x(θ ), y(θ ), I (θ ) and c(θ ). The positivity of c follows from Proposition 8.1.
Finally, we need the following statements regarding the convergence of h [t] (z) to h(z) in various situations.

Lemma 8.10
Assume that a sequence of points z t ∈ C converges to θ in a way that h [t] (z t ) is well defined for all t > 0. Then lim t→∞ where Λ 1 is the same constant as in the proof of Lemma 8.11. Applying the asymptotic expansion for the Gamma function to the right-hand side implies the claim.

Kernel estimates
Recall that the substitution of u = − exp(t I − t 1/3 cr ) into the integral kernel from Proposition 7.6 leads to the following kernel: For v, v = θ we interpret the integral above using its principal value as the contour L θ passes to the right of θ . In other words, we slightly deform the contour L θ to an upwards oriented contour L [t] θ consisting of vertical lines {θ ± Y i | Y ∈ [t −1/3 , ∞)} and a semicircle {θ + t −1/3 e φi | φ ∈ [−π, π]}, see Fig. 6.
The aim of this section is to give asymptotical bounds on behavior of the kernel. We use the following notation. For ε > 0 let L [t] θ,ε and C θ,ε denote the parts of the contours L [t] θ and C θ contained inside the ε-neighborhood of θ . Define the kernel Finally, define K(ε) = D ε ∪ A ε ∪ (C θ \C θ,ε ) as a contour consisting of the upwardoriented wedge D ε = {θ + Re ±2π i/3 | R ∈ [0, ε]}, a part of the contour C θ outside of the ε-neighborhood of θ , which we denote as C θ,≥ε , and a couple of upward-oriented arcs where ψ(ε) is chosen so that the result is connected. See Fig. 7 for a depiction of the contour K(ε). For the following statements we fix ε > 0.

Lemma 8.13
Let v, v ∈ K(ε) be a pair of points such that Re[h(v) − h(θ )] ≥ b for some fixed b ≥ 0. Then for t > ε −3 we have where the constants Λ, a > 0 depend only on ε and the parameters of the model.
Then the whole integral is bounded as follows: for some constant Λ 3 depending only on ε.
The claim follows by combining the bounds for P 1 and P 2 , . Recall that, in view of Proposition 7.6, it is enough to prove lim t→∞ det(I + K [t] ) L 2 (C θ ) = F GU E (r ), where we use the notation from the previous subsection.
Step 1. First we show that the limit is dominated by the behavior in a neighborhood of θ . Namely, in this step we show that for sufficiently small fixed ε > 0 we have It is more convenient to proceed in two steps, reducing first the contour C θ to D ε and then the kernel K [t] to K [t] ε . Note that we can deform C θ to K(ε) = D ε ∪A ε ∪(C θ \C θ,ε ) without changing the value of the determinant. Since C θ is a steep descent contour for − Re[h(z)], for each ε there exists a constant b 1 > 0 such that and using multi-linearity of determinants and Lemmas 8.13 and 8.14 we have the following upper-bound for some constants Λ 3 , a > 0 Step 2. After getting rid of the global parts of the integration contours, in this step we take the needed limit by scaling the neighborhood of θ by t 1/3 .
For now, fix ε > 0 which is sufficiently small so that the argument in the previous step works, but otherwise arbitrary. We would like to perform the following change of coordinates in the Fredholm determinant det(I + K [t] ε ) L 2 (D ε ) and the kernel K [t] ε : To take the limit as t → ∞ we start with the point-wise convergence of the kernel ε be the contour consisting of vertical lines {±Y i | Y ∈ [1, εt 1/3 ]} and a semicircle {e φi | φ ∈ [−π, π]}, oriented so that the imaginary part increases. Also, let L denote the limiting contour with the vertical lines {±Y i | Y ∈ [1, ∞]} and the same semicircle. Then Clearly −π t −1/3 for some Λ 3 > 0, which combined with Hadamard's bound is sufficient for an integrable and absolutely summable upper-bound on the kernel of the Fredholm determinant. Hence with the kernel K ∞ defined in (8.7).
Step 3. To finish the proof we need to show that the Fredholm determinant in the right-hand side is indeed equal to F GU E , which can be done by standard manipulations.
Rescaling v, z by c we can assume that Since the integrand has quadratic decay, we can deform the contour L to a contour going from ∞e −iπ/3 to ∞e iπ/3 with Re[ z] > 0. Now we can apply the same argument as in [6,Lemma 8.6]: For any z such that Re[ z] > 0 we use