The master T-operator for the Gaudin model and the KP hierarchy

Following the approach of [arXiv:1112.3310], we construct the master T -operator for the quantum Gaudin model with twisted boundary conditions and show that it satisfies the bilinear identity and Hirota equations for the classical KP hierarchy. We also characterize the class of solutions to the KP hierarchy that correspond to eigenvalues of the master T-operator and study dynamics of their zeros as functions of the spectral parameter. This implies a remarkable connection between the quantum Gaudin model and the classical Calogero-Moser system of particles.


Introduction
In this paper we discuss a correspondence between the following integrable systems: The link (i)-(ii) is a limiting case of the correspondence between quantum spin chains with the Yangian Y (gl(N ))-invariant rational R-matrices and the classical modified KP hierarchy based on the construction of the master T -operator [1,2]. The link (ii)-(iii) is a well-known story about dynamics of poles of rational solutions to soliton equations started by Airault, McKean and Moser [3] for the KdV equation, developed by Krichever [4] for the KP equation and extended to the whole KP hierarchy by Shiota [5]. The composition of (i)-(ii) and (ii)-(iii) implies the connection between the quantum Gaudin model and the classical CM model which is a limiting case of the connection between quantum spin chains and classical Ruijsenaars systems found in [1]. The link (i)-(iii) was also earlier established in [6] from a different reasoning.
The master T -operator was introduced in [1] 1 . It is a generating function for commuting transfer matrices of integrable vertex models and associated quantum spin chains which unifies the transfer matrices on all levels of the nested Bethe ansatz and Baxter's Q-operators in one commuting family.
It was also proven in [1] that the master T -operator, as a function of infinitely many auxiliary parameters ("times"), one of which being the quantum spectral parameter, satisfies the same hierarchy of bilinear Hirota equations as the classical tau-function does. This means that any eigenvalue of the master T -operator is the tau-function of a classical integrable hierarchy. For finite spin chains with Y (gl(N ))-invariant R-matrices this taufunction is polynomial in the quantum spectral parameter. The close connection of the spin chain spectral problem with integrable many-body systems of classical mechanics comes from dynamics of zeros of the polynomial tau-functions.
In this paper we obtain similar results for the quantum gl(N ) Gaudin model [8] with twisted boundary conditions and spins in the vector representation. It is known that this model can be treated as a certain limit of the integrable spin chain. However, the construction of higher integrals of motion emerging from the limiting procedure is not obvious. This makes the master T -operator for the Gaudin model a meaningful object interesting by itself.
The results of this paper can be outlined as follows. 1 A preliminary form of the master T -operator was previously discussed in [7].
• We construct commuting integrals of motion for the gl(N ) Gaudin model, with twisted boundary conditions and vector representations at the marked points in the quantum space, corresponding to arbitrary representations in the auxiliary space. They are presented in an explicit form using the matrix derivative operation. We also find functional relations satisfied by them.
• The master T -operator for the gl(N ) Gaudin model is the most general generating function for the commuting integrals of motion. It depends on an infinite number of auxiliary "time variables" t = {t 1 , t 2 , t 3 , . . .}, where t 1 can be identified with the spectral parameter x. The master T -operator is constructed explicitly using the matrix derivative.
• We show that the master T -operator satisfies the bilinear identity for the classical KP hierarchy and hence any of its eigenvalues is a KP tau-function. Here is a short dictionary of the Gaudin-KP correspondence:

Gaudin KP hierarchy
master T -operator ←→ τ -function spectral parameter ←→ the t 1 -variable higher transfer matrices ←→ Plücker coordinates Moreover, from the explicit form of the R-matrix and the Yang-Baxter equation it follows that this tau-function is a polynomial in x = t 1 . Therefore, according to [4,5], the dynamics of its roots in t i with i > 1 is given by equations of motion of the classical CM system of particles.
• The Gaudin-Calogero correspondence implies that the marked points x i in the Gaudin model (the inhomogeneities at the sites in the spin chain language) should be identified with initial coordinates of the CM particles while eigenvalues of the Gaudin Hamiltonians are their initial momenta. Eigenvalues of the Lax matrix for the CM model coincide with eigenvalues of the twist matrix (with certain multiplicities). Therefore, with fixed integrals of motion (action variables) in the CM model determined by invariants of the twist matrix, there are finite number of solutions for their values which correspond to different eigenstates of the Gaudin model. In other words, the eigenstates of the Gaudin Hamiltonians are in one-toone correspondence with (a finite number of) intersection points of two Lagrangian submanifolds in the phase space of the CM model. In short, the dictionary of the Gaudin-Calogero correspondence is as follows: Here is a more detailed summary of the Gaudin-Calogero correspondence. Let Φ be the 2n-dimensional phase space of the classical CM model with canonical coordinates {p 1 , . . . , p n , x 1 , . . . , x n }. As is known, the model is integrable and the complete set of integrals of motion in involution is given by H j = tr Y j (j = 1, . . . , n), where is the Lax matrix. In particular, −H 1 = i p i is the total momentum and is the CM Hamiltonian. Let L 1 ⊂ Φ be the Lagrangian hyperplane defined by fixing the x i 's and let L 2 ⊂ Φ be the Lagrangian submanifold defined as the level set of all the integrals of motion H j = N a=1 m a k j a with some fixed non-negative integer m a ∈ Z and real k a 's. Their intersection is a union of a finite number of points: L 1 ∩ L 2 = ∪ α∈I ψ α , with ψ α = (p Most of these results can be obtained by a limiting procedure from the corresponding results for quantum spin chains proved in [1]. The limit from spin chains to the Gaudin model is rather tricky, especially in the construction of higher integrals of motion, but in some cases it appears to be easier than the direct approach. In particular, the direct proof of the KP hierarchy for the master T -operator is more complicated for the Gaudin model because of a high degree degeneration of the latter. (Roughly speaking, it is technically easier to deal with a number of simple poles than with a multiple pole of high degree.) At the same time, along with the proof through the limit from spin chains (section 4.2.1 and Appendix C.1), we also give two direct proofs: one based on combinatorics of the symmetric group (Appendix C.2) and another one which exploits the matrix derivatives and matrix integrals technique (Appendix D).
2 The quantum Gaudin model as a limit of inhomogeneous spin chain Here x is the spectral parameter, η is an auxiliary parameter, I is the unity matrix. By e ab we denote the basis in the space of N×N matrices such that e ab has only one non-zero element (equal to 1) at the place ab: (e ab ) cd = δ ac δ bd . Note that P = ab e ab ⊗ e ba is the permutation matrix in the space C N ⊗ C N . A more general Y (gl(N ))-invariant R-matrix is e ab ⊗ π λ (e ba ) (2.2) which acts in the tensor product of the vector representation space C N and an arbitrary finite-dimensional irreducible representation π λ of the algebra U (gl(N )) with the highest weight λ. We identify λ with the Young diagram λ = (λ 1 , λ 2 , . . . , λ ) with = (λ) non-zero rows, where λ i ∈ Z + , λ 1 ≥ λ 2 ≥ . . . ≥ λ > 0. By e ab we denote the generators of the algebra U (gl(N )) with the commutation relations e ab e a b − e a b e ab = δ a b e ab − δ ab e a b . In this notation we have e ab = π (1) (e ab ), where π (1) is the N -dimensional vector representation corresponding to the 1-box diagram λ = (1).
Fix a matrix g ∈ GL(N ) called the twist matrix. For our purpose it is enough to consider diagonal twist matrices. A family of commuting operators (quantum transfer matrices or T -operators) can be constructed as where x i are arbitrary complex parameters. The trace is taken in the auxiliary space V λ where the representation π λ is realized. By R j0 λ (x) we denote the R-matrix (2.2) acting in the tensor product of the j-th local space C N of the chain and the space V λ (labeled by 0) 2 . The T -operators act in the physical Hilbert space of the model V = (C N ) ⊗n . It follows from the Yang-Baxter equation that the T -operators with the same η, g, x i commute for all x, λ and can be simultaneously diagonalized.
The normalization used above is such that T ∅ (x) = I ⊗n . It is convenient for the limit to the Gaudin model. Another useful normalization is In this normalization all T λ (x) and all their eigenvalues are polynomials in x of degree n.
At n = 0 the transfer matrix (2.3) is just the character: T λ (x) = χ λ (g), where χ λ (g) = tr π λ g is the character of the representation π λ given by the Schur polynomials s λ (y) of the variables y = {y 1 , y 2 , . . .}, y k = 1 k tr g k (the Jacobi-Trudi formula): with the complete symmetric polynomials h k (y) = s (k) (y) defined by Here the operator b i acts in the auxiliary space V λ labeled by 0.
where ξ(y, z) := k≥1 y k z k . It is convenient to set h k = 0 at k < 0. A "dual" form of (2.5) is where the elementary symmetric polynomials e k (y) Here and below λ is the transposed (reflected about the main diagonal) of the diagram λ. Let p 1 , . . . , p N be eigenvalues of g ∈ GL(N ) realized as an element of End (C N ). Then y k = 1 k (p k 1 + . . . + p k N ) and (see [9]). This formula implies that χ ∅ (g) = s ∅ (y) = 1.
A different but equivalent construction of quantum transfer matrices was suggested in [10]. It uses the special derivative operator on the group GL(N ) called there the coderivative operator. Let g be an element of the Lie group GL(N ) and f be any function of g. The (left) co-derivative is defined as 3 (2.10) The result of the action of D to a scalar function is a linear operator in C N , acting by D twice we get an operator in C N ⊗ C N and so on. In practice, we will specify in which space the operator is acting by adding a suffix, in the following way. Let V i ∼ = C N be several copies of C N , then, applying D to a scalar function k times we get an operator in In particular, applying the second D to (2.10), we can write and so on. If we work in the space (C N ) ⊗n , these operators can be realized as e cd ] = 0 for any i = j and any a, b, c, d ∈ {1, . . . , N }. 3 Originally, the co-derivative was defined [10] as This definition is equivalent to the one given here.
According to [10] the transfer matrix (2.3) can be represented as a sequence of coderivatives acting to the character: The Gaudin model can be obtained from here in the limit η → 0 provided that g is put equal to e ηh after taking the group derivatives. Here h is an element of the Lie algebra gl(N ) called the twist matrix for the Gaudin model. The standard way is as follows. Take λ = (1) (one box), then χ (1) (g) = tr g. Let us expand: We have D i tr g = g i , D j g i = P ij g i . Therefore, if g = e ηh , then and the η-expansion of the transfer matrix T (1) (x) reads: The O(1) and O(η) terms are central, so the operators at the order η 2 form a commutative family where are the commuting Gaudin Hamiltonians: [H i , H j ] = 0. Below we assume that h is diagonal, h = diag (k 1 , . . . , k N ), then the Gaudin Hamiltonians can be written more explicitly in the form be the decomposition of the Hilbert space of the Gaudin model V into the direct sum of eigenspaces for the operators M a with the eigenvalues m a ∈ Z ≥0 , a = 1, . . . , N . Then the eigenstates of H i 's are in the spaces V({m a }). Because a e aa = I is the unit matrix, a M a = nI ⊗n , and hence N a=1 m a = n. (2.17) Note also that

