Classical tau-function for quantum spin chains

For an arbitrary generalized quantum integrable spin chain we introduce a “master T -operator” which represents a generating function for commuting quantum transfer matrices constructed by means of the fusion procedure in the auxiliary space. We show that the functional relations for the transfer matrices are equivalent to an infinite set of model-independent bilinear equations of the Hirota form for the master T -operator, which allows one to identify it with τ -function of an integrable hierarchy of classical soliton equations. In this paper we consider spin chains with rational GL(N)-invariant R-matrices but the result is independent of a particular functional form of the transfer matrices and directly applies to quantum integrable models with more general (trigonometric and elliptic) R-matrices and to supersymmetric spin chains.


Introduction
The aim of this paper is to make precise the long anticipated connection between transfer matrices of quantum integrable models and τ -functions of classical integrable hierarchies of non-linear partial differential equations. We show that properly defined generating function for commuting Hamiltonians of quantum models is a τ -function in the sense that it obeys exactly the same bilinear identities of the Hirota form [1,2] as the classical τfunction does. Since the operator-valued generating functions commute for all values of the auxiliary parameters (one of which being the familiar spectral parameter), there is no problem with their ordering in non-linear equations. A similarity between quantum transfer matrices and classical τ -functions was first pointed out in [3] (see also [4,5]), where a discrete integrable dynamics in the space of integrals of motion of a quantum integrable model was introduced. This classical dynamics was identified with the discrete 3-term Hirota equation with special boundary conditions. It is equivalent to the so-called T -system which is a system of fusion relations among mutually commuting transfer matrices for bosonic [6,7] or supersymmetric [8,9] quantum integrable spin chains. Solutions to the Hirota equation carry full information about spectral properties of the quantum model. The diagonalization of transfer matrices by means of the nested Bethe ansatz technique was shown to be equivalent to an "undressing" chain JHEP09(2013)064 of Bäcklund transformations for the discrete Hirota equation. Later this approach was extended to supersymmetric integrable models [10,11].
An important further step was recently made in [12], where an operator realization of the Bäcklund flow describing the "undressing" process was suggested and a commuting family of transfer matrices (T -operators on different levels of the nested Bethe ansatz) depending on a number of discrete variables (labels of representation in the auxiliary space) and on the spectral parameter was constructed. These T -operators obey the discrete Hirota equation not only for usual but also for supersymmetric integrable models, as it was directly demonstrated in [13] for GL(M |N )-invariant spin chains.
The question which remained open after these works was how to embed these Toperators depending on discrete variables into an infinite classical integrable hierarchy with continuous time flows compatible with the discrete ones. In the present paper we suggest an answer which is rather simple and natural. Below in the introduction we describe it in a concise form.
In this paper we consider generalized quantum integrable spin chains with rational GL(N )-invariant R-matrix [14] R λ (u) = u I ⊗ I + N i,j=1 e ji ⊗ π λ (e ij ) (1.1) which acts in the tensor product of the N -dimensional vector representation π and an arbitrary finite dimensional irreducible representation π λ of the group GL(N ) labeled by a Young diagram λ. Here e ij are standard generators of the algebra gl(N ), I is the unit matrix and u ∈ C is the spectral parameter. A family of commuting operators (quantum transfer matrices or simply T -operators), shown in figure 1, can be constructed as T λ (u) = tr π λ R λ (u − ξ L ) ⊗ . . . ⊗ R λ (u − ξ 1 )(I ⊗L ⊗ π λ (g)) , (1.2) where ξ i are arbitrary parameters (inhomogeneities at the lattice sites) and g ∈ GL(N ) is called the twist matrix. The tensor product is taken over the quantum space (in the first space of (1.1)). The trace is taken in the auxiliary space where the representation π λ is realized. The T -operators act in the physical Hilbert space of the model H = (C N ) ⊗L which is the tensor product of L local spaces C N of the vector representations. Formally, our setting includes also models with higher representations at the sites because they can be obtained by fusing several vector representations at the sites with properly chosen parameters ξ i . The Yang-Baxter equation for the R-matrix implies that the T -operators with the same g commute for all u and λ and can be diagonalized simultaneously. They obey some functional relations which are given by the Cherednik-Bazhanov-Reshetikhin (CBR) determinant formulas [15,16]: where λ ′ 1 is the height of the first column of the Young diagram λ and ∅ is the empty diagram. These relations mean that the transfer matrices for arbitrary diagrams functionally depend on the transfer matrices corresponding to the 1-row diagrams. Figure 1. T-operator for an integrable, inhomogeneous spin chain with twisted boundary condition.

JHEP09(2013)064
We call such models "spin chains" in a rather broad sense, not implying existence of any local Hamiltonian of the Heisenberg type (integrable local interactions in general do not exist for inhomogeneous spin chains). However, even in the general case of arbitrary parameters ξ i the model still makes sense as a generalized spin chain with non-local interactions. The "spin variables" are vectors from the spaces C N at each site. Instead, one may prefer to keep in mind integrable lattice models of statistical mechanics rather than spin chains as such. In either case the final goal of the theory is the diagonalization of the transfer matrices of the type (1.2) which is usually achieved by the nested Bethe ansatz method in one or another form (see, e.g., [17][18][19]).
The key point of our approach is dealing with all the T -operators simultaneously by introducing their generating function of a special form which we will call the master Toperator. This generating function depends on an infinite number of auxiliary parameters t = {t 1 , t 2 , . . .} which we call times: 1 T (u, t) = λ s λ (t)T λ (u). (1.4) Here s λ (t) are the standard Schur functions (2.13) and the sum goes over all Young diagrams λ including the empty diagram. By construction, the master T -operators commute for different values of the times: [T (u, t), T (u ′ , t ′ )] = 0. (1.5) Our claim is that the master T -operator is the τ -function of the KP hierarchy with the times t 1 , t 2 , . . . . In other words, the transfer matrices T λ (u) are Plücker coordinates for the element of the infinite dimensional Grassmann manifold corresponding to the τ -function T (u, t). This means that T (u, t) obeys the Hirota bilinear equation 1 In a different form such an operator was introduced in [12].

JHEP09(2013)064
with respect to the t-variables with fixed u. Here the standard notation t + [z −1 ] = {t 1 + z −1 , t 2 + 1 2 z −2 , t 3 + 1 3 z −3 , . . . } is used. This functional relation encodes infinitely many partial differential equations for T obtained after expanding the left hand side in inverse powers of z i .
Moreover, when one incorporates the u-dependence of the master T -operator, the KP τ -function extends to the τ -function of the modified KP (MKP) hierarchy with the times t 0 = u, t 1 , t 2 , . . . . This means that T (u, t) obeys another Hirota bilinear equation which encodes infinitely many differential-difference equations of the MKP hierarchy. In fact this statement is very general and is independent of a specific functional form of T λ (u). It also does not depend on particular features of the quantum models in question: they can be lattice or continuous, bosonic or supersymmetric, with rational, trigonometric or elliptic R-matrices, etc. Most part of the paper is devoted to the proof of this statement. Our technical tools are the realization of the KP hierarchy as a dynamical system on an infinite dimensional Grassmann manifold [20] and the free fermionic construction of the τ -function [21,22] developed by the Kyoto school. The main idea of the proof is to identify the generating series (1.4) for the master T -operator with the Schur function expansion of the τ -function τ n (t) of the MKP hierarchy: τ n (t) = λ s λ (t)c λ (n). (1.8) The key fact proven in section 3 with the use of the free fermionic formalism is that the coefficients c λ (n) (Plücker coordinates) obey the determinant relations of the form c λ (n) = They mean that the Plücker coordinates for arbitrary diagrams are expressed through the basic ones, corresponding to the 1-row diagrams. Identifying n with the spectral parameter u and c λ (n) with the T λ (u), one can see that (1.3) coincides with (1.9). So, the statement that the master T -operator is the τ -function appears to be equivalent to the Cherednik-Bazhanov-Reshetikhin (CBR) determinant relations for T λ (u). For spin chains with rational R-matrices, one can say more. In this case, the master T -operator has an explicit realization in terms of the co-derivative introduced in [13]. It was also remarked in [13] that the basic bilinear identity for the gl(M |N )-characters 2 and their co-derivatives, which was the key part of the direct proof of the CBR formulas, might be related to the Hirota equation for a classical integrable hierarchy. Indeed, now we can see that the "master identity" from [12] is a particular case of the bilinear equation (1.7) for the master T -operator. Therefore, the master T -operator contains all Baxter Q-operators

JHEP09(2013)064
and T -operators on all levels of the nested Bethe ansatz. Presumably, this is also true for models with other types of R-matrices. The most important characteristic property of the master T -operators for finite spin chains with rational R-matrices is that they are polynomials in u: (1.10) The class of τ -functions with this property is well studied in the theory of soliton equations and can be characterized [23][24][25] in terms of algebraic geometry of curves with cusp singularities. Furthermore, the dynamics of the roots u j as functions of times t i , is known [26,27] to be given by the rational Ruijsenaars-Schneider model. Sections 2-4 of the paper contain the necessary material on the classical KP and MKP hierarchies and their τ -functions. Although the main result can be proved without using the formalism of free fermions, we found it instructive to present all arguments in the fermionic operator language of the Kyoto school. In section 2 the free fermion operators and the vacuum states are introduced and the vacuum expectation values are defined. In section 3 the τ -function is introduced as the vacuum expectation value of a group-like element and bilinear identities for it are derived. In most cases we use the standard notation and conventions from [21,22]. Section 4 is devoted to a detailed review of rational solutions to the MKP hierarchy. It is precisely the class of solutions relevant to quantum integrable models with rational R-matrices. The core of the paper is section 5, where the master T -operator is introduced and shown to satisfy the classical bilinear identities and Hirota equations.

Free fermions
In this section we recall the formalism of free fermions [21,22] used for construction of τ -functions of integrable hierarchies.

Fermionic operators
Let ψ n , ψ * n , n ∈ Z, be free fermionic operators with usual anticommutation relations [ψ n , ψ m ] + = [ψ * n , ψ * m ] + = 0, [ψ n , ψ * m ] + = δ mn . They generate an infinite dimensional Clifford algebra. We also use their generating series which can be regarded as free fermionic fields in the complex plane of the variable z.
Of particular importance is the operator where t k are arbitrary parameters (called times). It is convenient to denote the set of times by t = {t 1 , t 2 , . . .} and to introduce their generating function

JHEP09(2013)064
It is easy to check that the fields ψ(z), ψ * (z) transform diagonally under the adjoint action of the operator e J + : In terms of the symmetric Schur polynomials h k (t) defined by the corresponding formulas for ψ n , ψ * n can be written as

