Correlation Functions of XX0 Heisenberg Chain, q-Binomial Determinants, and Random Walks

The XX0 Heisenberg model on a cyclic chain is considered. The representation of the Bethe wave functions via the Schur functions allows to apply the well-developed theory of the symmetric functions to the calculation of the thermal correlation functions. The determinantal expressions of the form-factors and of the thermal correlation functions are obtained. The q-binomial determinants enable the connection of the form-factors with the generating functions both of boxed plane partitions and of self-avoiding lattice paths. The asymptotical behavior of the thermal correlation functions is studied in the limit of low temperature provided that the characteristic parameters of the system are large enough.


Introduction
The exactly solvable Heisenberg XXZ model is a prominent model describing the interaction of spins 1 2 on a chain. The integrability of the model via the algebraic Bethe Ansatz has led to important results, going from the spin dynamics up to the exact expressions for the correlation functions [1][2][3][4][5][6][7][8][9][10][11].
The XX0 Heisenberg chain is the zero anisotropy limit of the XXZ model, it also may be considered as a special free fermion case of the XY magnet [12,13]. It appears that XX0 model is related to many mathematical problems. It is related to the theory of the symmetric functions [14] and to the theory of plane partitions. Plane partitions (three-dimensional Young diagrams) [14][15][16] were then discovered to be connected with amazingly wide ranging problems in mathematics as well as theoretical physics. They are intensively studied, e.g., in probability theory [17,18], enumerative combinatorics [19], theory of faceted crystals [20,21], directed percolation [22], topological string theory [23], and the theory of random walks on lattices [16,[24][25][26].
The correlation functions of the XX0 chain are of considerable interest, and their behavior was intensively investigated for the system in the thermodynamic limit [9,[27][28][29]. In our paper we study the asymptotical behavior of the thermal correlation functions in the limit of low temperature provided that the chain is long enough while the number of flipped spins is moderate. Namely in this limit the thermal correlation functions are related to random matrix models [25]. This connection allows to uncover, in particular, the mapping between the correlation functions and the low energy sector of quantum chromodynamics [29].
We shall consider the XX0 Heisenberg model on the periodical chain. The representation of the Bethe wave functions via the Schur functions [14] allows to apply the well-developed theory of the symmetric functions to the calculation of the thermal correlation functions as well as of the form-factors. In the present paper we are interested in the correlation functions of two types: the correlation function of the states with no excitations on n consecutive sites of the chain that will be called persistence of ferromagnetic string, and the correlation function of the creation operator of the n excitations on the consecutive sites of the chain that will be called persistence of domain wall. Special attention will be paid to the combinatorial objects appearing in the calculations (the generating functions of plane partitions and random walks, the q-binomial determinants) and to the combinatorial interpretation of the obtained results. We will calculate the leading terms of their asymptotics, provided that the characteristic parameters of the system are large enough, including the critical exponents of these correlation functions in the low temperature limit, and the related amplitudes. These amplitudes are found to be proportional to the squared numbers of boxed plane partitions.
The paper is structured as follows. Section 1 is introductory. The XX0 model and its solution are presented shortly in Section 2, the considered correlation functions are defined and the amplitudes of the state vectors are written in terms of Schur functions. This representation allows to calculate the form-factors of operators in Section 3 applying the formulas of the Binet-Cauchy type. The persistence of ferromagnetic string as well as the persistence of domain wall are also calculated in this section. In Section 4 we deal with the combinatorial aspects of the problem. The q-binomial determinants are introduced and their connection with the generating functions of plane partitions is discussed. It is shown also that the form-factors, obtained in the previous section, under the special parametrization are expressed as the generating functions of boxed plane partitions and of the self-avoiding lattice paths. The asymptotical estimates of the correlation functions are obtained in Section 5. Discussion in Section 6 concludes the paper. In Appendix I we provide some notions concerning boxed plane partitions and their generating functions. The proof of the determinantal formulas crucial for this paper is given in Appendix II.

XX0 Heisenberg model and outline of the problem
The Heisenberg XX0 model on the chain of M + 1 sites is defined by the Hamiltonian Here the periodic boundary conditions σ # k+(M +1) = σ # k are assumed. The local spin operators σ ± k = 1 2 (σ x k ± iσ y k ) and σ z k obey the commutation rules: . The spin operators act in the space H M +1 spanned over the states M k=0 |s k , where |s k implies either spin "up", |↑ , or spin "down", |↓ , state at k th site. The states |↑ ≡ 1 0 and |↓ ≡ 0 1 provide a natural basis of the linear space C 2 . The sites with spin "down" states are labeled by the coordinates µ i , 1 ≤ i ≤ N. These coordinates constitute a strictly decreasing partition µ = (µ 1 , µ 2 , . . . , µ N ), where M ≥ µ 1 > µ 2 > . . . > µ N ≥ 0. The other important partition is λ = (λ 1 , λ 2 , . . . , λ N ) of weakly decreasing non-negative integers: L ≥ λ 1 ≥ λ 2 ≥ . . . ≥ λ N ≥ 0. The elements λ j are called the parts of λ. The length of partition l(λ) is equal to the number of its parts. The sum of all parts is the weight of partition, |λ| = N i=1 λ i . Partitions λ can be represented by Young diagrams. The Young diagram of λ consists of N rows of boxes aligned on the left, such that the i th row is right on the (i + 1) st row. The length of the i th row is λ i . The relation λ j = µ j − N + j, where 1 ≤ j ≤ N, connects the parts of λ to those of µ. Therefore, we can write: λ = µ − δ N , where δ N is the strict partition (N − 1, N − 2, . . . , 1, 0). There is a natural correspondence between the coordinates of the spin "down" states µ and the partition λ expressed by the Young diagram (see Fig. 1). Throughout the paper bold-faced letters are used as short-hand notations for appropriate N-tuples of numbers.
The N-particle state-vectors | Ψ(u N ) , the states with N spins "down", is convenient to express by means of the Schur functions [30,31]: where summation is over all partitions λ satisfying M ≡ M + 1 − N ≥ λ 1 ≥ λ 2 ≥ · · · ≥ λ N ≥ 0. The parameters u N = (u 1 , u 2 , . . . , u N ) are arbitrary complex numbers, and u 2 N ≡ (u 2 1 , u 2 2 , . . . , u 2 N ). The state |⇑ in (2) is the fully polarized one with all spins "up": |⇑ ≡ M n=0 |↑ n . The amplitudes in (2) are expressed in terms of the Schur functions S λ [14]: The conjugated state-vectors are given by If parameters u 2 j ≡ e iθ j (1 ≤ j ≤ N) satisfy the Bethe equations [27], then the state-vectors (2) become the eigen vectors of the Hamiltonian (1): Here (and throughout the paper) the notation θ N for N-tuple (θ 1 , θ 2 , . . . , θ N ) is reserved for the solutions to the Bethe equations (6), and e iθ N ≡ (e θ 1 , e θ 2 , . . . , e θ N ). The solutions θ j to the Bethe equations (6) can be parametrized such that where I j are integers or half-integers depending on whether N is odd or even. It is sufficient to consider a set of N different numbers I j satisfying the condition: M ≥ I 1 > I 2 > · · · > I N ≥ 0.
Then the eigen energies in (7) are equal to The ground state of the model is the eigen-state that corresponds to the lowest eigen energy E N (θ g N ). It is determined by the solution to the Bethe equations (8) at I j = N − j: and is equal to In the present paper, the two types of the the thermal correlation functions in a system of finite size will be considered. We call them the persistence of ferromagnetic string and the persistence of domain wall. The persistence of ferromagnetic string is related to the projection operatorΠ n that forbids spin "down" states on the first n sites of the chain: where H and θ g N are given by (1) and (10), respectively, and β ∈ C. Some results on this correlation function have been reported in [31].
The persistence of domain wall is related to the operatorF n that creates a sequence of spin "down" states on the first n sites of the chain: whereF + n is the Hermitian conjugated operator acting on the conjugated state-vectors (5), and θ g N −n is the set of ground state solutions to the Bethe equations (6) for the system of N − n particles: We assume thatΠ 0 andF 0 are the identity operators so that T (θ g N , 0, β) = 1 and The calculation of introduced correlation functions will be based on the Binet-Cauchy formula [32] adapted for the Schur functions: where the summation is over all partitions λ satisfying: L ≥ λ 1 ≥ λ 2 ≥ · · · ≥ λ N ≥ n. The Vandermonde determinant (4) is used in (14), and the entries T kj are given by:

The correlation functions
In this section we shall calculate the persistence correlation functions of ferromagnetic string and of domain wall. The calculation is natural to start with the derivation of the form-factors of appropriate operators.

The Bethe states and form-factors
Applying the orthogonality relation we reduce the calculation of the scalar products of the state-vectors (2) and (5) to the Bine-Cauchy formula (14): where summation is over all partitions λ with at most N parts, each of which is less or equal to M. The entries T o kj in (17) are given by (15) taken at n = 0: For v N = u N , this scalar product is equal to the squared "length" of the states (2): On the solutions u 2 N = v 2 N = e iθ N (8) to the Bethe equations (6), the entries (18) are equal to and for the square of the norm N 2 (θ N ) ≡ N 2 (e iθ N /2 ) we have: Notice, that if θ N and θ ′ N are different sets of the solutions to the Bethe equations, the related eigen vectors are orthogonal [11]: Ψ(e iθ N /2 ) | Ψ(e iθ ′ N /2 ) = 0. Being a complete and orthogonal set of states, these eigen vectors provide the resolution of the identity operator [3,4] I = where summation is over all different solutions to the Bethe equations (8).
The form-factor of the projectorΠ n (11) is defined by the ratio: Taking into account (2), (14) and (15), we find that the numerator of (22) is equal to: where the entries are equal to Taken at n = 0, Eq. (23) reproduces the answer for the scalar product (17).
On the solutions to the Bethe equations (8), the form-factor (23) takes the form with the entries of matrix K n (θ k , θ j ) equal to To obtain this answer we have subtracted and added the unity to the numerator of T kj (24) and used the equality (19). Finally, for the form-factor (22) we obtain In this form it is known as the emptiness formation probability, being the probability of formation of a string of n consecutive spin "up" states [6,8,10,11,27]. Let us consider then the form-factor of the domain wall creation operatorF n (12): To calculate this transition element, we first introduce an auxiliary operator D n (u) acting on an expectation · u considered as function of u: where The definition (28) implies that the expectation · u is first multiplied by the ratio of the Vandermonde determinants and then differentiated n times. Now we are ready to formulate the following Proposition 1 The action of operator D n (u) expressed by (28), (29), (30) on the scalar product Ψ(v N ) | Ψ(u N ) gives the form-factor of the domain wall creation operator F n (27): Proof From the definition of the state-vectors (2) and of the operatorF n (12) we obtain the representation of the form-factor: where summation runs over the partitions λ of the length N − n: The corresponding strict partitions are given bŷ µ =λ + δ N and µ = λ + δ N −n .
To derive (32), let us act byF n on the state | Ψ(u N −n ) given by (2) with the summation taken over λ ⊆ {(M + n) N −n }. The operatorF n acts non-trivially only on those vectors in |Ψ(u N −n ) , that do not contain spin "down" states on the first n sites: We have used the definition of the Schur function (3) to obtain the last equality. Applying the orthogonality relation (16) to right-hand side of (33), we find that right-hand side of (32) indeed holds true. Eventually, direct evaluation of right-hand side of (31) leads to right-hand side of (32) provided that the scalar product is represented through the Schur functions according to (17).
Proposition 1 enables us to obtain two summation rules for the products of the Schur functions, which are crucial in establishing of the combinatorial results for the correlation functions in question. So, we formulate the following Proposition 2 The following sums of products of the Schur functions take place: where the entries of the matrices (T kj ) 1≤k,j≤N and (T kj ) 1≤k,j≤N are: Proof Let us calculate right-hand side of (31) using the determinantal form of the scalar product (17): Taking into account the explicit form of the entries T o kj (18) we obtain: where the matrixT is given by (36). Since right-hand sides of (32) and (38) mutually coincide, the relation on the Schur functions (34) (which is of the type of (14)) holds true.
Repeating the arguments of Proposition 1, we find that the form-factor of the conjugated operatorF + n is equal to and to

Persistence of ferromagnetic string
Let us recall the main relations concerning the persistence of ferromagnetic string and its relationship with the problem of vicious walkers in the random turns model [33]. The problem of enumeration of the vicious walkers is actively investigated [34][35][36][37][38][39][40][41]. The random walks across one-dimensional periodic lattice are closely related to the correlation functions of the XX0 magnet [25,26,30].
Taking into account the explicit form of the state vectors (2) and (5), we obtain the following answer [31] for the matrix element parametrized by arbitrary u N and v N . The range of two independent summations over λ L and λ R is defined in (14), and µ L,R = λ L,R + δ N denote the corresponding strict partitions. The transition amplitude: is related to enumeration of paths of N vicious walkers moving across the sites of a onedimensional chain [25,26,30,31]. The expression (41) at β = 0 is in agreement with (23). The transition amplitude (42) satisfies the difference-differential equation derived in [26]. The corresponding solution has the determinantal form: where the entries are the transition amplitudes F k; l (β) = ⇑ |σ + k e −βH σ − l | ⇑ , which are (42) for N = 1. They may be considered as generating functions of single pedestrian traveling between l th and k th sites of periodic chain.
The substitution of (44) into (43) allows us to express the transition amplitude (42) through the Schur functions (3) and the Vandermonde determinants (4) [25]: where the summation is over (9). Substituting (45) into (41) and applying the Binet-Cauchy formula (14), we obtain [31]: where P L/n (y N , x N ) and F k; l (β) are defined by (14) and (44), respectively. At β = 0, the expression (47) transfers into (23). Notice, that for n = 0 the operatorΠ n = 1 and Eq. (47) gives the answer for the matrix element Ψ(v N ) | e −βH | Ψ(u N ) . Taking into account that where the eigen energy E N (θ N ) is given by (9), we obtain from (47) the answer for the persistence of ferromagnetic string (11), with the ground-state energy E N (θ g N ) given by (10). From the relation (46) follows the representation of the correlation function that we will use in the analysis of its asymptotical behaviour: where N (θ g N ) is the norm (20) of the ground state defined by (10), and P M/n (e −iθ N , e iθ g N ) is (14) on the solutions to Bethe equations (6).
The approach of the calculation of the persistence of ferromagnetic string used in this section allowed us to demonstrate the combinatorial nature of this correlation function. Naturally, the same answers could be obtained by insertion of the identity operator into the left-hand side of (41). In this way we shall calculate the persistence of domain wall in the next section.

Persistence of domain wall
To calculate the persistence of domain wall we insert the resolution of unity operator (21) into the numerator of (12) taken at arbitrary parametrization: The decomposition (52) transfers into (53) provided that (31) and (40) are used for each of form-factors in (52). The substitution of the equality (47), taken at n = 0, into (53) gives where D u N−n+1 ,u N−n+2 ,...,u N is given by (29); is defined analogously. After the differentiations the representation (54) takes the form: where A, B, C, D are the matrices with the entries: Finally, we obtain the answer for the persistence of domain wall (12): A (the same for B, C, and D). The explicit expression for the form-factor (32) allows us to express the persistence of the domain wall in terms of Schur functions starting with the relation (52): where the summation is over all solutions to the Bethe equation (6), and θ g N −n is the ground state solution of the system of N − n particles (13).

q-Binomial determinants and boxed plane partitions
Boxed plane partitions and q-binomial determinants are the important notions that will allow us to give the combinatorial interpretation of the asymptotical behavior of the correlation functions. Proposition 3 formulated in this section provides the determinantal formulas which enable the connection between the form-factors and enumeration of boxed plane partitions as well as of certain non-intersecting lattice paths.

q-Binomial determinants
The scalar product of the state-vectors (17), as well as the form-factors (23) and (38), are connected with the generating functions of boxed plane partitions (AI.1) and (AI.3) (see Appendix I). This connection takes place under special q-parametrization of the free variables u N and v N , and appropriate formulas are given by the statements of Proposition 3. Before turning to Proposition 3, we shall remind essential notions concerning the q-binomial determinants, [42].
For the arbitrary P and L ≤ N, these entries will take the form: This square matrix (T) 1≤j,k≤N consists of two blocks of the sizes L × N and (N − L) × N. When L = N, it consists of one block and is the matrix (18) under the q-parametrization. It seems appropriate to call the determinant, detT, as the Kuperberg-type determinant (see [43], where the problem of enumeration of alternating sign matrices has been investigated).
The q-binomial determinant a b q is defined by where a and b are ordered tuples: 0 ≤ a 1 < a 2 < · · · < a S and 0 ≤ b 1 < b 2 < · · · < b S .
The entries a j b i are the q-binomial coefficients (see Appendix II). In the limit q → 1, the q-binomial coefficients are replaced by the binomial coefficients a j b i . Then, the q-binomial determinant (58) is transformed to the binomial determinant: The binomial determinant (59) is non-negative and is positive provided b i ≤ a i , ∀i, [44]. Now we are ready to formulate the following Proposition 3 Let the square matrix (T) 1≤j,k≤N , consisting of two blocks of the sizes L × N and (N − L) × N, be defined by the entries (57) with P 2 < N < P . Then, the determinant of (T) 1≤j,k≤N is given by either of the following relations: where P ≡ P − N + 1, and Z q (L, N, P) is the generating function of plane partitions (AI.1) contained in a box B(L, N, P).
Proof Appendix II contains the proof of (60) and (61). The proof is based on the theory of the symmetric functions. The statements of Proposition 3 are valid for 1 ≤ L ≤ N. However, a formal relation can be written for L = 0 also: In this case, the q-binomial determinant is equal to q N 2 (P−1)P (its evaluation is in (AII.17)), and detT is nothing but the Vandermonde determinant.
The number of ways to travel from (0, 0) to (n, m) on a square lattice making elementary steps only to the north and to the east is equal to the binomial coefficient n + m m .
These ways are called the lattice paths. It was found in [44] that the binomial determinant (59) is equal to the number of self-avoiding lattice paths w 1 , w 2 , . . . , w S on a square lattice such that w i goes from In the considered case, Eq. (62), the binomial determinant is equal to the number of self-avoiding lattice paths starting at A i = (0, N +L+i−1) and terminating at B i = (L+i−1, L+i−1), 1 ≤ i ≤ P.
Because of the boundary conditions, this number of self-avoiding paths is equal to the number of self-avoiding paths starting at C i = (i, N + L + i − 1) and terminating at B i . The latter configurations are known as watermelons [33]. There exists bijection between watermelons and plane partitions confined in a box of finite size [36], and it provides the combinatorial proof of (62) (see Figure 2). Described lattice paths from (0, 0) to (n, m) contain in the rectangle n × m. If we place unit squares above and to the left of these lattice paths then the number of lattice paths is equal to the number of different ways to fit the Young diagrams with the largest part at most n and with at most m parts into an n × m rectangle. The q-binomial is the generating function of these Young diagrams (lattice paths), and each diagram comes with the weight q |λ| [45]. The weight of the lattice paths that terminate at B i = (L + i − 1, L + i − 1) and start at A i = (0, N + L + i − 1) or C i = (i, N + L + i − 1), respectively, differs on common factor q N 2 (P−1)P . Hence, the q-binomial determinant (60) is equal to the generating function of the plane partitions (AI.1) contained in B(L, N, P) multiplied on the common factor. The algebraic proof of this statement is given in Appendix II.
Due to Proposition 3, right-hand side of (63) is given by generating function of column strict plane partitions (61) with L and P replaced by N and M, respectively: and it coincides at q = 1 with the number of column strict partitions in B (N, N, M). Let us proceed with the form-factor of ferromagnetic string (23). Under the q-parametrization, the entries (24) are (57) with L = N and P = M − n. Therefore, due to Proposition 3, we may express the form-factor as the generating function of column strict plane partitions though in a smaller box B(N, N, M − n): The corresponding number of column strict plane partitions (AI.4) arises in the limit q → 1: Due to Propositions 1 and 2, the form-factor of the domain wall creation operator, Eq. (27), has the form in the q-parametrization: where the matrixT is given by (57) with L = N − n and P = M. The partitionsλ and λ are defined in (32). Applying (61), we obtain: Therefore, the form-factor is the generating function of plane partitions confined in the box B(N − n, N, M). In the limit we obtain the correspondent MacMahon formula (AI.2).

Low temperature asymptotics
Now let us turn to the main issue of the present paper -to the low temperature asymptotics of the correlation functions (11) and (12). We assume that our XX0 chain is long enough, M ≫ 1, while N is moderate: 1 ≪ N ≪ M. Besides, β in (11) and (12) is now inverse of the absolute temperature, β = 1 T (the Boltzmann constant is unity).

Persistence of ferromagnetic string
For large enough M, the summation over the solutions to the Bethe equations in the persistence of ferromagnetic string correlation function T (θ g N , n, β), Eq. (50), is replaced by integrations over the continuous variables. Under the same assumption, the elements of N-tuple θ g N of the ground state solutions (10) are such that cos θ g l ≃ 1. The approximate expression for T (θ g N , n, β) is of the form: (cos θ l −1) In the large β limit (low temperature limit), right-hand side of Eq. (70) can be re-expressed as: The power law decay in β of T (θ g N , n, β) (71) is governed by the critical exponent N 2 /2. The combinatorial factor P 2 M/n (1, 1) in A(N, n) (71) is, according to (66), the square of the number of column strict plane partitions A cspp (N, N, M − n).
The integral I N (72) is a special form of the Mehta integral related to the partition function of so-called Gaussian Unitary Ensemble [46,47]. Its value in terms of the gamma function [48] is known, and it is convenient for us to put it in the exponential form: In the considered limit the inverse of the square of the norm is equal to where ϕ N is defined in (73). Taking into account (73) and (74), we find out that the estimate (71) may be expressed as To study the asymptotical properties of the correlation functions in question, it is appropriate to represent ϕ N (73) through the Barnes G-function [49]: where γ is the Euler constant [48]. For every non-negative integer n G(n + 1) = (n!) n 1 1 2 2 . . . n n = n k=1 Γ(k) . (78) The theoretical aspects of the Barnes function are discussed in [50]. The equation (78) allows to re-express ϕ N (73): and thus the value of the integral (72) is given in terms of G-function: Asymptotically as z → ∞: where the constant A can be found in [49,50]. The behavior of ϕ N at N ≫ 1 results from (79) and (81): Hence, for the exponent (76) we obtain: where A is a constant. In order to estimate A cspp (N, N, M − n), we put P = M − n in (AI.4) and use (78): Taking into account (81), we find in the leading order: where B is a constant. Equation (85) gives us the asymptotical behavior of the number of column strict plane partitions in a high box with a square bottom B (N, N, M − n). Finally, taking into account (83) and (85), we put the asymptotic estimate (75) of the persistence of ferromagnetic string into the form: where C = AB 2 . If we assume that M and N are increasing, while the temperature T is decreasing, then from (86) it follows that T (θ g N , n, β) is decreasing provided that the restriction T < 1

Persistence of domain wall
Starting with the representation (55) for the persistence of domain wall correlation function, we can repeat the arguments of the previous subsection and deduce the following asymptotical expression: Furthermore, we find using (81): where D is some constant. Equation (88) Equation (89) enables us to state that F ( θ g N , n, β) is decreasing with increasing M and N provided that the temperature T is estimated analogously to the previous subsection.

Discussion
In our paper we have discussed the N-particle thermal correlation functions of the XX0 Heisenberg model on a cyclic chain of a finite length. We have considered the ferromagnetic string operatorΠ n (11) and the domain wall creation operatorF n (12). The calculations were based on the theory of the symmetric functions that allows us to express the answers in the determinantal form. The equations (23) and (38) for the form-factors which are the basic quantities in the above correlation functions are shown to be related to the generating functions of self-avoiding random walks and boxed plane partitions. The introduced q-binomial determinant plays an important role in the establishing of this relation.
The estimate of the asymptotical behavior of the persistence correlation functions of the operatorsΠ n andF n is done for low temperatures. The problem of calculating the asymptotical expressions leads to the calculation of the matrix integrals of the type of (70). These integrals were intensively studied in various fields of theoretical physics and mathematics, [51][52][53][54][55]. The low temperature approximation allows both to extract the combinatorial pre-factor and to reduce the matrix-type integrals to the partition function of the Gaussian Unitary Ensemble. Both correlations functions have a power law decay and have the same critical exponents, but their amplitudes are different: they depend on the squared number of plane partitions contained in a box of different size. These amplitudes are observable quantities. Expression for the corresponding exponent looks like the free energy appearing at small coupling for the large-N lattice gauge theory considered in [52]. This answer should be related to the third order phase transition [52], a possibility of which for the XX0 spin chain is discussed in [29]. The results obtained in [56] for the thermal correlation functions of the XXZ Heisenberg chain with the infinite coupling constant allows to argue that the third order phase transition is possible in this model as well.

Acknowledgement
This work was partially supported by RFBR (No. 13-01-00336) and by RAS program 'Mathematical methods in non-linear dynamics'. We are grateful to A. M. Vershik for useful discussions.

Appendix I
Here we provide some notions concerning boxed plane partitions and their generating functions while more details can be found in [14].
An array (π ij ) 1≤i,j of non-negative integers that are non-increasing as functions both of i and j is called plane partition π. The integers π ij are called the parts of the plane partition, and |π| = i,j π ij is its volume. Each plane partition has a three-dimensional diagram which can be interpreted as a stack of unit cubes (three-dimensional Young diagram). The height of stack with coordinates (i, j) is equal to π ij . It is said that the plane partition corresponds to a box of the size L × N × P provided that i ≤ L, j ≤ N and π ij ≤ P for all cubes of the Young diagram. If π ij > π i+1,j , i.e., if the parts of plane partition π are decaying along each column, then π is called the column strict plane partition.
We shall denote the box of the size L × N × P as the set of integer lattice points: An arbitrary plane partition π contained in B(N, N, P ) may be transferred into a column strict plane partition π cspp corresponding to B(N, N, P +N −1) by adding to π the N ×N matrix  which corresponds to a minimal column strict plane partition. The volumes of the column strict plane partition and correspondent plane partition are related: The generating function of plane partitions contained in B(L, N, P ) is equal to According to the classical MacMahon's formula, there are plane partitions contained in B(L, N, P ). It is clear that right-hand side of (AI.1) tends to A(L, N, P ) in the limit q → 1.
The generating function of the column strict plane partitions placed in B(N, N, P ) is equal to The limit q → 1 gives the number A cspp (N, N, P ) of the column strict partitions placed in B(N, N, P ): Γ(j) Γ(j + P + 1) Γ(j + N) Γ(j + P + 1 − N) , where the expression in terms of the gamma-functions [48] is appropriate to study the asymptotics. Notice that  In the limit q → 1, the q-binomial coefficient N r becomes the binomial coefficient N r .
Two analogues of the Pascal formula exist for the q-binomial coefficients The q-Vandermonde convolution for the q-binomial coefficients has the form The r th order elementary symmetric function e r = e r (x N ) of N variables, x N = (x 1 , x 2 , . . . , x N ), is defined by where 1 ≤ r ≤ N. The generating function identity can be used to define e r as the coefficients in the product: (1 + tx 1 )(1 + tx 2 ) . . . (1 + tx N ) = 1 + e 1 t + e 2 t 2 + · · · + e N t N .
It is convenient to introduce special notations for e r at x N = q N and x N = q N /q [14]: Consider the weakly decreasing partitions λ, P ≥ λ 1 ≥ λ 2 ≥ · · · ≥ λ N ≥ 0, where P ≡ P − N + 1. The conjugate partitionλ is the partition with the Young diagram consisting of columns of height λ i : N ≥λ 1 ≥λ 2 ≥ · · · ≥λ P ≥ 0. The corresponding strict partitionμ is given byμ =λ + δ P . The parts ofμ respect: P ≥μ 1 >μ 2 > · · · >μ P ≥ 0. The Schur functions (3) related to the partition λ are expressed through the elementary symmetric functions (AII.6) related to the conjugate partitionλ [14]: (AII. 8) In order to express the determinant of the matrix (57) as the q-binomial determinant, Eq. (60), we will use the statement of Proposition 2 under the q-parametrization (56): where the entries are: (AII.10) Let us denote the sum over partitions in (AII.9) by Σ S . Applying (AII.8) we bring it into the form: det eλ j −j+k (q N ) 1≤j,k≤P det eλ p−p+l (q L /q) 1≤l,p≤P , (AII.11) where summation runs over the conjugate partitionsλ, L ≥λ 1 ≥λ 2 ≥ · · · ≥λ P ≥ 0. In the explicit form, where the notations (AII.7) are taken into account. By definition, R 0 = L 0 = 1, while both R r and L r are equal zero for r > N. The Binet-Cauchy formula allows to express (AII.12) as the determinant of the P × P matrix:   where 0 ≤ s ≤ P − 1, we express Σ S as the determinant of the matrix with the entries given by the q-binomials To bring this determinant to the q-binomial form, Eq. (58), we will use the Pascal formulas (AII.4). As a first step, we combine the rows in (AII.15) in the following way: the P th is multiplied by q N +1 and the (P − 1) th is added to it, the (P − 1) th is multiplied by q N +1 and the (P − 2) th is added to it; . . . ; the 2 nd is multiplied by q N +1 and the 1 st is added to it. The second step is concerned with the matrix obtained: the P th row is multiplied by q N +2 and the (P − 1) th one is added to it, the (P − 1) th row is multiplied by q N +2 and the (P − 2) th one is added to it; . . . ; the 3 rd is multiplied by q N +2 and the 2 nd is added to it. After P − 1 steps, compensating the obtained factor by Since det Q = 1, the q-binomial determinant in (AII.17) is equal to q N 2 (P−1)P , and Σ S becomes the double product thus justifying the double product in (61).