Commuting operators for the Gaudin model
A general algebraic construction of higher members of the Gaudin commutative family was proposed in [11,12]. In principle, their explicit form can be found from the ηexpansion of the spin chain transfer matrices T λ (x) for general λ. But there is a problem how to extract non-trivial integrals of motion from this expansion. For fundamental representations (one-column diagrams) this problem was solved by Talalaev [13]. Basically, Talalaev's idea was to consider special linear combinations of the transfer matrices T λ (x) such that their η → 0 limits start from higher degrees of η. The coefficients in front of the leading terms as η → 0 commute and can be regarded as higher transfer matrices (generating functions of integrals of motion) of the Gaudin model. Here we extend this procedure to transfer matrices associated with diagrams of arbitrary shape.
Let us modify the definition of the transfer matrices (2.12) as follows: Here g is put equal to e ηh after taking all the derivatives. To avoid misunderstanding, we stress thatT λ (x) is not the T -operator for the twist matrixg = g − I because the co-derivatives act to g rather thang. Since χ λ (g − I) is a linear combination of characters χ µ (g) with different µ and g-independent coefficients,T λ (x) is a linear combination of the T µ (x)'s. More precisely, using (2.9) and the Cauchy-Binet formula, one can prove that where is the binomial coefficient (see [9, page 47, example 10]). The sum is taken over all Young diagrams µ that are contained in λ including the empty diagram and λ itself. Therefore, In particular,T which agrees with Talalaev's prescription for one-column diagrams. As we shall see, the leading behavior ofT λ (x) as η → 0 is O(η |λ| ), so we can define the higher Gaudin transfer matrices as To represent them in a more explicit form, we need to modify the definition of the co-derivative by passing from the group derivative to the Lie algebra derivative. Let h be an element of the Lie algebra gl(N ) and f be any function of h. We define In fact this is the usual matrix derivative used in the theory of matrix models. For example: Similarly to (2.11), we have and so on. We remark that our modified co-derivative is commutative (d 1 d 2 = d 2 d 1 ), although the co-derivative for the original spin chain is non-commutative (D 1 D 2 = D 2 D 1 ).
The following lemma is crucial for our construction. Sketch of proof. We illustrate the idea of the proof by the example k = 2. The general case is proved in the same way but requires too bulky formulas. Using (2.11), we can write .
After the rescaling ε i → ηε i the leading η → 0 term can be written in the form The homogeneity of the function f and definition (3.8) then imply that The character χ λ (g) is a homogeneous function of degree |λ|. Therefore, by virtue of the lemma, the family of commuting operators for the (twisted) Gaudin model can be constructed as follows: (3.10) or, in the polynomial normalization, For example: where H(x) is given by (2.13). Note that in this approach H(x) emerges from T (1 2 ) (x) or T (2) (x) rather than T (1) (x).
Hence it is clear that e x tr h T G (x, t) depends on x, t 1 only through their sum x + t 1 . In particular, we will use the relation The expansion in the Schur functions is The T -operators T G λ (x) can be restored from the master T -operator according to the formula In particular, With a given z ∈ C, the following special shift of the time variables is often used. As we shall see below, T G (x, t ± [z −1 ]) regarded as functions of z with fixed t contain an important information. Here we only note that equation (4.4) implies that T G (x, 0 ± [z −1 ]) is the generating series for T -operators corresponding to the one-row and one-column diagrams: We are going to prove that the master T -operator (4.1) satisfies the bilinear identity for the KP hierarchy [14,15] which states that Here ξ(t, z) := k≥1 t k z k and the integration contour is chosen in such a way that all singularities coming from the T G 's are inside it while those coming from e ξ(t−t ,z) are outside it.
Remark. The bilinear identity is invariant under the change of signs of all times: if and taking the residues, one obtains from it the 3-term Hirota equation [15,16,17,18] (the Fay identity) (4.8) Tending z 0 → ∞, one obtains a simpler looking 3-term equation (4.9) It appears to be equivalent to its z 3 → ∞ limit (the differential Fay identity) 4 (4.10) In fact it was proved in [19,20] that all the Hirota equations of the form (4.8), (4.9), (4.10) are equivalent to each other and to the bilinear identity (4.7).
This means that each eigenvalue of the master T -operator is a tau-function of the KP hierarchy. Equation (4.3) is the expansion of the tau-function in Schur polynomials [16,21,22].