Dirac vacua and excited states
Next, we introduce a vacuum state |0 which is a "Dirac sea" where all negative mode states are empty and all positive ones are occupied: (For brevity, we call the modes with indices n ≥ 0 as positive). With respect to this vacuum, the operators ψ n with n < 0 and ψ * n with n ≥ 0 are annihilation operators while the operators ψ * n with n < 0 and ψ n with n ≥ 0 are creation operators. The dual vacuum state has the properties 0| ψ * n = 0, n < 0; 0| ψ n = 0, n ≥ 0.
We also need the "shifted" Dirac vacua |n and the dual vacua n| defined as Similarly to the properties of the Dirac vacuum, for the shifted Dirac vacuum we have: We will also use ψ n |n = |n + 1 , ψ * n |n + 1 = |n , (2.8) n + 1| ψ n = n| n| ψ * n = n + 1| .
(2.9) Excited states are obtained by filling some empty states and creating some holes. Let us introduce a convenient basis of states in the fermionic Fock space F. The basis states JHEP09(2013)064 |λ, n are parameterized by Young diagrams λ and numbers n of the Dirac vacua in the following way. Given a Young diagram λ = (λ 1 , . . . , λ ℓ ) with ℓ = ℓ(λ) nonzero rows, let ( α| β) = (α 1 , . . . , α d(λ) |β 1 , . . . , β d(λ) ) be the Frobenius notation for the diagram λ [28]. Here d(λ) is the number of boxes in the main diagonal and where λ ′ is the transposed (reflected about the main diagonal) diagram λ. In other words, α i is the length of the part of the i-th row to the right from the main diagonal and β i is the length of the part of the i-th column under the main diagonal (not counting the diagonal box). Then (2.10) The definition of the vacua implies that e J + |n = |n while the dual coherent states n| e J + are expanded as linear combinations of the basis vectors as follows: where the sum runs over all Young diagrams λ including the empty diagram, and s λ (t) is the Schur polynomial corresponding to the diagram λ. It can be expressed through the functions h k (t) defined in (2.4) (which are Schur polynomials for one-row diagrams) with the help of the Jacobi-Trudi formula [28]: (2.13)

The expectation values and normal ordering
The vacuum expectation value n| . . . |n is a Hermitian linear form on the Clifford algebra defined on bilinear combinations of fermions by the properties n|n = 1, n| ψ i ψ j |n = n| ψ * i ψ * j |n = 0 for all i, j and n| ψ i ψ * j |n = δ ij for j < n, n| ψ i ψ * j |n = 0 for j ≥ n.
One can see that the basis vectors (2.10) are orthogonal with respect to the scalar product induced by the expectation value: Expectation values of general products of fermionic operators are given by the Wick's theorem. Let v i = j v ij ψ j be linear combinations of ψ j 's only and w * i = j w * ij ψ * j be linear combinations of ψ * j 's only, then the standard Wick's theorem states that

JHEP09(2013)064
One may define the normal ordering • • (. . .) • • with respect to the Dirac vacuum |0 : all annihilation operators are moved to the right and all creation operators are moved to the left taking into account their anti-commutativity under mutual permutations. For example, More generally, for any linear combinations f 0 , f 1 , . . . , f m of the fermion operators ψ i , ψ * j we have the recursive formula where f j means that this factor should be omitted.
In fact a normal ordering can be defined with respect to any vacuum state. Another useful normal ordering is the one defined with respect to the bare vacuum |∞ , which is the absolutely empty state. We denote this normal ordering by × × (. . .) × × . It means that all ψ * 's are moved to the left and all ψ's to the right, with taking into account the sign factors appeared each time one operator is permuted with another. For example, In general, for this normal ordering equation (2.14) holds for the vacuum expectation taken with respect to the bare vacuum.

Group-like elements
Neutral bilinear combinations mn b mn ψ * m ψ n of the fermions, with certain conditions on the matrix b = (b mn ), generate an infinite-dimensional Lie algebra [22]. Exponentiating these expressions, one obtains an infinite dimensional group (a version of GL(∞)). Elements of this group can be represented in the form The inverse element can be written in the same way with the matrix (−b mn ). An important example is provided by the operator e J + (2.2). It is straightforward to see that the group elements of the form (2.15) obey a rather special property that the adjoint action of such elements preserves the linear space spanned by the fermion operators ψ n and the dual space spanned by ψ * n . More precisely, we have: where the matrix R = (R nl ) of the induced linear transformation is given by R = e b . Let us note, following the works of the Kyoto school [22,29], that the group elements can be equivalently represented as normal ordered exponents of bilinear forms. For example, G given in (2.15) can be written as

JHEP09(2013)064
with the matrices A, B determined by the matrix b in (2.15) according to the formulas [29] B = e b − I , Here I is the unity matrix and P + is the projector on the positive mode space ((P + ) ik = δ ik for i, k ≥ 0 and 0 otherwise). The product of two group-like elements is also a group-like element: where C = A + B + AB. As a slight generalization of (2.17), one can obtain the formulas where B (n) = (B ik ) i,k≥n is the half-infinite matrix obtained from B by truncation up to the n-th row and n-th column and Note that the r.h.s. of the last formula is the (rs)-minor of the matrix I + B (n) (up to a sign). The general theory of normal ordering and the group elements in the the Clifford algebra was given in [29]. What is important here, is that the normal ordering allows one to represent in the form (2.17) not only the group elements but also some non-invertible elements of the Clifford algebra that satisfy the same commutation relations (2.16) with some matrix R.
We call the elements G of the Clifford algebra, such that the commutation relations (2.16) hold with some (not necessarily invertible) matrix R, the group-like elements. If the matrix R fails to be invertible, so does G, as an element of the Clifford algebra. In this case it can not be represented in the exponential form (2.15). However, in general it still admits a representation as a normal ordered exponent of bilinear forms in the fermionic operators.
Let us give an example. Let Ψ, Φ * be arbitrary linear combinations of the fermion operators ψ n , ψ * n , respectively, and consider the element where γ := ∞| ΨΦ * |∞ and α, β are related as e γβ = 1 + γα. For general values of α (and for γ = 0) this element is invertible and the two representations are equivalent. However, at α = −1/γ the invertibility breaks down and the element G becomes which can not be written in the exponential form (2.15) but can be represented as the normally ordered exponent.

JHEP09(2013)064
From (2.16) it easily follows that any group-like element obeys the commutation relation which is equivalent to the bilinear identity valid for any states |V , |V ′ , U | , U ′ | from the space F and its dual. This is the basic identity for deriving various bilinear equations for the τ -function (see below).

The generalized Wick's theorem
Let v i = j v ij ψ j be linear combinations of ψ j 's only and w * i = j w * ij ψ * j be linear combinations of ψ * j 's only, as before. A more general version of the Wick's theorem given in [21] is where G is any group-like element. Here we imply that n| G |n = 0 for all n but in fact this restriction can be removed by multiplying the both sides by an appropriate power of n| G |n . The similar formula with the exchange v j ↔ w j also holds. This statement can be proved by induction using the bilinear identity (2.22) (see appendix A). The following particular case is often useful in applications. Writing n − m| = n| ψ n−1 . . . ψ n−m with m > 0, we get from (2.23): For our purpose, it is important to note that there exists an alternative determinant representation of the expectation value in the l.h.s. of (2.24), which is another form of the generalized Wick's theorem: Indeed, the first columns of the two matrices in the r.h.s. of (2.24) and (2.25) coincide and one can show that the j-th column of the matrix in (2.24) is equal to the j-th column of the matrix in (2.25) plus a linear combination of the first j − 1 columns (see appendix A). We also note an equivalent form of (2.25): which is obtained from it by conjugation of the states and operators.

JHEP09(2013)064
3 The τ -function and the Baker-Akhiezer functions As it has been established in the works of the Kyoto school, the expectation values of grouplike elements are tau-functions of integrable hierarchies of nonlinear differential equations. In particular, is the tau-function of the KP hierarchy depending on the times t = {t 1 , t 2 , . . .}, which means that it obeys an infinite set of Hirota bilinear equations [1,2]. More generally, introducing a "discrete time" n via one obtains a τ -function which solves the modified KP (MKP) hierarchy [30][31][32][33]. 3 Equations of the MKP hierarchy are differential-difference equations which include shifts of the variable n. At each fixed n, the τ -function (3.2) is a τ -function of the KP hierarchy, so the n-evolution represents an infinite chain of Bäcklund transformations. The meaning of the discrete variable n depends on the class of solutions. For some classes of solutions n can be regarded as a continuous, or even complex, variable. The last point will be of a prime interest for us since in section 5 n will be identified with the spectral parameter u of the quantum transfer matrix of a spin chain.

The bilinear identity
The τ -function obeys a bilinear identity which is a direct consequence of (2.22). Setting |V = |n , |V ′ = |n ′ with n ≥ n ′ , where |n and |n ′ are two shifted Dirac vacua, we have because either ψ k or ψ * k kills the state in each term in the r.h.s. of (2.22). Now, setting U | = n + 1| e J + (t) , U ′ | = n ′ − 1| e J + (t ′ ) , we can write where we have used the commutation relations of the fermion operators with e J + (t) and the bosonization rules

JHEP09(2013)064
established by the Kyoto school. Here and below we use the standard notation In this way we arrive at the bilinear identity for the τ -function valid for all t, t ′ and n ≥ n ′ . It encodes all the PDE's of the MKP hierarchy.
The choice of the integration contour C depends on the type of the time evolution. Formally, the evolution factor e ξ(t−t ′ ,z) z n−n ′ has an essential singularity at ∞, and the contour should then encircle the whole complex plane (only the singularity at ∞ stays outside the contour). This is the usual prescription and if only a finite number of times t k are nonzero, then it is correct. However, in general, when the values of the times are such that the factor e ξ(t−t ′ ,z) z n−n ′ has singularities for finite values of z, then the prescription is as follows: the contour C must encircle all singularities of the function , leaving all the singularities of e ξ(t−t ′ ,z) z n−n ′ outside the contour.
Remark. There is a freedom to multiply the τ -function of the MKP hierarchy by an exponent of any linear form of times with constant coefficients and by an arbitrary function of n: Clearly, this transformation preserves the form of the bilinear identity. By noting that we see that the transformation τ n (t) → C(n)τ n (t) means G → G × × e j d j ψ * j ψ j × × for the group-like element. In the MKP hierarchy, there are no restrictions on the form of C(n). Such restrictions appear in the more general 2D Toda lattice hierarchy, where this function is constrained to be of the form C(n) = Ce C 0 n with constant C, C 0 .

3-term bilinear equations
, we see that the essential singularity at ∞ splits into 3 simple poles at z 1 , z 2 , z 3 : .
Taking the residues, we arrive at the 3-term relation