The quantum Giambelli formula
All the bilinear equations for T G (x, t) follow from the quantum Giambelli formula for the Gaudin model: where the notation T G l,k (x) := T G (l+1,1 k ) (x) is used and d(λ) is the number of boxes in the main diagonal of the Young diagram λ. Note that the quantum Giambelli formula for the Gaudin model has the same form as the one for the original spin chain 5 before the η → 0 limit. In the polynomial normalization we have ) and sum such equations for the pairs (z 1 , z 2 ), (z 2 , z 3 ) and (z 3 , z 1 ). The result is equation (4.9). 5 The quantum Giambelli formula for U q (B n )-invariant vertex models was proposed in [23].
If (4.12) holds, then the Jacobi identity for this determinant produces the 3-term bilinear identities (the Plücker relations) for the coefficients of the Schur function expansion (4.3). This implies the KP hierarchy for T G (x, t) (see [15, example 2, page 959]). A more direct proof is given in the appendix. Set The proof of the quantum Giambelli formula is based on the following lemma: The operator Q(z, ζ) obeys the "exchange relation": This is our basic identity. It is proved in the appendix. We give two different proofs: one through the limit from spin chain (Appendix C.1) and a direct proof (Appendix C.2). Note that (4.14) is equivalent to the very special case of the Fay identity (4.8) at . . = x n (in this case the third term vanishes). We will use the notation T G,n λ (x) for the T -operator acting on n sites. Proof of the corollary. We will write T G,n α,β (x) := T G,n α,β for brevity. It is enough to prove that d n+1 T G,n α 1 β 1 d n+1 T G,n At n = 0 T G,0 αβ is just the hook character χ α,β := χ (α+1,1 β ) (h) and the assertion can be easily proved by passing to the generating function of hook characters [9] E(z, ζ) = α,β≥0 Let us multiply the determinant by z α 1 1 z α 2 2 (−ζ 1 ) β 1 (−ζ 2 ) β 2 and sum over all α i , β k . We then see that hence the statement follows. More generally, the assertion of Lemma 4.1 means that 6 Now we claim that vanishing of the determinant for any n ≥ 0 is equivalent to vanishing of for any n ≥ 1. This is clear from the identity and the "initial condition" D 0 = 0 established above. In a similar way, one can show that the assertion is equivalent to vanishing of for any n ≥ 0. But this follows from (4.16) and the corollary is thus proved.
After these preliminaries, the proof of the quantum Giambelli formula is easy. Suppose that it holds for some n ≥ 0 (for example, it holds for n = 0, in which case it is the usual Giambelli formula for characters; see, e.g., [9]). Let us apply 1 to the both sides. In the l.h.s. we get T G,n+1 λ (x) -the T -operator for the model on n + 1 sites. In the r.h.s. we get det 1≤i,j≤d(λ) Using the rule of differentiating determinants and the fact that the matrix being a submatrix of the (d n+1 T G,n α,β (x)) α,β≥0 , has rank 1, we see that the r.h.s. is equal to This proves the quantum Giambelli formula for the model on n + 1 sites.

The master T -operator and matrix derivatives
Another proof of the bilinear identity (4.7) can be given using the technique of matrix derivatives and matrix integrals.
Let us first consider the special case of the master T -operator for the Gaudin model (4.1) at x = x 1 = x 2 = . . . = x n : The operator d is just the matrix derivative with respect to the transposed matrix h t , which is well known in the theory of matrix models (for more details see, e.g., [24]): One can introduce the generating function for the master T -operators which depends on an auxiliary external matrix A: Clearly, this function generates the master T -operators for any n: (4.19) injn . In a similar way, one can introduce a generating function for products of the master It has the form This expression can be regarded as a formal series in A with coefficients which are functions of t and t . Equivalently, it can be regarded as a formal series in t and t with coefficients depending on A.
We have, for example: In general, the product (4.20) is just the n'th matrix derivative of this function with respect to the (transposed) matrix A: Therefore, the bilinear identity (4.7) for (4.17) and any n is equivalent to the following scalar identity: The proof of this identity, as well as its generalization to arbitrary x and x i 's, is given in Appendix D.
The CBR determinant formulas for the original spin chain are: They are spectral parameter dependent analogs of the Jacobi-Trudi formula for characters of the group element g (the twist matrix). In fact, (4.26) and (4.27) reduce to (2.5) and (2.7) in the limit |x| → ∞.
We have found the following analogs of the CBR formulas for the Gaudin transfer matrices: In the same ways as (4.26), equation (4.28) is also an analogue of the Jacobi-Trudi formula (2.5) and reduces to it in the limit |x| → ∞. Equation (4.29) is a "dual" determinant formula (an analogue of (2.7)).
The determinant formulas (4.28) and (4.29) follow from the differential Fay identity 7 . A sketch of proof is given in Appendix B. Here we present an auxiliary determinant formula obtained by iterations of the differential Fay identity.
For any positive integer N and any subset Solving (4.10) for (4.30) recursively, we obtain the determinant formula Note that (4.9) is a Plücker identity and (4.10) is the Jacobi identity for minors of the ,j≤m entering (4.31). As is shown in Appendix B, (4.31) is a generating function for (4.28).

The Baker-Akhiezer functions
According to the general scheme, the Baker-Akhiezer (BA) function and its adjoint corresponding to the tau-function (4.1) are given by the formulas For brevity, we will refer to both ψ and ψ * as BA functions. In terms of the BA functions, the bilinear identity (4.7) can be written as Using the definition (4.1), we have: 36) 7 We also remark that the quantum Giambelli formula (4.11) follows from (4.28).
where we write simply det(z − h) instead of det(Iz − h). From these formulas we see that e −xz−ξ(t,z) ψ(x, t; z) is a polynomial in z −1 of degree N while e xz+ξ(t,z) ψ * (x, t; z) is a rational function of z with poles at the points z = k i (eigenvalues of the matrix h) of at least first order (because of the det(z − h) in the denominator). Moreover, since each co-derivative raises the order of the poles, these poles may be actually of a higher order, up to n + 1. Also, as is seen from the second formula, this function has the N -th order zero at z = 0. (We assume that k i = 0.) Regarded as functions of x, both e −xz ψ(x, t; z) and e xz ψ * (x, t; z) are rational functions of x with n zeros and n poles which are simple in general position. From (4.1) and (4.35), (4.36) it follows that The BA functions for any solution of the KP hierarchy satisfy the following differential equations: We also note the formulas for the stationary BA functions ψ(x, z) := ψ(x, 0; z), ψ * (x, z) := ψ * (x, 0; z) which directly follow from (4.35), (4.36): Below we will also need the relation ) This relation follows from the bilinear identity (4.7) by applying ∂ tm and putting t k = t k afterwards.

Zeros of the master T -operator as the Calogero-Moser particles
The eigenvalues of the master T -operator are polynomials in the spectral parameter x: The roots of each eigenvalue have their own dynamics in the times t k . This dynamics is known [4] to be given by the rational CM model [27]. The inhomogeneity parameters of the Gaudin model play the role of coordinates of the CM particles at t i = 0: . Comparing with the third equation in (3.12), we conclude that the initial velocities are proportional to eigenvalues of the Gaudin Hamiltonians: This unexpected connection between the quantum Gaudin model and the classical CM model was observed in [6] using a different reasoning.

The Lax pair for the CM model
Following Krichever's work [4], let us derive equations of motion for the t 2 -dynamics of the x i 's. It is convenient to denote t 2 = t and put all other times to zero since they are irrelevant for this derivation. From (4.41) we see that The method of [4] is to perform the pole expansion of the linear problem (4.39) for the BA function ψ. From (4.35) we have (4.37)). Plugging this into (4.39), we obtain The l.h.s. is a rational function of x with first and second order poles at x = x i (possible poles of the third order cancel automatically) vanishing at infinity. Therefore, to solve the linear problem it is enough to cancel all the poles. Equating the coefficients in front of each pole to zero, we get the following system of equations for i = 1, . . . , n: These equations can be rewritten in the matrix form: The compatibility condition of the system (5.5) iṡ It is the Lax representation for the equations of motion of the CM model: The matrices Y, T form the Lax pair for the CM model.
In a similar way, plugging the adjoint BA function to the adjoint linear problem (4.40), we get where c * = (c * 1 , c * 2 , . . . , c * n ) t . Using (5.5), (5.11) and recalling that c 0 (z) = z −N det(z − h), we find the solutions for the vectors c, c * : . For the functions ψ, ψ * themselves we then have: Here we list some properties of the matrices X, Y, T to be used in what follows.
As is well known (and easy to check), the matrices X, Y satisfy the commutation relation [X, Y ] = I − 1 ⊗ 1 t (5.14) (here 1 ⊗ 1 t is the n×n matrix of rank 1 with all entries equal to 1).
but the last trace is 0 as trace of a commutator.
Let E ik be the basis matrices of gl(n) having just one non-zero element (equal to 1) at the place ik: (E ik ) i k = δ ii δ kk . The following relations are easy to verify directly:

Integrals of motion
The matrix Y is the Lax matrix for the CM model. As is seen from (5.8), the time evolution preserves its spectrum, i.e., the coefficients J k of the characteristic polynomial are integrals of motion. The highest integral, J n , was found explicitly in [28], where a recurrence procedure for finding all other integrals of motion was also suggested. In fact this procedure is equivalent to the following explicit expression for the characteristic polynomial: Note that this expression is well-defined because the sum obtained after expansion of the exponential function in the r.h.s. contains a finite number of non-zero terms.

Eigenvalues of the Lax matrix
The singularities of the c(z, t), c * (z, t) as functions of z are the same as singularities of the functions ψ, ψ * in the finite part of the complex plane. From (4.35) we see that c(z, t) has the pole of order N at z = 0 and no other poles. At the same time the first equation in (5.12) states that there are possible poles at eigenvalues of the matrix Y (t) (which do not depend on time). Therefore, they must be canceled by zeros of det(z − h) which are at z = k i and are assumed to be simple. If all eigenvalues of Y are distinct, such cancellation is possible only if n ≤ N . However, the most interesting setting for the quantum spin chains and for the Gaudin model in particular is quite opposite: n > N or even n N (large chain length at a fixed rank of the symmetry algebra). Therefore, we conclude that in this case On the first glance, a multiple eigenvalue k i might lead to an unwanted pole of ψ at z = k i coming from the higher order pole of the matrix (z − Y ) −1 which now cannot be compensated by the simple zero of det(z − h). In fact the higher order poles do not appear in the vector (z − Y ) −1 1 because 1 is a special vector for the matrix Y which can be decomposed into N Jordan blocks of sizes m i ×m i . However, they do appear in the vector (z − Y t ) −1 1 and the function ψ * has multiple poles at z = k i (with multiplicities m i + 1).
Another way to see that eigenvalues of the Lax matrix Y are the same as eigenvalues of the twist matrix h (with appropriate multiplicities) is to compare expansions of (4.42) and (5.13) at large x. From (4.42) we have: The expansion of (5.13) at t = 0 gives (the lemma from the previous subsection is to be used): where we set Y 0 := Y (0). Therefore, we conclude that where M a is the operator (2.16). Hence we see that the M a is the "operator multiplicity" of the eigenvalue k a . In the sector V({m a }) the multiplicity becomes equal to m a .

Equations of motion in the hamiltonian form
The hamiltonian form of equations of motion is with the Hamiltonian This result was generalized to the whole hierarchy in [5]: The H k 's are higher integrals of motion (Hamiltonians) for the CM model. They are known to be in involution [28,29,30]. This agrees with commutativity of the KP flows.
The integrals H k are connected with the integrals J k introduces in (5.17) by Newton's formula [9] n k=0 J n−k H k = 0 (5.23) For completeness, we give a short derivation of (5.22) which is a version of the argument from [5]. The main technical tool is equation (4.44) which states that Matching coefficients in front of the poles at x i , we get Using (5.12), we have: Equations (5.16) and (5.15) allow us to rewrite the result as follows: This proves the first equality in (5.22). Next, applying ∂ t 2 to (5.24), we get: We continue the chain of equalities using the second formula in (5.16): This is the second equality in (5.22).

Determinant formula for the master T -operator
The results of [5] imply an explicit determinant representation of the tau-function. It is easy to adopt it for the master T -operator T G (x, t) (5.1). Let X 0 = X(0) be the diagonal matrix X 0 = diag(x 1 , x 2 , . . . , x n ), where x i = x i (0) and Y 0 be the Lax matrix (5.6) at t = 0, with the diagonal elements being the Gaudin Hamiltonians H i = −p i (0) (see (5.2)). Then Plugging this into (4.32), (4.33) we find formulas for the stationary BA functions: Let us show that these formulas are equivalent to the stationary versions of (5.13). Using commutation relation (5.14), we have: which shows that (5.13) and (5.26), (5.27) are indeed equivalent at t = 0.

Spectrum of the Gaudin Hamiltonians from the classical CM model
It follows from the above arguments that eigenvalues of the Gaudin Hamiltonians H i , i = 1, . . . , n (2.15), can be found in the framework of the classical CM system with n particles as follows. Consider the matrix i.e., given the initial coordinates x i and the action variables H j = tr Y j 0 one has to find possible values of the initial momenta p i = −H i . Taking into account equations (5.17) and (5.18), we can represent the equations for H i in the form of the equality which has to be satisfied identically in z. This identity is equivalent to n algebraic equations for n quantities H 1 , . . . , H n .
In other words, the eigenstates of the Gaudin Hamiltonians correspond to the intersection points of two Lagrangian submanifolds: one obtained by fixing the x i 's and the other obtained by fixing the H i 's, with values of the latter being determined by eigenvalues of the twist matrix. This purely classical prescription appears to be equivalent to the Bethe ansatz solution and solves the spectral problem for the quantum Gaudin Hamiltonians.