JHEP09(2013)064
where the last two terms are obtained from the first one by the cyclic permutations of the indices. Another choice is n ′ = n and t ′ .
This leads to a slightly more general equation , we see that the essential singularity at ∞ splits into simple poles at z 1 , z 2 : .
Besides the residues at these points, there is also a contribution from the residue at ∞, so we obtain the 3-term relation It can be formally regarded as a particular case of (3.8) in the limit z 3 → ∞, z 0 → 0. The limit z 3 → ∞ is smooth while the other one requires to assign a meaning to lim The form of the evolution multiplier, z n e ξ(t,z) , suggests to set formally which actually holds if the function τ n (t − [z −1 ]) for arbitrary n can be analytically continued from a neighborhood of ∞ to a neighborhood of the point z = 0 and is regular there.
If it is singular, as in the case of solutions relevant to quantum spin chains, (3.10) should be substituted by a more general prescription (for an example see (3.14) below) which, however, also allows one to regard (3.9) as a particular case of (3.8).

Examples of τ -functions
Here we consider some particular cases of the general fermionic construction of τ -functions which will be important for us in section 5. Let us start from a very simple example. Set Ψ = ψ(p) + αψ(q) and consider the operator • • e βΨψ * (r) • • (with |r| < min (|p|, |q|)). It is a group-like element for any β but at β = β 0 = ( 0| Ψψ * (r) |0 ) −1 it is not invertible and equals β 0 Ψψ * (r). Therefore, τ n (t) = n| e J + (t) Ψψ * (r) |n is a τ -function of the MKP hierarchy. Then we would like to consider the limit r → 0. The limit is singular and requires some care. From we see that in order to get a well-defined limit, one should multiply the τ -function by r n−1 before tending r → 0 (this is just a transformation of the form (3.6)). Then we come to the conclusion that τ n (t) = n| e J + (t) Ψ |n − 1 = n| e J + (t) ψ(p)+αψ(q) |n−1 = p n−1 e ξ(t,p) + αq n−1 e ξ(t,q) is a τ -function. Indeed, it is the τ -function for the 1-soliton solution of the MKP hierarchy. This example can be generalized. Fix a finite set of distinct points p i ∈ C, i ∈ I, and construct a (finite) linear combination of the operators ψ(z) and their z-derivatives at the points p i : Examples of the operators Ψ I are ψ(p), ψ(p) |n is a τ -function for any fixed β. Let us take β = ( 0| Ψ I ψ * (r) |0 ) −1 , then the group-like element is non-invertible but n| e J + (t) Ψ I ψ * (r) |n is still a τ -function. The limit r → 0 can be performed as in the previous example and we conclude that This example can be further generalized by considering products of operators of the form Ψ I . Take N operators of the form Ψ I for different sets I 1 , . . . , I N : Ψ I 1 , . . . , Ψ I N . As before, we can construct the group-like elements Note that if the auxiliary points r j are all distinct from the points p k labeled by the sets I i , r j = p k for all j, k ∈ I 1 ∪ . . . ∪ I N , then these operators commute. Again, we choose β j = ( 0| Ψ I j ψ * (r j ) |0 ) −1 , then each operator becomes the product β j Ψ I j ψ * (r j ) (non-invertible). Thus the above τ -function reduces to a τ -function of the form (β 1 · · · β N ) n| e J + (t) Ψ I 1 · · · Ψ I N ψ * (r N ) · · · ψ * (r 1 ) |n . In order to implement the limit r j → 0, we redefine r j → εr j with ε → 0, then it is not difficult to verify that is the Vandermonde determinant. The first equation says that in order to get a well-defined limit, one should multiply the τ -function by ε (n−1)N − 1 2 N (N −1) before tending ε → 0. The r i -dependent factors (r 1 . . . r N ) −n+1 ∆ N (r i ) in the r.h.s. of (3.12) are irrelevant because they can be eliminated by a transformation of the form (3.6). We thus conclude that is a τ -function.

JHEP09(2013)064
In this case the prescription (3.10) is not directly applicable. The behavior of the function τ (t ± [z −1 0 ]) at z 0 → 0 can be found by using the bosonization formulae (3.4) in the opposite direction and moving the fermion operators ψ(z 0 ) or ψ * (z 0 ) to the very right position. In this way we obtain a more general version of the prescription (3.10): which is equally enough to deduce the bilinear equation (3.9) from (3.8) as a limiting case.

Plücker coordinates
One may expand the τ -function in the Schur polynomials: The coefficients c λ (n) ("Plücker coordinates") can be determined from the formula (3.2) by inserting the complete set of states in between the operators e J + and G and using (2.11): The Schur function expansions of τ -function are also discussed in [34][35][36][37][38]. The coefficients c λ (n) generalize the characters of representations of the linear groups. Remarkably, they admit determinant representations of the Giambelli type as well as of the Jacobi-Trudi type which express c λ (n) for arbitrary λ through c λ (n) for the hook and the one-row diagrams, respectively.

Formulas of the Giambelli type
Applying to (3.16) the Wick's theorem in the form (2.23), we obtain the Giambelli-like formula for c λ (n): JHEP09(2013)064

Formulas of the Jacobi-Trudi type
Alternatively, we may apply to (3.16) the Wick's theorem in the form (2.25). The details are given in appendix B. The result is where c s (n) := c (s−1|0) (n) = n−1| ψ * n+s−1 G |n are the expansion coefficients for one-row diagrams and λ ′ 1 is the height of the first column of the diagram λ. Therefore, the Plücker coordinates of any MKP τ -function satisfy the determinant relations (3.18) of the Jacobi-Trudi type sometimes called quantum Jacobi-Trudi formulas. Note that the transformation c λ (n) → C(n)c λ (n) with arbitrary C(n) preserves the determinant formula (3.18). Clearly, this freedom corresponds to the possibility of multiplying the τ -function by any C(n), see (3.6).
In fact the inverse statement is also true: given any set of quantities c λ (n) that satisfy the determinant relations (3.18), the function λ c λ (n)s λ (t) is a τ -function of the MKP hierarchy. In order to prove this, it is enough to show that for any given set of quantities Let us give a sketch of proof. First of all we note that the element G, if it exists, is not unique because × × e i>k A ik ψ * i ψ k × × |n = |n for any strictly lower triangular matrix A. This means that we can consider only upper triangular matrices. For simplicity we assume that c ∅ (n) = 0 for all n. Then we can (and from now on will) put c 0 (n) = 1 without any loss of generality because this is the matter of normalization due to the freedom to multiply the MKP τ -function by an arbitrary function of n as in (3.6). It is convenient to formally put c k (n) ≡ 0 for k < 0. Consider the infinite matrix C nm = c m−n (n + 1) (m, n ∈ Z). It is an upper triangular matrix with 1's on the main diagonal. One can see that for any matrix C of this form there exists a strictly upper triangular matrix B such that Indeed, according to (2.20),

Example: characters
As an example, let us consider the τ -function corresponding to the coherent states

JHEP09(2013)064
where J −k are given by equation (2.2), x k = 1 k tr g k , g ∈ GL(N ). The τ -function (3.2) is explicitly given by where p i are eigenvalues of the matrix g. This τ -function does not depend on n. Therefore, the coefficients c λ (n) are n-independent and coincide with GL(N )-characters evaluated for the element g ∈ GL(N ): The formulas (3.17) and (3.18) become the well-known determinant formulas for characters [28] (Giambelli and Jacobi-Trudi respectively).
Remark. This simple example may be somewhat misleading because in this case the τ -function (3.20) is simultaneously the τ -function with respect to the "negative" times t −k = x k . For more complicated solutions constructed below from the group element g ∈ GL(N ) with x k = 1 k tr g k this is not true in general.

Baker-Akhiezer functions and wave operators
The Baker-Akhiezer (BA) function and its adjoint are given by the Sato formulas In terms of the BA function and its adjoint, the bilinear identity (3.5) acquires the form An important role in the theory [39] is played by the wave operator (or dressing operator ) W (n, t) directly related to the BA functions. It is a pseudo-difference operator of the form and its formal inverse

JHEP09(2013)064
The Baker-Akhiezer functions are obtained by applying the wave operators to the bare exponent z n e ξ(t,z) : . .} is used. These formulas allow us to express the wave operator and its inverse in terms of the τ -function as follows (cf. [40]): The (pseudo-difference) Lax operator of the MKP hierarchy is obtained by dressing the shift operator e ∂n with the wave operators: The Lax operator obeys the Lax equations where (L k ) ≥0 means the (finite) part of the series containing the shift operators e j∂n with j ≥ 0.

Rational solutions of the KP and MKP hierarchy
In this section we study rational solutions of the KP and MKP hierarchies. For these solutions, the τ -function is a polynomial function of the times n and t i (possibly multiplied by an exponential function of a linear combination of the times). Stressing the fact that n for solutions of this class is not necessarily integer but can be any complex number u ∈ C, we change the notation n → u.

The construction of rational solutions
The general construction of rational solutions to the KP hierarchy was given by Krichever in [23], see also [24,25]. This construction can be directly extended to the MKP hierarchy.
In this section we show how it works starting from the bilinear identity (3.5). The solutions are obtained in the Casorati determinant form which is a discrete analog of the Wronskian determinant (see, e.g. [41]).

Baker-Akhiezer functions and bilinear identity
First of all, let us specify the bilinear identity for the case at hand. Since in the solutions we are going to consider u is not necessarily integer, 0 and ∞ are in general branch points of the BA function. In this case the appropriate form of the bilinear identity is where [0, ∞] is an arbitrary cut between 0 and ∞ and the contour C [0,∞] is such that it encircles the cut but does not enclose any other singularity of the BA functions.
According to the Krichever's theory of rational and more general algebro-geometric solutions, they can be characterized and explicitly constructed by fixing certain analytic properties of the BA function on a Riemann surface. For rational solutions, the Riemann surface is the Riemann sphere (compactified complex plane), which represents a genus zero algebraic curve with singularities. Correspondingly, the BA function is, in this case, a rational function on the complex plane multiplied by power-like and exponential factors which give the required asymptotics (3.27) (the essential singularity at infinity). Let us assume the following ansatz for the BA function: where the coefficients w i depend on the times u, t j . In fact this is just the truncated series (3.27). The multiplicity N of the pole at z = 0 of the function z −u e −ξ(t,z) ψ u (t, z) is a discrete parameter characterizing the class of solutions to be constructed. Given the BA function of the form (4.2), an important dynamical information is contained in the adjoint BA function. In general, it may have poles of arbitrary order at some points p i ∈ C (we assume that p i = 0). As we shall see soon, the bilinear identity is consistent when the number of these points does not exceed N , and, in case of general position their number is just equal to N . Let us adopt this and fix N distinct points p i ∈ C, p i = 0, and also assume that the adjoint BA function has poles of orders M i + 1 at p i . Therefore, its general form is whereã im (u, t) is set to be 0 for m ≥ M i + 1. As we shall see in section 5, the τ -functions of our principal interest (master T -operators for spin chains with rational R-matrices) are, by construction, polynomials in u. This means that the coefficientsã im (u, t) are rational functions of u, so we assume this in the present section.