Concluding remarks
We have shown that the most general generating function of commuting T -operators in the Gaudin model (the master T -operator) satisfies bilinear equations of the classical KP hierarchy. This implies that each eigenvalue of the master T -operator is a classical taufunction. By construction, these tau-functions appear to be polynomials in the spectral parameter x of the Gaudin model which can be identified with the KP time t 1 . The dynamics of zeros of polynomial tau-functions leads to a close connection with integrable many-body problems.
This result immediately leads to the important conclusion that the eigenstates of the Gaudin Hamiltonians are naturally labeled by intersection points of two Lagrangian submanifolds in the phase space of the classical Calogero-Moser system of particles. This is a degenerate case of a more general correspondence between quantum spin chains (of the XXX and XXZ types) and the classical Ruijsenaars-Schneider model outlined earlier in [1,2]. Presumably, this "quantum-classical" correspondence extends to models with elliptic R-matrices and to supersymmetric integrable models.
Recently, the link between quantum spin chains and integrable many-body problems of classical mechanics has been discussed [31,32] in the context of supersymmetric gauge theories. Physical consequences of this "quantum-classical correspondence" are yet to be recognized and articulated while its mathematical roots lie deeply in quantum theory of Hitchin's integrable systems [33,31].
In the main text the Planck constant in the Gaudin model was set to 1. Here we would like to remark that if one introduces the Planck constanth in the Gaudin model, then the master T -operator satisfies theh-version of the KP hierarchy [19] which is obtained by the transformation t k →h −1 t k for k ≥ 1 and x →h −1 x. The coupling constant of the CM model (the coefficient 2 in the numerator of the second term in the r.h.s. of (5.21)) becomes 2h 2 .
Because of the space limitation, in this paper we have not addressed the construction of Baxter's Q-operators for the Gaudin model. In fact the master T -operator unifies them in one commuting family with the higher T -operators. Namely, taking residues of T G (x, t + [z]) at the poles (at the points k i which are eigenvalues of the twist matrix) one can introduce a family of 2 N Q-operators. The bilinear identities among them define the Bäcklund transformations from which the Bethe equations follow. For the original spin chains, the Bäcklund transformations were discussed in [34,35,36,37] on the level of the eigenvalues of the transfer matrices, and in [7,1,2] on the level of operators. We plan to extend this approach to the Gaudin model in a separate publication.

Note added
After the arXiv version of the current paper [arXiv:1306.1111] appeared, a paper [A. Gorsky, A. Zabrodin, A. Zotov, Spectrum of Quantum Transfer Matrices via Classical Many-Body Systems, arXiv:1310.6958] appeared on the arXiv. They clarified the relationship between inhomogeneous quantum spin chains and classical integrable many-body systems.

Appendix: Proofs and details
In this appendix, we present some details omitted in the main text. Some of the facts related to the KP hierarchy are consequences of the general theory known in the literature (see, e.g., [14,15,16,19,22]). However, to make the text self-contained, we give here some details of the calculations in our notations.