Krichever's conditions from the bilinear identity
Since the poles at p i are the only singularities of the product ψ u (t, z)ψ * u ′ (t ′ , z) outside the contour C [0,∞] , the bilinear identity (3.24) then implies that for all u, t j , u ′ , t ′ j . Let us rewrite this identically in the form where the functions R i are From the polynomiality of the τ -function it is clear that they are rational functions of u ′ . If so, the only way to satisfy (4.4) for all u ′ is to put each term equal to zero separately (because it is a sum of different exponential functions multiplied by rational functions which can be identically zero if and only if all rational functions vanish identically). Therefore, the bilinear identity is equivalent to N conditions 5) or, more explicitly, where a im ≡ k iãim (0, 0) with arbitrary non-zero k i are parameters of the solution in the Krichever's approach. We call equations (4.7) the Krichever's conditions [23].
Remark. In the Krichever's approach, the parameters a im can be taken arbitrary while in our caseã im (0, 0) appear to be constrained by N conditions of the form (4.23), as we shall see below. That is why we prefer to work with the set of independent parameters a im .
The initial values of the adjoint BA function coefficientsã im (0, 0) are not independent and differ from the a im by the properly chosen N factors k i . The set of conditions (4.7) yields a system of N linear equations for N coefficients w k which demonstrates the consistency of the bilinear identity for BA functions of the assumed form and allows one to fix them. Then the τ -function associated with the ψfunction according to (3.22) solves the MKP hierarchy. The points p i and the entries of the matrix a im are free parameters of the solution. The coefficients w k appear to be rational functions of their arguments u, t i while the τ -function is a polynomial (multiplied by the exponential function of a linear combination of times). From the algebro-geometric point of view this solution is associated with a highly singular algebraic curve which is the Riemann sphere with cusp singularities at the points p i .

Solving the linear system
It is easy to see that conditions (4.7) are equivalent to the system of linear equations These functions are polynomials in u, t i multiplied by the exponential factor p u i e ξ(t,p i ) . In the case when all t i are equal to 0, they admit a more explicit representation: The system (4.8) can be solved by applying the Cramer's rule. This results in the following determinant representation for the BA function: .
From the definition (4.9) we have (4.12) Using this, it is straightforward to verify that equation (4.11) agrees with the general relation (3.22), with the τ -function being given by the determinant in the denominator: .
p u i e ξ(t,p i ) . It follows from (4.11) that the last coefficient in (4.2), w N , is given by (4.14)

JHEP09(2013)064
We also note the property from which one can see that the first coefficient in (4.2), w 1 , is given by (4.16)

The adjoint BA function
Similarly to (4.12), we have: Note that this function regarded as a function of z has a pole of order M k + 1 at z = p k . The principal term is We also see from (4.17) that the function A k (u, t + [z −1 ]) has a simple zero at z = 0.
In order to obtain a more detailed information about singularities of this function, let us add and subtract the term z u e ξ(t,z) in the numerator in (4.17) and separate the nonsingular part as z → p k from the singular one: (4.19) The first sum gives the multiple pole structure at the point p k while the second term is obviously regular at p k and has possible (essential) singularities and branching only at 0 and ∞. Note that A k (u, t + [z −1 ]) is, by construction, a rational function of z with simple zero at z = 0. This is obvious from (4.17) but not from (4.19). On the other hand, the representation (4.19) explicitly shows the multi-pole structure of this function which is rather obscure in (4.17).
Rewriting (4.12) in the form it is straightforward to check that (4.20) Expanding (4.17) in powers of z, we get

JHEP09(2013)064
and so the expansion of τ u (t + [z −1 ]) around ∞ reads . (4.21) We thus see that the adjoint BA function has the determinant representation .

(4.22)
It has multiple poles at the singular points p i of the algebraic curve which is the Riemann sphere with complicated cusp-like singularities. The principal parts at the poles are easily extracted from the determinant representation and equation (4.19).
) has a simple zero at z = 0, it is clear from (4.13) that the function τ u (t + [z −1 ]) and thus the function z u e ξ(t,z) ψ * u (t, z) have zero of order N at z = 0, whence the coefficientsã im (u, t) in (4.3) appear to be constrained by N relations 4.1.5 The multi-pole structure of the adjoint BA function Let us introduce the notation Then, expanding the determinant in the numerator of (4.22) in the first column, we obtain: Using (4.19) we can extract the poles:

JHEP09(2013)064
Therefore, for all n ≥ 0, or, in terms of the τ -function, Note that the expression in the r.h.s. has the same structure as (4.13) and, therefore, is a polynomial τ -function. This fact will be used below for the Bäcklund transformations. We also note the factorization formula which follows from (4.25). At u ′ = u = t ′ j = t j = 0 this formula gives the relation betweeñ a k0 ≡ã k0 (0, 0) and a k0 :ã

Undressing Bäcklund transformations for the rational solutions
As it will be clear in section 5, the main relations of the Bethe ansatz method are naturally built into the construction of rational solutions to the MKP hierarchy. In particular, the functionsĀ k (u, t) and A k (u, t) will be identified, up to some irrelevant factors, with the (eigenvalues of) T -operators on the first and the last levels of nesting in the nested Bethe ansatz scheme. The parameters p i will be identified with eigenvalues of the twist matrix g. We will also see that all these formulas can be understood in the operatorial sense (in the quantum space of the spin chain) since they involve only commuting T -operators. The nested Bethe ansatz scheme itself appears to be equivalent to a chain of some special Bäcklund transformations of the initial rational MKP solution with Krichever's data p i , a im that "undress" it to the trivial solution by reducing the number of singular points in succession. Here we present the idea solely in terms of the MKP hierarchy leaving the precise identification with objects from quantum integrable spin chains for the next section.
Adding/removing a point p i to/from the Krichever's data of a rational solution is a Bäcklund transformation. It sends a polynomial τ -function to another one. The key equation that allows one to implement such a transformation is (4.27). Let us start with the n = 0 case. Specifically, consider the function (4.30) According to the n = 0 case of (4.27), it is equal to i.e., to the minor M i,1 of the matrix A i (u + 1 − j) multiplied by a i0 . Therefore, it is a τ -function, i.e., it satisfies the same 3-term Hirota equation as τ u (t) does and τ → τ [i] is indeed a Bäcklund transformation. Below we assume that a i0 = 0 for all i = 1, . . . , N , then the transformation is non-trivial.
Having at hand τ . (4.33) It has multiple poles at the points p 1 , . . . , p N except the point p i . This process can be continued by taking the residue of z −u−1 e −ξ(t,z) τ which leads to In the first line d n = i 1 + . . . + i n + 1 2 n(n + 1) and