A KP hierarchy from Giambelli relation
Here we show that if the coefficients c λ of the expansion τ (t) = λ c λ s λ (t) satisfy the where c (l|k) := c (l+1,1 k ) are the coefficients for the hooks, then τ (t) is a KP tau-function, i.e., it solves the bilinear identity for the KP hierarchy. We assume c ∅ = 0.
Using this notation, we can represent (A1) and the Giambelli identity for the Schur functions as d × d determinants: c λ = det 1≤i,j≤d(λ) c (α i |β j ) , s λ (t) = det 1≤i,j≤d(λ) s (α i |β j ) (t). Plugging these into the expansion and separating the contribution of the empty diagram, we get: The sums over β j in each term can be handled with the help of the Cauchy-Binet formula: Next, it is not difficult to notice that the sums over α i and d represent the expansion of the determinant det ij (δ ij + A ij ) in terms of diagonal minors of the matrix A ij . Therefore, we have To proceed, we need the identity [9] k≥0 It means, in particular, that the lower-triangular matrices (with 1's on the main diagonal) h i−j (t) and h i−j (−t) are inverse to each other. (Recall that h k (t) = 0 at k < 0.) This is obvious from the fact that the product of their generating functions is 1. Using this identity, we rewrite (A2) as follows: Finally, we conclude that equation (A2) can be represented in the form where the Z×Z ≥0 matrix S is given by One can see that this formula is the general solution to the KP hierarchy given in the form of determinant of a semi-infinite matrix [14,16]. The unusual argument of h i−l (−t instead of t) does not spoil anything because the KP hierarchy is invariant under changing signs of all times (see the remark after equation (4.7)). Therefore, (A2) provides the general solution to the KP hierarchy and thus solves the bilinear identity for any c (j|k) .

B Details on functional relations for higher Gaudin T -operators
From the determinant (4.31) to CBR formula (4.28) The equivalence of (4.31) and (4.28) can be shown in the following way 8 . First, let us rewrite (4.31) as follows 8 See also discussions on the generation function of the transfer matrices in [7,37].
Using (2.6) and expanding the determinant in both side of (B1), we obtain where S m is the permutation group on {1, 2, . . . , m} and sgn(σ) is the signature of σ ∈ S m . Then, comparing the coefficient of in the both side of (B2), we get where we used (2.5). This reduces to (4.28) at t = 0 after a renormalization 9 and the identification m = (λ).
The Jacobi identity (B11) for the matrix for j 1 = k 1 = m − 1 and j 2 = k 2 = m corresponds to (B8). In fact we obtain the relations where I = {i 1 , i 2 , . . . , i m−2 } and ∆({i 1 , i 2 , . . . , i m }) = 1≤a<b≤m (z −1 ia − z −1 i b ). Substituting the above relations to (B11), we obtain (B8). Among (B12)-(B17), equations (B14) and (B17) are rather non-trivial. Equations (B14) and (B17) can be proved using the Leibnitz rule for the derivative with respect to x, taking linear combinations of columns in determinants and using the following identity valid for any matrix (A jk ) 1≤k,j≤m and parameters z k : In a similar way, one can show that the Plücker identity (B9) for the matrix M = z j−m+2  C Two proofs of Lemma 4.1

C.1 Proof by means of the limit from spin chain
The proof is based on the limit from the spin chain, where the relation follows from the Hirota equations for the master T -operator (or from the quantum Giambelli formula for the transfer matrices) proved in [1].
Applying the general formula (3.3) for the hook characters χ α,β (g) = χ (α+1,1 β ) (g), we get Alternatively, this expansion can be derived by comparison of the generating functions. The generating function of the hook characters is It is easy to obtain the following relation between generating functions for the characters χ α,β (g) and χ α,β (g − I):

This implies the relation for characters (C2).
The important thing is that the coefficient in front of χ α ,β (g) in (C2) factorizes into a product of α and β -dependent factors. Using Lemma 3.1 we can write Substituting (C2) into (C3) and taking note on the relation (C1), we see that which is (4.16).

C.2 Direct proof
We will now provide an alternative proof of Lemma 4.1 which does not rely on the results proven in [10,7,1], although it uses some ideas similar in spirit to the ones used in these papers. The same argument will also prove the commutation of the operators Q(z, ζ). More precisely, we will show that The S z case is the assertion of Lemma 4.1. The S z,ζ case is the commutativity of these operators that can also be proven from the Yang-Baxter equation. Our argument is universal and works for both cases, thus, later on we denote both S z and S z,ζ by the same symbol S.
In order to prove (C5), one can first notice that where S n is the permutation group over the set {1, 2, . . . , n}; and P σ denotes the permutation operator associated to the permutation σ, i.e. the operator such that for any For instance, for the transposition τ (i,j) (which permutes indices i and j, τ (i,j) (i) = j, τ (i,j) (j) = i, and does not act on all other indices, τ (i,j) (k) = k if k = i, j), we have P τ (i,j) = P i,j . The expression (C8) is easily proven by recurrence over the number of spins, by noticing that in d n d ⊗(n−1) w(z), the derivative d n can (due to the Leibnitz rule for derivatives) either act on w(z) to give (d ⊗(n−1) w(z)) ⊗ z 1−hz which corresponds, in equation (C8), to the permutations such that σ(n) = n, or it can act on any factor 10 z to produce σ∈S n−1 P i,n P σ z 1−hz ⊗n , which produces, in equation (C8), all the terms where σ(n) = n.
Moreover, for arbitrary numbers α 1 , α 2 , . . . α m and z 1 , z 2 , . . . z m the Leibnitz rule for derivatives of products gives (C9) 10 We remind that h (i) denotes the operator I ⊗i−1 ⊗ h ⊗ I ⊗n−i , so that z Noticing that for any operators A, B, the tensor product A ⊗ B can be written as (A ⊗ I) (I ⊗ B), we can also write this last equality as In what follows we will often use this notation. For instance, it is convenient for the following generalization of (C9): where C(σ) denotes the set of the cycles 11 of the permutation σ: for instance if n = 2, then the identity permutation has two cycles C(1) = {{1}, {2}}, whereas the cyclic permutation has only one cycle C(τ ) = {{1, 2}}. The expression (C10) can be proven by recurrence by the same argument as (C8).
Thus, (C10) allows one to express Q(z, ζ) as As an example of manipulations with this operator, one can check that for arbitrary k, l ∈ {1, 2, . . . n} we have Q(z, ζ) = P k,l Q(z, ζ)P k,l : for each σ, we have P k,l P σ P k,l = P τ (k,l) •σ•τ (k,l) , where τ (k,l) is defined above as the transposition which exchanges k and l, and where • denotes the composition of permutations. Moreover, σ → τ (k,l) • σ • τ (k,l) is a bijection of S n , and the set {i 1 , i 2 , · · · , i m } is a cycle of τ (k,l) • σ • τ (k,l) if and only if {τ (k,l) (i 1 ), τ (k,l) (i 2 ), · · · , τ (k,l) (i k )} is a cycle of σ. Hence the relation Q(z, ζ) = P k,l Q(z, ζ)P k,l is clear as soon as we notice that h (i) P k,l = P k,l h (τ (k,l) (i)) . This relation, which can also be obtained directly from the commutation relations d i d j = d j d i , also implies by recurrence that For notational simplicity, we introduce the operators Relation (C5) is then equivalent to where S(O 1 O 2 ) is a polynomial of degree at most two in each h (i) . Of course, one can note that due to relations like h (i) P i,j = P i,j h (j) , the degree in each individual h (i) is well defined only if the position of the operators P σ is specified. We will define degrees by putting all permutation operators to the left of all operators h (i) .
We will first see that the term with degree 0 in all of the h (i) 's vanishes, then that the terms having degree 2 in any of the h (i) 's vanish, and finally that the terms having degree 1 in any of the h (i) 's also vanish. To this end we will proceed by recurrence assuming that n ≥ 2 and that for all 1 ≤ n < n, we have S(O where #c denotes the number of elements in c. Interestingly enough, one can show from the expression [38,39] of the characters of the representations of S n , that (C16) can be rewritten as where χ α,β (σ) denotes the character of the permutation σ in the hook representation (α + 1, 1 β ) of the symmetric group S α+β+1 . Indeed, the Murnaghan-Nakayama rule [38,39] says that for σ ∈ S n , the character χ α,β (σ) is obtained by summing contributions from each Young tableau like the following one (which contributes to χ α,β (σ) where α = 7, β = 4, and σ is for instance the permutation 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 → 1, 4, 2, 3, 5, 7, 8, 9, 10, 6, 12, 11)  i.e. for each tableau such that each number i appears c i times (where the permutation σ has l cycles of respective size c 1 , c 2 , . . . , c l ), in such a way the these numbers increase in each row or column, and that for each number i, the set of boxes having the label i is convex 13 . The contribution of each such tableau is l i=1 (−1) h i −1 where h i is the number of rows where the label i is present at least once. Hence (by denoting as I the set of all labels i ≥ 2 present in the first column) one gets where |I| denotes the number of elements in I.
where the equality is obtained by noticing that in the l.h.s., each term where I 1 cancels with the term where I is substituted with I \ {1}. This proves that the relations (C16) and (C17) are identical.
One can then note that up to a numerical factor a α,β , each combination σ P σ χ α,β (σ) is a projector onto the space associated to the irreducible representation (α + 1, 1 β ) of S α+β+1 . Hence the orthogonality of these spaces implies the relation 14 This allows us to compute which clearly shows that S(O 1 O 2 | h=0 ) = 0, which is the vanishing of the term with degree 0 in all the h (i) 's.
Terms of degree 2. In order to understand the terms of higher degree in h (i) of the product O 1 O 2 , we will first investigate the terms of higher degree in h (i) of the operator 14 One way to show (C21) is by introducing the Young symmetrizer c λ associated to a representation λ, which obey λ ≡ σ P σ χ λ (σ) ∝ σ P σ c λ P σ − 1 (see for instance (3.11) on page 126 of [40]). In the case λ = λ , the Lemma (4.3.H) on page 126 of [40] says that λ λ = 0, whereas the relation written as = right before (3.11) on page 126 of [40] is the statement that if λ = λ , then λ λ ∝ λ .
O m : for instance, if we denote by (O m ) (h (n) ) 1 the coefficient 15 whereσ ∈ S n−1 is defined bŷ where τ (i,j) is the transposition of i ↔ j defined previously. The first line in (C23) comes from the fact that if σ(n) = n, then {n} is a cycle of σ and h (n) can only appear through the factor i∈{n} z where we used (C12) and (C25). Using the recurrence hypothesis, which contains the statement that S(O 2 ) (which we will denote as (S(O 2 )) (h (n) ) 2 ) does vanish.
Some notations. The analysis of the terms with degree 1 requires a couple of notations that we will now introduce. We have already seen that the coefficient of O m with degree 1 in h (j) is (O m ) (h (j) ) 1 = (−ζ m z m P j,n n−1 k=1 P k,n ) (O [n−1] m ⊗ I) P j,n , where one can notice 15 In this notation, (. . . ) (h (k) ) 1 denotes the coefficient of degree one in h (k) , where the ambiguity in the degree of each individual h (k) (due to P i,j h (i) = h (j) P i,j ) is fixed as explained above: by moving all permutation operators to the left of all h (i) . that n−1 k=1 P j,n P k,n = k∈{1,2,...,n}\j P k,j P j,n . This allows to rewrite this expression as (−ζ m z m k =j P k,j ) (O (C27) One can also notice, using the same arguments as in the proof of (C12), that it is equal w(ζm) . We will also use natural generalizations of this notation, for instance Another important remark is that all the operators that we deal with are of the form σ P σ c σ , where c σ is a coefficient which only contains some scalars and the operators h (j) . We denote this coefficients c σ as [· · · ] σ . For instance, we have where where We also have where we identify the permutations σ ∈ S n−1 with the permutations σ ∈ S n such that σ(n) = n. For any operator O on (C N ) ⊗n and any σ, τ ∈ S n with τ 2 = 1, we also have 16 Using this notation (and at intermediate steps the convention P n,n = I), one can see 16 Multiplying P τ from the right of the expansion OP τ = σ P σ OP τ σ , we obtain O = σ P σ OP τ σ P τ = σ P σ•τ P τ OP τ σ P τ = σ P σ P τ OP τ σ•τ P τ , where we used the fact that P 2 τ = 1, P σ P τ = P σ•τ and τ 2 = 1. Taking note on the fact that P τ OP τ σ•τ P τ does not effectively contain the permutation operator, and comparing the above equation with an expansion O = σ P σ [O] σ , we obtain (C33). (C34) follows from a similar argument. ⊗ I) only contains Pσ wherẽ σ(n) = n. Hence, one can only get P k,n Pσ = P σ if k = σ(n). The next line is an elementary manipulation with the permutation operators (using the fact that in I−P n,σ(n) the term I contributes only if σ(n) = n), and the last equality comes out by the same argument as the first one.
For simplicity, the sequence of equalities (C35) will be written as −ζ m z m n−1 , where 17 F (σ) for the operator F (σ) of σ ∈ S n by definition equals to σ P σ [F (σ)] σ .
As a first example of using the above notations, let us note that for any i < n, can also be written as P σ P i,n P −1 σ O 1 O 2 P i,n , which is also equal to P σ P i,n P −1 σ P i,n O 1 O 2 (due to (C12)). A consequence is that, since we have already proven that in O 1 O 2 , the coefficient of (h (n) ) 2 vanishes, we deduce that it also vanishes in P σ P i,n P −1 σ P i,n O 1 O 2 . It means that in O 1 O 2 , the coefficient of (h (i) ) 2 vanishes just like the coefficient of (h (n) ) 2 .
An important remark about the equalities (C35) is that the left-hand side is the This expression allows to conveniently find that for instance On the other hand, the right-hand side −ζ m z m (O [n−1] m I (σ(n)) ) k =σ(n) P k,σ(n) is convenient in order to express (O m P i,j ) (h (n) ) 1 : for instance one sees that where (C33) is used. More generally we get (O m P i,j ) (h (n) ) 1 = −ζ m z m (O [n−1] m I (σ(n)) ) k =σ(n) P k,σ(n) P i,j , for arbitrary i, j. 2 . Using the arguments above (cf. eqs. (C12), (C25), (C35), and the last paragraph above), we see that this term is equal to where we used the notations defined above 18 . In what follows we will show, with the necessary details, how this term vanishes (when S is applied) for permutations σ such that σ(n) = n. We will also discuss very shortly how the argument can be transposed to the simpler case of the permutations such that σ(n) = n.
In order to prove that it vanishes (after the functions S is applied), we can use the recurrence hypothesis S Q w(ζ) . If we take the derivative (i.e. apply d⊗ from the left) of this relation and use the definition (C13) to rewrite Q(z m , ζ m ) in terms of O m , we get Hence, we also have S[f (A)] = 0 for any linear transformation f which commutes with S. In particular we have 19 Assuming that σ(n) = n, this identity reads where I (σ(n)) ) k =n,σ(n) P k ,n P n,σ(n) , ⊗ I), The identity (C40) is just an explicit reformulation of (C39), as discussed in more details in the very end of the present appendix. I (σ(n)) ) P n,σ(n) k =n,σ(n) P k,σ(n) σ = P n,σ(n) k =n,σ(n) P k,σ(n) (O where the first equality comes from the fact that we take the coefficient [. . . ] σ , which implies that in ( k =σ(n) P k,σ(n) ) 2 = k,k =σ(n) P τ (k ,σ(n)) •τ (k,σ(n)) , we should only keep the terms where (τ (k ,σ(n)) • τ (k,σ(n)) )(n) = σ(n), i.e. the terms where k = n and k = n. From (C42), we see that S [B 1 ] σ and S [B 8 ] σ sum up to which is zero by the same argument. We can also express the sum S [B 4 + B 6 ] σ as For simplicity, the argument was here written for permutations σ such that σ(n) = n.
On the other hand, if σ(n) = n then one obtains the same outcome (C46), by writing the condition S n−1 k=1 P k,n A + A k =σ(n) P k,σ(n) σ (h (n) ) 2 = 0, and following the same steps as above, the main difference being that B 3 and B 5 are then absent from (C40) 20 .
By summing up over all permutations, we hence obtain that 2 ) the term of degree one in h (n) (cf. (C37)) does vanish as we wished to prove.
To conclude the proof, we have seen that in S(O 2 ), the terms of degree 2 in any h (i) vanish, as well as the terms of degree 1 in h (n) . Multiplying by permutations, we deduce that the terms of degree 1 in any h (i) vanish, and we have also shown that the term of degree 0 in all the h (i) 's vanishes. Hence we proved (by recurrence) that S(O 2 ] σ ) (h (n) ) 2 . In the case when σ(n) = n, there is no term k = σ(n) in this sum.
More details for equation (C40). In this paragraph we will show that (C40) is equivalent to (C39). First, we can see that the first two terms in this expression are where we used (C25). Next, we can see that [B 3 2 ]σ, whereσ = τ (n,σ(n)) • σ obeysσ(n) = n, and we decompose it as in (C32) then we see that C only involves [O 2 ] σ where σ (n) = n. But from (C31) one sees that for such permutations, one has [O where the first equality relies on the condition σ(n) = n and the fourth equality rewrites (O The fifth equality in (C48) uses the same argument as in the derivation of (C35), namely that when we pick the coefficient [. . . ] σ all terms withk = σ(n) vanish. Then we find 4 k=1 [B k ] σ = ([ n−1 k=1 P k,n A] σ ) (h (n) ) 2 . In order to conclude the proof that (C40) is an explicit rewriting of (C39), we will show that we also have 8 k=5 [B k ] σ = ([P n,σ(n) A P n,σ(n) k =σ(n) P k,σ(n) ] σ ) (h (n) ) 2 , which will require to use a slightly different notation for the degrees in the variables h (i) . So far, we have been using the notation X (h (k) ) p to denote the coefficient in front of the term which has degree p in h (k) in X . In this definition, the ambiguity which arises in the definition of the degree in each individual h (i) (due to relations like P i,j h (i) = h (j) P i,j ) was fixed by conventionally moving all permutation operators to the left of the operators h (i) . We will now also use the notation X (h (k) ) p to denote the coefficient in front of the term which has degree p in h (k) in X when the ambiguity in the degree is fixed by moving all permutation operators to the right. For instance, when i = j, we have Let us now show that 8 k=5 [B k ] σ = ([P n,σ(n) A P n,σ(n) k =σ(n) P k,σ(n) ] σ ) (h (n) ) 2 , by splitting the right-hand side into pieces. The first piece is ([P n,σ(n) O ⊗((1−hz 2 )(1− hζ 2 ))) P n,σ(n) k =σ(n) P k,σ(n) ] σ ) (h (n) ) 2 , which splits into the term k = n on the one hand, and the other terms in the other hand. The term where k = n is given by: by the same argument as in the derivation of B 3 . Next the terms with k = n are given by where the equality (C57) arises because only the terms withk = n contribute to (C56). Finally, the second piece in ([P n,σ(n) A P n,σ(n) k =σ(n) P k,σ(n) ] σ ) (h (n) ) 2 is given by ([P n,σ(n) 2 P n,σ(n) k =σ(n) P k,σ(n) ] σ ) (h (n) ) 2 , which is given by where the equality (C61) uses the fact that [O ((1 − hz 1 )(1 − hζ 1 )) (σ(n)) ]σ vanishes unlessσ(σ(n)) = σ(n).

D The bilinear identity and matrix integrals
In this appendix we give an independent proof of the fact that the master T -operator satisfies the bilinear identity for the KP hierarchy. We use the technique of matrix derivatives (see section 4.2.3) and matrix integrals (see [41]). As it follows from that section, it is enough to prove the bilinear identity in the form (4.25) with an auxiliary matrix A. Here we assume that the matrix A is positively defined and Hermitian 21 .

A complex matrix integral
Let A be a positively defined (Hermitian) matrix and B, C be arbitrary N × N matrices. We need the following Gaussian matrix integral over N × N complex matrices W : with the flat measure Proof of the bilinear identity in the special case x = x 1 = . . . = x n Let us start with the special case of the master T -operator at x = x 1 = . . . = x n given by (4.17). The matrix integral representation (D1) suggests to consider the modified generating function Φ G (t, t ; A) = (det A) N Φ G (t, t ; A) 21 The type of the matrix A is actually not very important for us. We take it to be Hermitian for simplicity reasons. For our purpose it is sufficient to assume that this matrix has enough number (namely, N 2 ), of independent entries, so that one could take derivatives with respect to them in (4.21).
instead of (4.20). Obviously, this does not affect the bilinear identity (4.25) and so it is enough to prove it for Φ G (t, t ; A). The combination of (4.21) and (D1) gives: Here we have used invariance of the complex matrix integral with respect to the shifts of the integration variables W → W + B, W † → W † + C with arbitrary complex matrices B and C. (This invariance clearly holds at N = 1 while the N > 1 case is reduced to a multiple integral of the same type.) Remark. Note that in general the matrix integrals in the last two lines of (D3) diverge, so they make sense only as formal series in t k and t k , each coefficient of this series being well-defined. This is precisely the meaning that we need from the generating function, so one should not worry about convergence of the integrals of this type ((D3) and below) since their integrands have to be understood as power series in t k and t k .
The commutativity of the master T -operators implies that Φ G (t, t ; A) = Φ G (t , t; A). Therefore, we have two different expressions for the same generating function: Then, using (D2), we obtain another integral representation of the same generating function: A complex matrix W can be decomposed as with unitary U , diagonal w and strictly upper-triangular R. The elements of the matrix w are eigenvalues of W : w = diag (w 1 , . . . , w N ). The flat measure on the space of complex matrices is where [dU ] is the Haar measure for the unitary group, ∆(z) = i>j (z i − z j ) is the Vandermonde determinant and c N is an N -dependent constant (see, e.g., [41,43]). There is an equivalent but different decomposition: where V is unitary and Q is strictly upper-triangular so that Q † is strictly lowertriangular. For this decomposition the measure is Let us use decomposition (D7) for the integral from the last line of (D3), and decomposition (D9) for the integral (D6): (D12) Let us re-denote the integration variables in the last integral as V → U (a unitary matrix) and Q → R (an upper-triangular matrix). Then the half-sum of (D11) and (D12) gives which is again well-defined as a formal series in t, t . Then, for each term of this series, one should take the residue at z = ∞ in the l.h.s. of the bilinear identity (4.25). This means that if we treat the integral as a formal series, the z-integration should be performed first: This function is regular at w i = w j for all i and j. Thus the bilinear identity (4.25) is equivalent to the following relation: In fact each term in the sum over k is equal to zero. To see this, let us consider only the dependence of the integrand on the corresponding variable w k . It is easy to see that the integrand is antisymmetric with respect to the interchange w k ↔w k (complex conjugation), hence the integral over the complex plane vanishes. Indeed, the second and the third lines are obviously symmetric while in the first line we have where N (w) does not depend on w k . This expression is obviously antisymmetric w.r.t. the conjugation w k ↔w k , so the integral vanishes and the bilinear identity is thus proved.
Example: N = 1. Let us consider the simplest possible example N = 1. Then A ∈ R + is just a real positive number. Then where the h i 's are the elementary Schur polynomials (2.6) and is a series in z −1 . Then the right-hand side of (D13) reduces to ∞ i,j=1 h i (t)h j (t ) ∞ dz e ξ(t−t ,z)P i,j (z) where the third line follows from the second one and the relation (A3). The integral in the last line is anti-symmetric with respect to the change i ↔ k, thus the whole sum vanishes.
As it follows from the arguments given in section 4.2, the special case of the bilinear identity proven above is already enough for the proof of the general case. However, a sketch of the direct proof by means of the matrix integrals is given below for completeness.

Proof of the bilinear identity in the general case
To prove the bilinear identity for arbitrary x−x i let us consider a generalization of (4.21). This generating function depends on n external matrices A 1 , . . . , A n (all of which are assumed to be positively defined) through The master T -operator (4.1) is a coefficient of the expansion of (D22) in front of the term linear in all A α , thus for the proof of (4.7) it is enough to show that ∞ e ξ(t−t ,z) Φ G (t − z −1 , t + z −1 ; A 1 , . . . , A n )dz = 0.
The prove is a generalization of the one given above. First of all, let us simplify the notation and denote x α − x by x α . As above, we also slightly modify the definition of the generating function: Φ G (t, t ; A 1 , . . . , A n ) = (det A 1 . . . A n ) N Φ G (t, t ; A 1 , . . . , A n ).
(D25) Then Φ G (t, t ; A 1 , . . . , A n ) = n α=1 d 2 W α e − tr n α=1 A −1 α WαW † α −Wα ∂ ∂g t −xα −W † α( ∂ ∂h t −xα) e (tk tr h k +t k tr g k ) × e tr( k≥1 (tk W †k (D26) where we have used the change of variables Again, using the commutativity of the master T -operators and symmetry (D2) of the complex matrix integral, we can obtain another matrix integral representation for the same generating function: Now we apply the same argument as in the previous proof to the integral over the matrix W 1 . Namely, we take the sum of the last lines of (D26) and (D27) with different parametrization of complex matrix W 1 ((D7) and (D9) respectively). After taking the residue in z in (D24) the same argument as before goes through.

E More about commutation relations
As stated in the main text, the commutation of the T -operators can be derived from the Yang-Baxter relation, or can as well be obtained from the commutation of T -operators of the spin chain, by taking the limit η → 0. The derivation of this commutation is a bit technical, and in [10] it was first obtained for T -operators corresponding to symmetric representations (using the Yang-Baxter equation), and then for general representations using the Cherednik-Bazhanov-Reshetikhin formula.
In this appendix, we give another proof of this commutation, based on the results of Appendix C.2, where the commutation relation [ Q(z 1 , ζ 1 ) , Q(z 2 , ζ 2 ) ] = 0 (E1) was proved. The present proof uses the polynomial normalization (3.11), but the normalization does not affect the commutation relations, hence the result holds in the normalization (3.10) as well.
Let us deduce from (E1) that the T -operators corresponding to different hook-representations, but to the same value of x, do commute with each other. At the level of generating functions, we should prove that S z,ζ (y n + d n ) . . . (y 1 + d 1 ) w(z 1 ) w(ζ 1 ) · (y n + d n ) . . .
where y i = x − x i and the notation S z,ζ was introduced in Appendix C.2. This relation (E2) is easily proven by recurrence over the number n of spins, if we notice that its l.h.s. is polynomial of degree at most two in each of the variables y i . The term with degree 2 in y n is exactly the relation with one less spins while the term with degree 1 in y n is its derivative. Hence the l.h.s. is independent of y n . Using the relation d i d j = d j d i , we deduce that it is also independent of all the y i 's. But when all the y i 's are set to zero, (E2) reduces to the relation (E1), which was already proven. This proves the relation (E2), by recurrence over the number n of spins.
We have proven that the T -operators corresponding to hook-representations commute with each other. Then the Giambelli formula (4.11) allows us to deduce that all the operators (x − x n + d n ) . . . (x − x 1 + d 1 )f (h) commute with each other, for all functions of h which are arbitrary linear combinations of characters, if they have the same value of x. Noticing that one deduces that the successive derivatives ∂ k x (x − x n + d n ) . . . (x − x 1 + d 1 )f (h) also commute with all operators (x − x n + d n ) . . . (x − x 1 + d 1 )f (h), which allows one to conclude that the commutation also holds for two operators corresponding to different values of x.