JHEP09(2013)064
is the Vandermonde determinant. In particular, on the second to last level we have and on the last level τ [12..
Note that a i0 = A i (0, 0). By taking residues of the Hirota equations in the form (3.7), (3.9), one can obtain bilinear relations for the τ -functions on the neighboring levels of this chain. In particular, we have: We can also rewrite the determinant (4.36) in terms of (4.38):   In a similar way, one may define the "undressing chain" of such transformations as in (4.35) and (4.36): which leads to The construction of the undressing based on the highest poles is somewhat more complicated because it explicitly uses the multiplicities of the highest poles but the advantage is that it always gives the non-vanishing result (because the coefficients in front of the highest poles cannot vanish by their definition).

Fermionic realization of rational solutions
Here we suggest a realization of the rational solutions as vacuum expectation values of some fermionic operators. As before, we fix N points p i and the coefficients a im . Let us introduce the following fermionic operators: The operators Ψ i (p i ) are particular cases of Ψ I i (3.11) when each set I i contains just one element N + 1 − i. Therefore, we already know that the matrix element is a τ -function.
On the other hand, we can calculate it using the Wick's theorem in the form (2.26). We get: which coincides with the τ -function (4.13) after the change n → u. Let us also note the formulas which directly follow from the bosonization rules (3.4). The simplest example is M i = 0 when Ψ i (p i ) = ψ(p i ) (we take a i0 = 1) and

JHEP09(2013)064
where m is a non-negative integer (0 ≤ m ≤ N ). The τ -functions corresponding to (4.43) have the same form except for the prefactor (· · · ). These τ -functions reduce to (4.46) for m = 0. The determinant formulas (4.36) for these τ -functions follow from the Wick's theorem (2.26) in the same way as m = 0 case (4.47). The BA functions for (4.50) are defined by (4.51) The BA functions for (4.43) are defined exactly in the same way. In addition, both of them coincide under the condition that the prefactor of the τ -functions are not zero. By using a formula similar to (4.48), we obtain (4.52) Then the Krichever's conditions for the Bäcklund transformations follow straightforwardly from (4.52) and the relation (Ψ k (p k )) 2 = 0.
5 The τ -function in quantum integrable models

T -operators
We consider generalized quantum integrable spin chains with rational GL(N )-invariant R-matrix. The R-matrix acting in the tensor product of the N -dimensional vector representation of GL(N ) and an arbitrary finite dimensional representation π λ (labeled by a Young diagram λ) has the form where u ∈ C is the spectral parameter 4 and e ij are standard generators of the algebra gl(N ) with matrix elements (e ij ) αβ = δ iα δ jβ . In the second term, the vector representation (the one for the Young diagram with one box) π (e ji ) is abbreviated to e ji . A large family of commuting operators (quantum transfer matrices or simply Toperators) can be constructed as where ξ i are arbitrary parameters (inhomogeneities at the lattice sites) and g ∈ GL(N ) is called the twist matrix. For technical simplicity, we assume below that the twist matrix

JHEP09(2013)064
has distinct, non-zero eigenvalues. The trace is taken in the auxiliary space where the representation π λ acts. The T -operators act in the physical Hilbert space of the model H = (C N ) ⊗L . The Yang-Baxter equation implies that the T -operators with the same g commute for all u and λ and can be diagonalized simultaneously. It follows then from (5.1) and (5.2) that all their eigenvalues are polynomials in u of degree ≤ L.
It is natural to define the T -operator for the trivial representation (empty Young diagram) as where the multiplication by the identity operator in the r.h.s. is implicit. We introduce a special notation for the polynomial in the r.h.s. This polynomial characterizes the model through its zeros ξ i . Another normalization, sometimes more convenient, is to set In this normalization all T -operators are rational functions which may have poles only at the points ξ i and T ∅ (u) = 1.

Determinant formulas for T -operators
The T -operators for different representations in the auxiliary space are known to obey an infinite number of functional relations which follow from the fusion procedure and, eventually, from the Yang-Baxter equation. First, there are determinant formulas which express T λ (u) for arbitrary λ through T -operators for symmetric representations T s (u) := T (s) (u) corresponding to one-row diagrams (s) [15,16]. They are called Cherednik-Bazhanov-Reshetikhin (CBR) formulas or quantum Jacobi-Trudi identities. In our normalization they have the form A direct proof "from the first principles" for the models with the rational GL(N |M )invariant R-matrix was given in [13]. The proof uses the operation of co-derivative with respect to the twist matrix g (see below). When the dependence on the spectral parameter disappears (say, in the limit u → ∞) one obtains the standard Jacobi-Trudi formulas (2.13) for characters. We note in passing that the determinant formula for the generalized Schur functions (constructed in terms of the Gauss decomposition of a matrix (resp. orthogonal polynomials)) obtained in [42] (resp. [43]) has a structure similar to (5.6).

JHEP09(2013)064
Second, there are determinant formulas of the Giambelli type which express T λ (u) for arbitrary λ through the T -operators T l,k (u) := T (l+1,1 k ) (u) corresponding to the hook diagrams (l + 1, 1 k ) [44]: Note that there are no shifts of the spectral parameter here. In appendix C we prove that the determinant representation of the Giambelli type (5.7) follows from the quantum Jacobi-Trudi identities (5.6). In the normalization (5.5) the pre-factors cancel and the determinant formulas look simpler: Below we will argue that these formulas allow one to establish a close connection with classical MKP integrable hierarchy and to recover a hidden free fermionic structure in the auxiliary space.
We emphasize that the notion of the master T -operator is independent of the specific form of the R-matrix. 6 The key point is that T λ (u) in (5.9) is the transfer matrix obtained by the standard fusion procedure in the auxiliary space (the horizontal leg) labeled by the Young diagram λ. The master T -operator is thus well-defined for spin chains with trigonometric and elliptic R-matrices as well and the statement that it is a τ -function still holds. The statement is also applicable to the T -operator corresponding to any higher representation in the quantum space of the model. It can be also applicable to the qcharacters of representations of quantum affine algebras [47].
R be the universal R-matrix associated with Uq(ĝl(N )). R is considered to be an element of B+ ⊗ B−, where B± are Borel subalgebras. Let π λ (u) be an evaluation representation of Uq(ĝl(N )) based on a tensor representation of Uq(gl(N )) labeled by the Young diagram λ. This representation depends on the spectral parameter u since the evaluation map does. We can define the universal T -operator (cf. [45]) as T λ (u) = (Tr⊗id)(π λ (u)⊗id)((g⊗1)R), where g is a twist element (a function of the Cartan elements). Then our universal master T-operator is obtained by substituting this into (5.9). Now (5.9) is an element of B−. One can obtain a master T -operator on a particular quantum space V if we consider a representation of B−. For example, if we consider V = π (1) (ξ1) ⊗ π (1) (ξ2) ⊗ · · · ⊗ π (1) (ξL), where {π (1) (ξj)} L j=1 are fundamental representations we obtain the master T -operator for a lattice model in section 5.3.3 (as a limit q → 1). If we consider a space of the vertex operators of CFT for V , we will obtain a master T -operator for CFT. It is tempting to rewrite the universal master T -operator as a kind of determinant over some function of (π (1) (u) ⊗ id)(R), where π (1) (u) is the fundamental representation. This should be a generalization of a generating function of transfer matrices (cf. [46]).

Hirota equations for the master T -operator
As any τ -function, the master T -operator satisfies the bilinear identity for all u, u ′ , t, t ′ . Here the integrand should be interpreted as a formal power series on z and the contour C should encircle z = ∞. By standard manipulations, one can derive from (5.16) the infinite KP and MKP hierarchies of differential equations. Making the shift of variables t → t + a, t ′ → t − a in (5.16), where a = (a 1 , a 2 , . . . ) and setting u − u ′ ∈ Z ≥0 , one obtains, in the usual way, the following Hirota-type bilinear differential equation for the MKP hierarchy: where D t k is the Hirota derivative with respect to t k . For example, the first equation of the MKP hierarchy reads The bilinear relations (3.7), (3.9) for T (u, t) can be written as

The master T -operator in models with rational R-matrices
For models with rational R-matrices, the T -operators T λ (u) (5.2) admit a more explicit realization in terms of the co-derivative operators acting on functions of the twist matrix g ∈ GL(N ) [13]: where χ λ (g) = tr π λ (g) is the character of g in the representation π λ . All necessary information about the co-derivativeD is given in appendix D. Set As immediately follows from its definition, the master T -operator (5.9) can be written in terms of the co-derivative as follows Note that each term in the direct product introduces a spin into the underlying quantum spin chain. Alternatively, using the fact that (det g) −uD (det g) u =D + u, we can also write it as hiding the spectral parameter into the trivial τ -function (with no spins) In this form it is clear that the spectral parameter is naturally combined with t i into a hierarchy of times u = t 0 , t 1 , t 2 , . . .: The explicit form of the simplest T -operators for a chain of length L = 1, 2 is: JHEP09(2013)064 where P is the permutation operator (see appendix D for definitions). For clarity, the unit matrix I is everywhere written explicitly.
Using the method of [12], it is possible to show that the master T -operator contains the whole family of Baxter Q-operators as well as the T -operators that emerge on intermediate levels of the nested Bethe ansatz. This will be explained in section 5.5.
It follows from (5.26) that . (5.31) Using these formulas, one can explicitly verify that Comparing with the prescription (3.14), we see that there is an extra factor (det g) ±1 which is due to a slightly different definition of the τ -function: in our case it is a pure polynomial in u while in (3.14) it is a polynomial multiplied by an exponential function. DenotingD L (u) := L i=1 (u−ξ i +D) and w(z) = (det(I − zg)) −1 for brevity, we can rewrite (5.19), (5.20) in the form In the second equation, one can recognize the "master identity" from [12] (to identify them, the redefinition z i → 1/z i is required). The relation (5.32) implies that equation (5.34), which is the Hirota equation for MKP hierarchy, follows from (5.33), which is the Hirota equation for the KP hierarchy, in the limit z 3 → 0. Note that at t i = 0 (i.e., F g = 1) and u = ξ 1 = · · · = ξ L = 0, the third term in equation (5.34) vanishes sincê D L · 1 = 0 and the equation becomes the "basic identity" (5.9) from [13] where we put z = z −1 1 , z ′ = z −1 2 .

JHEP09(2013)064
When z → 0, (5.35) reduces to the following identity [13] linear in w, where we used w(z) = 1 + z tr g + O(z 2 ) and Here P 1...L = P 12 P 23 · · · P L−1,L is the full cyclic permutation in the quantum space stemming from (D.3) and from the repeated use of (D.7) (see appendix D). The subscript in g 1 means that it acts, as a matrix, in the space number 1 of the spin chain. Note that the order of operators in (5.35) is irrelevant due to the commutation of transfer-matrices with arbitrary auxiliary representations. Similarly, for z ′ → 0, in the first order in z ′ we obtain: It is a direct consequence of (5.36) if we recall that the two factors in (5.38) commute: both are particular cases of the same family of commuting transfer-matrices. Eqs. (5.36) and (5.38) were proven in a combinatorial way in [13]. Multiplying the left hand sides and the right hand sides of (5.36) and (5.38) we restore again the relation (5.35) which served there as the main tool for the direct proof of the CBR formulas. As is explained above, equation (5.38) was proven in [13] in a combinatorial way. Let us now show how the case F g = 1 of (5.34), i.e., can be deduced from it. First we notice that when L ≥ 1, equation (5.35) is the particular case of (5.39) when ξ i = u for all i = 1, . . . L. Next we notice that although (5.35) only holds if L ≥ 1, the third term in (5.39) is such that (5.39) also holds (trivially) when L = 0.
Since the left hand side of (5.39) is polynomial in the variables u i := u − ξ i , this allows for a recurrence procedure which proves that the term with an arbitrary degree in each u i = (u − ξ i ) in (5.39) is equal to zero. See appendix E for a detailed proof on this.

Proof of the determinant representation and the Hirota equation for the master T -operator
A chain of arguments in [12,13] allows one to give a self-contained proof "from the first principles" of the central claim of this paper: the master T -operator for rational spin chains is a tau-function of the MKP hierarchy. Here we give a direct proof of this statement using the representation of the master T -operator through the co-derivatives. The master identity in the form (5.20) or (5.34) directly follows from the determinant formula (5.22) by applying the Jacobi identity. So in order to prove the master identity it JHEP09(2013)064 is enough to prove (5.22). In fact, for this proof we can assume, without loss of generality, that t = 0. Indeed, the product 7 can approximate F g (t) for any set of t with any needed accuracy, if K is sufficiently large, by an appropriate choice of the parameters z 1 , z 2 , . . . , z K . In other words, we change the variables from t to z 1 , z 2 , . . . , z K (which is the Miwa change of variables): or, for a subset {i 1 , i 2 , . . . , i n } ⊂ {1, 2, . . . , K}, Then, in order to prove (5.34) for arbitrary t, it is sufficient to prove (5.34) for arbitrary t = 0 + K k=1 [z k ], which is equivalent to proving (5.22) when t = 0. Therefore, we only have to prove the identity , (5.43) This formula expresses the general T {i 1 ,i 2 ,...,in} (u) through T {i} (u) with only one index (the generating operator for transfer matrices in symmetric representations). The proof goes along the same lines as the proof of the CBR formula in [13]. It consists of two steps: 1. We should prove that T {i 1 ,i 2 ,...,in} (u) is a polynomial operator of the power L, as it is clear from the definition; 2. We should prove that the coefficients of polynomials in the r.h.s. and the l.h.s. of (5.43) are equal.
In order to prove the first statement, it is enough to show that the determinant in the numerator vanishes at all zeros of the denominator, namely at u = ξ l +j, l = 1, . . . , L, j = 1, . . . , n − 1. For that, let us compare two consecutive columns, the j'th and the (j + 1)'th, of the matrix in the numerator. If the 2 × 2 minor determinant

JHEP09(2013)064
built out of any two arbitrary rows k and k ′ of these two columns are zero, then the whole determinant vanishes. 8 Indeed, the corresponding operatorial identity appeared in [13] in relation to the proof of the CBR formula. This identity can be deduced from (5.39), where the third term is equal to zero at position u = ξ l . The identity (5.39) itself is, as shown above, deduced from the relation (5.35) (which is proven diagrammatically in [13]) by recurrence over the size L of the spin chain. The proof of the second statement, about the equality of coefficients in polynomials on both sides of (5.43), can be achieved by comparing their large u i = u − ξ i asymptotics. First, we can rewrite (5.43) as follows: The r.h.s. must be a linear polynomial in each of the variables u l . We can then reconstruct this polynomial from the large u l asymptotics. For example, for large u 1 we find where we used the fact that at large u 1 the co-derivative in the bracket 1 + 1 u 1 +1−jD ≃ 1 + 1 u 1D can occur only once in each of L! terms of determinant and the Leibniz rule to give the right overall asymptotics u 1 × 1/u 1 = 1. Expanding in this way for each of the remaining u l 's we recover in the r.h.s. the original definition (5.42).
In this way we have proved the determinant formula (5.22). It amounts not only to the proof of the master identity (5.20) but also to the fact that the master T -operator satisfies the whole general Hirota bilinear identity (5.16) for the MKP hierarchy.
Let us note that the commutativity of two general master operators (1.5) follows, in virtue of the determinant representation (5.22), from the commutativity of generating functions of T -operators for the symmetric tensor representations A proof of this last relation can be also found in [13].

Analyticity properties of the BA functions related to the master Toperator
Equations (5.31) also allow us to characterize the class of solutions of the classical Hirota equations relevant to the quantum spin chain. First of all, it is obvious that the master

JHEP09(2013)064
T -operator is a polynomial in u of degree L. Second, consider the (operator-valued) BA function and its adjoint: From (5.31) it is clear that the function ψ is of the form (4.2) while the function ψ * is of the form (4.3) with the poles at the points p i which are the eigenvalues of the twist matrix g (assumed to be all distinct and nonzero). This means that in the appropriate basis g can be diagonalized: The maximal possible order of each pole is L + 1. Indeed, the T -operator for the empty spin chain (at L = 0) gives first order poles at p i : and each co-derivative increases the order of each pole by 1. Note also that (5.51) (and thus the adjoint BA function) has zero of order N at z = 0.

The wave operators as non-commutative generating series
The wave operators ( Expanding both sides of (5.31) in powers of z, we get where χ k (g), χ k (g) are characters of symmetric and antisymmetric representations respectively. From (5.55) it is clear that h a (−∂)T (u, t) = 0 for a > N because the antisymmetric GL(N ) character χ a (g) vanishes. Therefore, the series (5.52) is truncated:

JHEP09(2013)064
Moreover, since for GL(N ) χ N (g) = det g, we have from (5.55): and so the coefficient in front of the last term in (5.56) is equal to (−1) N det g T (u+1,t) T (u,t) as it should, as soon as the points p i are identified with eigenvalues of the matrix g (compare with (4.14) and recall that the τ -function from section 4 and the master T -operator from this section should be identified as τ u (t) = (det g) u T (u, t)).
We conclude that the τ -function (5.9) corresponding to the GL(N )-invariant spin chain on L sites must obey the following two conditions: • It is a polynomial in u of degree L multiplied by exponential function of a linear function of times, • The series for the wave operator (or, equivalently, the expansion of the Baker-Akhiezer function) truncates at the N -th term.
At last, we remark that the wave operator (5.52) at t i = 0 is nothing else than the "non-commutative generating series" of T -operators for antisymmetric (fundamental) representations T a (u) := T (1 a ) (u): Similarly, W −1 (u) is the "non-commutative generating series" 9 of T -operators for symmetric representations T s (u) := T (s) (u) (cf. [3]): The operator T N (u) = det g T ∅ (u+1) = det g φ(u+1) is known as the quantum determinant of the quantum monodromy matrix. These formulas become more simply looking in the normalization (5.5): It would be interesting to clarify the role of the classical Lax operator (3.32) in the quantum spin chain.

Preliminary remarks
We will construct a set of Q-operators which can be identified with the Baxter Q-operators in the sense that: 1. The Q-operators are building blocks for the T -operators meaning that the latter are written in terms of the former in the Wronskian form. The Wronskian formulas resolve functional relations for the T -and Q-operators known as T Q-relations.

The Q-operators satisfy bilinear functional relations (the QQ-relations).
With the additional input of analytical properties of the Q-operators (which are polynomials in u in our case), they imply the system of Bethe equations for their roots.
Our construction of the Q-operators is directly related to the undressing chain of Bäcklund transformations from section 4.2. In our approach, the nested Bethe ansatz is represented by this undressing chain, where the size of the matrices involved in the Wronskian-like solution is successively reduced, as in equation (4.36). As it was shown in section 4.2, the undressing procedure basically consists in picking the singular terms of the function T (u, t + [z −1 ]) (in this section it is an operator) at the poles at z = p 1 , p 2 , . . . , p N . Since the poles are in general of a high order, there are different versions of this procedure according to the possibility to consider different terms in the singular part of the Laurent expansion. As it was argued in section 4.1.5, they are essentially equivalent and differ by certain normalization factors only. Below in this section we consider in some detail two distinguished choices: picking the coefficients in front of the highest and the lowest poles.
As it was already mentioned, the master T -operator introduced in section 5.3.3 is a polynomial τ -function in the variable u while the τ -functions from section 4 are polynomials multiplied by the overall exponential factor N i=1 p u i , so that the definition of the master T -operator differs from the τ -function (4.13) by the factor (det g) −u . More precisely, let |ω be a basis of eigenvectors of the master T -operator and let ω| be the dual basis such that ω|ω ′ = δ ω,ω ′ ). Then where the τ -function in the r.h.s. depends on the state |ω . According to this, the undressing procedure should be slightly modified in order to preserve polynomiality of the τ -functions at each level. In this section we will work in the basis e α where the twist matrix g is diagonal: ge α = p α e α , α = 1, 2, . . . , N :

JHEP09(2013)064
A basis of eigenvectors of g ⊗L = g ⊗ g ⊗ . . . ⊗ g L is then given by where M i is the number of indices α l in the vector |α 1 α 2 . . . α L equal to i. These numbers can be regarded as eigenvalues of the operator 10 In [12,13], it is shown that the T -operators can be explicitly written through diagrammatic expressions, and T (u, t) is then a combination of permutations and diagonal operators. As a consequence, the operators M i commute with all the T -operators. In fact this follows from the GL(N )-invariance of the R-matrix and means that the eigenstates |ω of the T -operators can be classified according to the set of quantum numbers For what follows it is important to note that M i (|ω ) is the multiplicity of zero at z = p i of the eigenvalue of the operator (z I − g) ⊗L on the state |ω . More precisely, the following operator identity holds: (5.62)

The undressing procedure
For our purpose in this section the most convenient version of the undressing procedure is picking the coefficient in front of the most singular term of the Laurent expansion of T (u, t + [z −1 ]) around its poles. As is shown in section 4.2, this version differs from the one based on taking residues by a normalization factor (which is an operator in the present section).
Applying the undressing procedure based on the highest poles, we note that there are N different ways to define the T -operator at the first level of nesting by choosing i ∈ {1, . . . , N } and picking the coefficient in front of the highest pole at p i : 10 Note that mi is a Cartan element eii of gl(N ) in the fundamental representation and Mi is its co-product ∆ (L−1) (eii). In addition, g can be written as with respect to the above eigenbasis.

JHEP09(2013)064
where M i is the operator defined in (5.61). This "operator multiplicity" reflects the fact that the order of the highest pole depends on the sector of the Hilbert space where the diagonalization is performed and is equal to M i + 1. This expression (5.63) is nothing but the equation (4.42) with an extra factor which ensures the polynomiality of T (i) (u, t). Here we introduced a notation | |ω to denote the state dependence of the function.
As the simplest example, consider the case L = 1. From (5.29) we find: As z → p i , the r.h.s. has a double pole in the subspace spanned by the basis vector e i and a simple pole when restricted to the subspace spanned by e α with α = i. In the former subspace the eigenvalue of (z − p i ) M i is z − p i while in the latter it is 1. The explicit calculation of (5.63) gives: The sign factor is (−1) (ij) = (−1) i+j+1 ε ij , where ε ij = 1 if i < j and −1 if i > j.
With this sign convention T (ij) (u, t) is symmetric in the indices i, j, or, equivalently, the transformations T → T (i) → T (ij) and T → T (j) → T (ji) commute for any i, j, namely, T (ij) = T (ji) . The expression (5.67) is nothing but the n = 2 case of (4.43) with an extra factor This factor is necessary to make T (ij) (u, t) a polynomial of the variable u.

JHEP09(2013)064
"nesting level" is defined to be the number n of indices in the set {i 1 , i 2 , . . . , i n }. In complete analogy with the cases n = 1, 2 considered above, we define T -operators T (i 1 i 2 ...in) (u, t) on all intermediate levels by the formula (cf. (4.43)): which leads to where ∆ n is the Vandermonde determinant (4.37) and d n is defined in section 4.2. The T -operator at level n = 0 is the master T -operator defined in (5.26). According to this definition, the correspondence of eigenvalues of the T -operators at level n and the τ -functions at level n given by (4.44) is as follows: The τ -function in the r.h.s. is obtained by the undressing procedure from τ u (t; |ω ) according to the rules presented in section 4. Note that the eigenstate |ω is common for all the T -operators. One can also introduce "Plücker coordinates" for the T -operators on intermediate levels by the expansion whose form is identical to (5.9). The diagrams λ correspond here to representations of the subalgebra gl(N − n) ⊂ gl(N ). It is clear that the bilinear relations of the form (4.39), (4.40) (cf. (5.19), (5.20)) hold for the T -operators T (i 1 ...in) (u, t) at any level and relate T (i 1 ...inij) (u, t) to T (i 1 ...ink) (u, t) and T (i 1 ...in) (u, t): where i, j, k ∈ {1, 2, . . . , N } \ {i 1 , i 2 , . . . , i n }, i = j, i = k, j = k. We thus have the undressing chain that terminates at the level N :

JHEP09(2013)064
It can be shown that the last (the completely undressed) member of this chain is an invertible operator which does not depend on u, t: T (12...N ) (u, t) = T (12...N ) . An advantage of the undressing at the highest poles is that this operator has a simple explicit form: which follows from the comparison of the large u asymptotics of (4.13) and (5.26). Namely, from (4.13) and (4.9) it follows that In the notation of section 4 (see eq. (4.44)), we can write In particular, at J = {1, 2, . . . , N } this yields

JHEP09(2013)064
where we have used (5.76). These Wronskian-like determinants can also be obtained by taking residues of (5.23) at z i k → p i k . The Jacobi identity for the determinant (5.82) gives the so-called QQ-relations, from which the Bethe equations follow. In fact the operators Q J (u, t) introduced in the way explained above are the (specially normalized) T -operators on the intermediate levels of the undressing procedure. The Baxter Q-operators Q J (u) are then defined as their restrictions to zero values of t: We see from the above formulae that Q J (u) is expressed as a residue of the master T-operator. Moreover, if we restrict (5.80) to t = 0, and expand it similarly to (4.10) then we see that the coefficients a im in (4.9) can themselves be expressed from the residues of the master T-operator.
Comparing (5.11) to the equation (4.21) (restricted to t = 0) and using the definition (5.72), we conclude that for the symmetric representation λ = (s) = (s, 0, 0, . . .) the following formula holds: The Jacobi identity for this determinant is known as the T Q-relation. For all other representations this result can be plugged into an analogue of the CBR-formula (1.3): (5.86) to get that the determinant formula (5.85) actually holds for arbitrary Young diagrams with at most N − n rows (if it has more than N − n rows, then the T -operator vanishes identically). In the particular case I = ∅ equation (5.85) yields the Wronskian determinant representation of the T -operator at fixed representation: At L = 0 the pre-factors and the operators Q k (u) become simple (equal to 1/p j ) and this formula coincides with the determinant representation (3.21) for characters. By construction, all the Q-operators are polynomials in u and the determinant formula (5.87) at λ = ∅ must give the fixed polynomial φ(u) = L j=1 (u − ξ j ). This requirement implies the system of Bethe equations for roots of the Q-operators (or rather their eigenvalues).
The undressing procedure is schematically illustrated in figure 2 for the case N = 2. The green double arrow corresponds to the Bäcklund transformation (5.91), which "undresses" the T -operators by gradually decreasing the size of the matrices at each site of the lattice (which is the maximal possible height for the Young diagrams giving rise to nonzero T -functions). The blue arrow shows how the original T -operators are reconstructed by "dressing" with the Q-operators.

On other versions of the undressing procedure
The undressing procedure and the Q-operators can be defined in a number of equivalent ways. Here we mention some of them.
Lowest poles. Instead of picking the highest order pole, one might consider another version of the undressing procedure, based on taking residues at other poles (not of the highest order), as in section 4.2. For example, at the first level (cf. (4.30)), we may define: This definition does not use the operators M i , at the price of less transparent normalization factors. In particular, the definition (5.88) with the unit pre-factor leads to a more complicated operator T (12...N ) .
Construction of [12]. At last, let us compare our present construction with the one introduced in [12]. On the first level, noticing that (z − p i ) M i = The normalization factor N i is independent of u and t, and commutes with all the Toperators. For this reason it was disregarded in [12]. On higher levels (5.70) can as well be written as: where N i 1 ...in is an operatorial normalization factor which we don't specify here. Let us now show that (5.91) coincides with equation (32) from [12].
1. In [12], the set of times t was not introduced. But as it can be seen from (5.31), the transformation t → t + [z −1 ] is equivalent to the insertion of w(z −1 ) = (det(I − z −1 g)) −1 to the right of the co-derivatives.
2. The limit z j → 1/x j used in [12] has the same effect as taking residue at the point z j = p j if we change notations z → 1/z, p j → x j , and notice that res z=p f (z) = lim z→p (z − p)f (z) for a function f with a simple pole.
3. In [12], depending on the representation, an overall factor is necessary to convert characters of g = diag p i | 1≤i≤N ∈ GL(N ) into characters of g J = diag p j | j∈J ∈ GL(|J|). In our current notations, this amounts to changing F g (t) to F g J (t). One can easily see that this is exactly what the factor e − n α=1 ξ(t,p iα ) does. Indeed, e − n α=1 ξ(t,p iα ) = F g J (t)/F g (t).
This comparison shows that the undressing procedures of [12] and the one of the present paper are essentially equivalent. But in the setting of our present paper this undressing procedure has a transparent connection to the classical integrable hierarchy and its polynomial τ -functions. Our construction clarifies the combinatorial structure of the set of the 2 N Q-operators on the Hasse diagram found in [68]. Moreover, it naturally takes into account the fact that all these Q-operators, as well as all the T -operators, are polynomials in the spectral parameter u. As a consequence, the Bethe equations are automatically contained in this analysis, as it was explicitly demonstrated in [12].

Dynamics of zeros of the master T -operator
where the roots u j depend on all other times t i (we have omitted an irrelevant exponential function in front of the product). In some special cases 11 one or more roots u j can tend to infinity and so the degree of the polynomial can be less than L. Note that we treat equation (5.92) as an operator relation and the roots u j (t) are (commuting) operators.
Another possibility is to treat (5.92) as a relation for eigenvalues of the master T -operator, then the roots of each eigenvalue have their own dynamics in the times t i . This dynamics is known [26,27] to be given by the rational Ruijsenaars-Schneider model. A solution of the spectral problem for the master T -operator in terms of the integrable dynamics of its zeros and a connection with Bethe equations will be given elsewhere. As is well known, the Ruijsenaars-Schneider model contains the Calogero-Moser system of particles as a limiting case. In this connection let us note that a similar relation between the quantum Gaudin model (which can be regarded as a degeneration of quantum spin chains or vertex models with rational R-matrices) and classical Calogero-Moser system was found in [69,70] from a different reasoning. An important remark is in order. As is seen from (5.3), the inhomogeneity parameters ξ i are coordinates of the Ruijsenaars-Schneider particles at t i = 0: ξ i = u i (0). Clearly, they provide only a half of the initial data necessary to fix the dynamics. The other half is initial values of momenta which are implicitly determined by the above mentioned conditions on the τ -function (5.9). Since the quantum transfer matrix for the spin chain is of the size N L ×N L , it has not more than N L different eigenvalues. This means that the system of equations that fix the initial momenta has only a finite number of solutions, corresponding to common eigenstates of the transfer matrices.

Concluding remarks
Our main result in this paper is the identification of the generating function of commuting integrals of motion in generalized quantum spin chains (called here the master T -operator) with the τ -function for integrable hierarchies of classical soliton equations. Among other things, this result means that quantum integrable models have a hidden free fermionic (Grassmann) structure in the auxiliary ("horizontal") space. It would be very interesting to relate it to the hidden Grassmann structure uncovered in [71][72][73][74] in the quantum ("vertical") space of quantum integrable models. Another context where τ -functions of classical integrable hierarchies enter the theory of quantum integrable models and associated 2D lattice models of statistical mechanics is calculation of scalar products and partition functions with domain wall boundary conditions [75][76][77].
Many things remain to be done in order to further clarify and generalize this classicalquantum correspondence and, possibly, to use it for inventing new methods to solve both kinds of problems. Here we list some obvious directions of further development: 11 For example, when the boundary twist parameters tend to 1.

JHEP09(2013)064
-To construct the Bethe vectors within this formalism (which does not use the standard algebraic Bethe ansatz) and to clarify their role in the classical MKP hierarchy and its polynomial reduction according to Krichever.
-To extend the results to quantum integrable models with trigonometric R-matrices. This is more or less straightforward and clear how to do. Whether the approach based on the co-derivative is applicable is not quite clear at the moment, however. Instead of rational solutions of the MKP hierarchy one should consider solutions with trigonometric dependence on the spectral parameter. The Krichever's construction is directly extendable to this case. Zeros of the master T -operator move as particles of the trigonometric Ruijsenaars-Schneider model.
-To extend the results to quantum integrable models with elliptic (Belavin's) Rmatrices. This is the most general and most interesting case. Zeros of the master T -operator move as particles of the elliptic Ruijsenaars-Schneider model. The relevant class of MKP solutions is a special subclass of general elliptic (two-periodic) solutions. One of the interesting problems here is to characterize this subclass in algebro-geometric terms through the corresponding algebraic curves. Elliptic solutions in general position correspond to smooth curves which are coverings of elliptic curves. Presumably, the solutions relevant for quantum spin chains are not of general position and the corresponding algebraic curves are singular.
-To extend the results to supersymmetric GL(N |M )-invariant integrable spin chains. The general definition of the master T -operator and the Hirota equation remain the same. The realization through the co-derivative is known from [12], and can certainly be related to Krichever's rational τ -function.
-The fact that the Baxter Q-operators and T -operators on intermediate levels of the nested Bethe ansatz are contained in the master T -operator seems to be general and not restricted to models with rational R-matrices. One can try to find a general procedure to obtain these objects from the master T -operator applicable in all the cases listed above. The example of rational solutions suggests that such a procedure should be formulated in terms of Baker-Akhiezer functions on singular algebraic curves.
-A development of this work might also contribute to representation theory of quantum affine algebras and, in particular, to the theory of q-characters introduced by Frenkel and Reshetikhin [47] as a spectral parameter dependent version of usual characters.
-We also expect that our result can be extended to quantum models associated with more general Lie algebras. This is in accordance with the fact that the KP hierarchy has analogs for Lie algebras other than A N −1 [22].
-Our classical interpretation of the quantum integrability should be also applicable to the 2-dimensional integrable quantum sigma models, where the operatorial formalism is almost unexplored. The most interesting case of such applications could be JHEP09(2013)064 the construction of an operatorial form of the AdS/CFT Y-system (and the associated T-system) [49,[78][79][80] and clarification of its relation to the classical integrable hierarchies.
A more ambitious goal is to suggest an alternative approach to diagonalization of quantum transfer matrices which would not use the Bethe ansatz at any step. This does not seem hopeless because the quantum problem appears to be reducible to a purely classical integrable dynamics of Calogero-like systems of particles.
In a similar way, one can prove that the determinant formulas of the Jacobi-Trudi type hold for arbitrary Young diagrams. Let us rearrange (3.16) as follows: By using the relations ψ * m ψ m = 1 − ψ m ψ * m and n−β 1 −1| ψ m = 0 for m ∈ {n−β j −1} Taking into account the following relations for the Frobenius notation: , we can rewrite this as follows: Applying the Wick's theorem in the same way as for the λ = (α|β) case, we finally obtain the desired result (3.18).

JHEP09(2013)064
C Proof of the quantum Giambelli formula The quantum Giambelli formula was proposed in [44] in relation to the analytic Bethe ansatz for finite dimensional U q (B (1) N ) modules. It was also mentioned there that this formula follows from Sylvester's theorem without detailed calculations. To make the text self-contained, we will give a proof of (5.7) based on the CBR formula (5.6). Here we follow [28, page 88, example 22], where the classical case (without spectral parameter) is discussed. In this section, we will use Frobenius notation for a Young diagram mentioned in section 2.2. We introduce the following notation: We also introduce tuples I = (−α 1 , . . . , −α d , 1, 2, . . . , n), J = (1, 2, . . . , n). We will use the following relations for the Frobenius notations: Then (5.6) can be written as Next we introduce the following matrices: where I 1 = (−α 1 , . . . , −α d ), J 1 = (1, 2, . . . , d), and I d×d is a d × d unit matrix and (0) d×n is a d × n zero matrix. Then we find Let us replace j-th row of A with i-th row of B and write the determinant of this matrix as c i,j . We will consider a matrix C = (c ij ) 1≤i,j≤d+n . Taking note on the relations c ij = δ ij det A for d + 1 ≤ i, j ≤ d + n and c ij = 0 for d + 1 ≤ i ≤ d + n, 1 ≤ j ≤ d, we obtain det 1≤i,j≤d+n We will use det A = n j=1 t j,j . Expanding the j-th row of the determinant c i,j , we obtain Then we obtain an identity det 1≤i,j≤d

D The co-derivative
In this appendix we summarize formulas on the co-derivative introduced in [13]. The (left) co-derivative is the operator that acts on any function f of a group element g ∈ GL(N ). Before we provide its complete definition (in terms of components), let us introduce it symbolically as follows: φ α β e αβ is a general element of the Lie algebra gl(N ) parametrized by φ α β , with e αβ being the matrix with only one nonzero (equal to 1) matrix element at the place (αβ): (e αβ ) µ ν = δ αµ δ βν . To put this more precisely, let us start from scalar functions. In this case, we do not need the symbol ⊗ and get:

JHEP09(2013)064
where I is N × N unit matrix. If f is a N × N matrix valued function, the result of (D.1) belongs to N 2 × N 2 matrix. In particular, if f (g) = g, we have in components: In invariant form, this formula can be written aŝ where P is the permutation operator defined as P(A ⊗ B) = (B ⊗ A)P for any N × N matrices A, B. Explicitly, we have P = N a,b=1 e ab ⊗ e ba . In the main text, we considered actions of the co-derivative on functions in (C N ) ⊗L . In this case, it is convenient to define the co-derivative in the following way. Let us consider an embedding of any matrix A on C N into (C N ) ⊗L , and introduce the following notation: (D.5) The action ofD i on any function of g ∈ GL(N ) in (C N ) ⊗L is defined aŝ For example for f (g) = g j , we obtain a generalization of (D.4): where P ij = ab e (i) ab e (j) ba . For i = j case, the permutation operator reduces to the second order Casimir, and in this case, it is P ii = ab e (i) ab e (i) ba = ab e (i) aa = N I ⊗L since we are considering the fundamental representation. Then we obtainD i g i = N g i . In particular for L = 2 case,D 1 g 2 = P 12 g 2 stands for (D.4), whileD 2 g 1 = P 21 g 1 stands for (I ⊗D)(g ⊗ I) = P(g ⊗I). Now the Leibniz rule for any functions A(g), B(g) of g ∈ GL(N ) can be written as where I = (i 1 , i 2 , . . . , i n ), J = (i α 1 , i α 2 , . . . , i α k ), K = (i β 1 , i β 2 , . . . , i β n−k ) (α 1 < · · · < α k ; β 1 < · · · < β n−k ) are, as sets, subsets of {1, 2, . . . , L}, which satisfy I = J ⊔K, J ∩K = ∅.

JHEP09(2013)064
The ordered product on any indexed operator {A i } i∈I is defined as − → i∈I A i = A i 1 · · · A in . The Leibniz rule can be used, for example, {i,j,k} (D i g j )(D i g k ) + C (n−1) {i,j,k}D iDj g k = C (n−2) {i,j,k} P ij g j P ik g k + C (n−1) {i,j,k} P jk P ik g k (D.12) Let S (n) (J) be a subgroup of the permutation group over the set J whose elements have exactly n cycles: any element σ of S (n) (J) has a cycle decomposition σ = τ 1 · · · τ n ({τ i } are mutually disjoint cyclic permutations on J). Then we have the following formula Here P σ is a matrix in the space (C N ) ⊗L , with the matrix elements [P σ ] i 1 ...i L j 1 ...j L = δ i 1 j σ(1) · · · δ i L j σ(L) .
The action of σ ∈ S (n) (J) on i ∈ {1, 2, . . . , L} \ J is formally defined as σ(i) = i. The summations are taken over all the decompositions of {1, 2, . . . , L} into two disjoint sets I and J. The above identity can be proven by induction in L and using the Leibniz rule.
Yang-Baxter equation and commutation of co-derivatives. One can check that these co-derivatives obey the following commutation relation D iDj −D jDi = (D i −D j )P ij = −P ij (D i −D j ), (D.14) where the last equality comes fromD j P ij = P ijDi .

JHEP09(2013)064
This relation can be understood from the following Yang-Baxter equation : (u τ (k) +D τ (k) ) (u j − u i + P ij ) (D. 15) where τ (k) = For instance, if we write (D.15) when L = 2, and u 2 = 0, then the coefficient of degree one in u 1 gives exactly (D.14). Multiplying (D.15) by P ij from the left, we can rewrite (D.15) in the notation of the tensor product used in the main text: The equation (D.15) can either be proven "from scratch" by proving explicitly the property (D.14) of co-derivatives, or it can be deduced (when the co-derivatives act on linear combinations of characters, which is sufficient for this paper) from the usual Yang-Baxter Identity.

E Proof of (5.39)
The left hand side of (5.39) is a polynomial of u 1 , u 2 , . . . , u L . The degree of this polynomial with respect to each u i (1 ≤ i ≤ L) is 0, 1 or 2. Taking note on this fact, we will prove (5.39) using the induction with respect to the length of the chain L. The proof that (5.39) holds for L = L 0 , consists in the following three steps: (i ) First we will see that the terms of degree two in one of the variables u i are zero under the hypothesis that (5.39) holds for a chain with length L = L 0 − 1.
This step is easy since the term of degree two in u j is obtained by the substitutionŝ D L 0 (u) . We see that the operator obtained from the action of these co-derivatives is a trivial operator on the j th space times an operator acting on a spin chain with L 0 − 1 spins (corresponding to all sites except the site j). We can then recognize that the coefficient of u 2 j in the left-hand side of (5.39) at L = L 0 is equal to the left-hand side of (5.39) at L = L 0 − 1 (up to a relabelling of the spaces), which is zero by the recurrence hypothesis.
(ii ) The other terms have degree one in some variables u i 1 , u i 2 , . . . , u in and degree zero in the other variables u j 1 , u j 2 , . . . , u j L 0 −n (where {j 1 , . . . , j L 0 −n } denotes the complement of {i 1 , . . . , i n }). We can show that all the terms where 1 ∈ {i 1 , . . . , i n } vanish, assuming that (5.39) holds for a chain with length L = L 0 − 1.

JHEP09(2013)064
To perform this step, we have to investigate the coefficient of degree one in u 1 , in each term of the left-hand-side of (5.39). For instance, in the operator D L 0 (u + 1)w(z −1 1 ) D L 0 (u)w(z −1 2 ) , the term of degree one in u 1 is given by where the right-hand-side is obtained by the Leibniz rule. The same analysis is easily performed for the other terms D L 0 (u + 1)w(z −1 2 ) D L 0 (u)w(z −1 1 ) and D L 0 (u + 1)w(z −1 1 )w(z −1 2 ) D L 0 (u) , (appearing in (5.39)) and we obtain that the coefficient of degree one in u 1 is obtained by the action of (I+D) on the left-hand-side of the identity with L = L 0 − 1 spins. Hence the recurrence hypothesis shows that they vanish.
(iii ) Finally, we will deduce the vanishing of the term of degree one in the variables u i 1 , u i 2 , . . . , u in and degree zero in the variables u j 1 , u j 2 , . . . , u j L 0 −n , where {j 1 , . . . , j L 0 −n } denotes the complement of {i 1 , . . . , i n }. We will prove this vanishing by recurrence over n (i.e. we will assume that it holds for n − 1, and we will deduce it for n).
First, let us assume that the terms with degree one in u k and degree zero in u k+1 (i.e. the terms where k ∈ {i 1 , . . . , i n } and k + 1 ∈ {j 1 , . . . , j L 0 −n }) do vanish, and let us deduce that the term with degree zero in u k and degree one in u k+1 (i.e. the terms where k ∈ {j 1 , . . . , j L 0 −n } and k + 1 ∈ {i 1 , . . . , i n }) also vanish. To this end, we can use a Yang-Baxter equality (D.17) for each term in the identity (5.39). For instance for the first term (i.e. for z 2 D L 0 (u + 1)w(z −1 1 ) D L 0 (u)w(z −1 2 ) ) we write

JHEP09(2013)064
If we sum this Yang-Baxter identity for each term in (5.39) and we keep only the terms of degree one in the variables u i 1 , u i 2 , . . . , u in and degree zero in the variables u j 1 , u j 2 , . . . , u j L 0 −n , then we get an identity for the sum of the following terms • First, the terms which contain P k,k+1 . As the total degree in u k is one and the total degree in u k+1 is zero, these terms do not have any occurrence of u k or u k+1 in the factors containing co-derivatives. One can easily see that (up to multiplication by P k,k+1 ) these terms are exactly the terms of (5.39) with degree zero in the variables u j 1 , u j 2 , . . . , u j L 0 −n and u k , and degree one in the other variables. By the hypothesis of the recurrence over n, these terms vanish.
• Next, the terms of the first line in (E.2) which do not contain P k,k+1 . If we sum them for the three terms of (5.39), we exactly obtain the terms of degree one in the variables u i 1 , u i 2 , . . . , u in and degree zero in the variables u j 1 , u j 2 , . . . , u j L 0 −n of (5.39) (with k ∈ {i 1 , . . . , i n } and k + 1 ∈ {j 1 , . . . , j L 0 −n }). By hypothesis these terms do vanish.
• Finally, the terms of the second line in (E.2) which do not contain P k,k+1 . If we sum them for the three terms of (5.39), we exactly obtain the coefficient of degree one in u k+1 and degree zero in u k (and with degree one, resp zero in the other variables i m resp j m ). These last terms vanish (because the sum of expressions of the form (E.2) is zero, and because the other terms vanish). Here we used the fact that the coefficient of u 1 k u 0 k+1 in a polynomial f (u σ(1) , u σ(2) , . . . , u σ(L) ) coincides with the coefficient of u 0 k u 1 k+1 in a polynomial f (u 1 , u 2 , . . . , u L ), where f (u 1 , u 2 , . . . , u L ) is any polynomial in u 1 , u 2 , . . . , u L . This shows that if the terms with degree one in u k and degree zero in u k+1 (i.e. the terms where k ∈ {i 1 , . . . , i n } and k + 1 ∈ {j 1 , . . . , j L 0 −n }) do vanish, then the term with degree zero in u k and degree one in u k+1 (i.e. the terms where k ∈ {j 1 , . . . , j L 0 −n } and k + 1 ∈ {i 1 , . . . , i n }) also vanish.
The same result is easily obtained if replace k + 1 with k − 1, and by iterating this argument, we obtain that the term of degree one in any set of n variables u i 1 , u i 2 , . . . , u in (and degree zero in the other variables) do vanish. Indeed, the above argument allows to reduce these term where 1 ∈ {i 1 , . . . , i n }, treated above. This proves the equation (5.39) by recurrence.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.