Canonical formulation for the thermodynamics of $sl_n$-invariant integrable spin chains

Integrable quantum spin chains display distinctive physical properties making them a laboratory to test and assess different states of matter. The study of the finite temperature properties is possible by use of the thermodynamic Bethe ansatz, however at the expense of dealing with non-linear integral equations for, in general, infinitely many auxiliary functions. The definition of an alternative finite set of auxiliary functions allowing for the complete description of their thermodynamic properties at finite temperature and fields has been elusive. Indeed, in the context of $sl_n$-invariant models satisfactory auxiliary functions have been established only for $n \leq 4$. In this paper we take a step further by proposing a systematic approach to generate finite sets of auxiliary functions for $sl_n$-invariant models. We refer to this construction as the canonical formulation. The numerical efficiency is illustrated for $n=5$, for which we present some of the thermodynamic properties of the corresponding spin chain.


Introduction
Almost a century has now passed since Bethe's seminal work in which he first transformed the solution of the one-dimensional Heisenberg spin chain into a system of nonlinear algebraic equations [1].The understanding of this model and of other similar systems has developed since.Lieb identified the same eigenvectors of the Hamiltonian as belonging to a transfer matrix of a classical statistical model [2], and not long after Baxter put forward the idea of commuting transfer matrices [3] which naturally become generators of conserved charges.By now, Yang-Baxter integrable spin chains and other exactly solvable models form a field of their own.The question remains to what extent quantum integrability is an asset in determining physical and mathematical properties.
In this respect, the treatment of fully interacting integrable systems is a much more intricate task than for the case of free fermion models.Integrability imposes factorization of the multi-particle S-matrix into 2-particle S-matrices, however with non-trivial scattering.To mention from the many physically relevant properties a most important example, consider the evaluation of correlation functions for arbitrary distances, fields and temperature which is, even from a wide point of view, successful so far only for a few integrable systems [4,5].Each of the most prominent integrable models inspired the development of a range of different (and sometimes model-specific) techniques, which is very much exemplified by the several formulations of Bethe ansatz that emerged so far [6][7][8][9][10][11][12].On one hand it evinces ingenuity in the field, on the other it demonstrates a remaining lack of a universal comprehension.
Similarly, the evaluation of thermodynamical properties, like other physical properties, has not escaped this situation.In principle, once a model is identified as Yang-Baxter integrable and a Bethe ansatz is established, all what is left is to solve the Bethe equations in order to obtain the eigenvalues for evaluating the partition function.This is the basis of the so-called Thermodynamical Bethe Ansatz (TBA) [13][14][15].Such a program leads to the study of all patterns of Bethe roots corresponding to bound states, which culminates into the string-hypothesis for integrable quantum chains [15].An alternative procedure, known as the Quantum Transfer Matrix (QTM) method, allows for the evaluation of the thermodynamical potential through a single eigenvalue of a new object [16][17][18][19].
Among the criticisms to the former method are the lack of proof of the string hypothesis and the fact that one needs to solve an infinite number of nonlinear integral equations (NLIEs) to obtain the desired properties.Therefore, the numerical evaluation relies on a truncation procedure that is hard to justify a priori.The second program has been highly successful in a number of cases [20][21][22][23][24][25][26].Besides the fact that no string hypothesis is required, the resulting nonlinear integral equations are finite in number, which allows for an efficient numerical evaluation.The drawback of this approach is that no systematic way to produce such sets of nonlinear integral equations has been found so far, since one needs to define suitable auxiliary functions whose construction can be very puzzling.
In this paper we apply the QTM method for sl n -invariant models and propose a set of auxiliary functions which are solutions to a finite set of NLIEs (Section 4).Our choice is based on bilinear relations similar to the T-system [27,28] with one additional advantage: we find factorization of the symmetric fused terms, which implies the truncation to a finite number of independent functions (Section 3).We prove these relations through manipulations of the Yangian version of Young tableaux (Section 2).In the special case of sl 5 (Section 5), we present a numerical evaluation of certain thermodynamical quantities for a model of hard core bosons resp.itinerant pairs of electrons on a chain.

Setting the stage 2.1 The Quantum Transfer Matrix
The transfer matrix is a fundamental object in classical statistical models.For a homogeneous system with degrees of freedom on the lattice one can compute the partition function of the model by arranging its Boltzmann weights into a matrix (Lax operator).Then, by performing simple algebraic operations like taking products of successive copies followed in case of periodic boundary conditions by a trace in the auxiliary space, one obtains the transfer matrix.The largest eigenvalue of this yields the bulk properties of the model replacing the calculation of an infinite sum over all configurations.More specifically, in the case of two-dimensional lattice models (such as vertex models, for example) one assigns independent degrees of freedom to all bonds of the lattice, thus attaching a tensor product structure to the Hilbert space.The Lax operator L is defined as where e (j) are tensorized Weyl matrices, acting only on the quantum (or vertical) space V j , j = 1, . . ., L.
Then, the transfer matrix is given by where we took the transfer direction from row to row and V A denotes the auxiliary (or horizontal) space.
Integrability is realized by the commutativity of row-to-row transfer matrices.One way to see this is to construct a set of Boltzmann weights parameterized by a (complex) spectral parameter such that for arbitrary values λ, µ [T (λ), T (µ)] = 0, (2.3) and hence T (λ) serves as a generating family for the conserved charges through its derivative.The latter condition is global, while one only needs to know local Boltzmann weights to produce transfer matrices.For Eq.(2.3) to hold a sufficient condition is the celebrated Yang-Baxter equation Under the assumption that the model is fundamental, i.e. all local spaces are isomorphic, we make no distinction between the intertwiner (R-matrix) and the Lax operator.Furthermore, for simplicity, in what follows we will assume a set of desirable properties, namely Time reversal: Parity: [L 12 (λ, µ), P 12 ] = 0, (2.7) where P ij denotes the permutation operator and t k the transposition in the k-th space.A consequence of Eq.(2.4) in case of the standard initial condition (regularity) for isomorphic spaces is the unitarity relation Also thanks to the standard initial condition, one can find a local Hamiltonian through the logarithmic derivative of the row-to-row transfer matrix, . (2.9) A set of solutions satisfying all properties above corresponds to the Perk-Schultz model [29].The Lax operator in this case reads which satisfies where G is the first fundamental representation of any generator of the sl n algebra.
In principle, one can use the Bethe ansatz to obtain the eigenvalues and the universal eigenstates of the integrable manifold.However, in order to evaluate the partition function of the spin chain one needs to sum over all states.This task is best accomplished by following an alternative route.
The statistical operator e −βH appears in the expansion of the row-to-row transfer matrix (2.12) Now define a conjugated transfer matrix as so that it admits a similar expansion with the shift operator reversed Consequently, By choosing λ = − βJ N , terms O(λ 2 ) become small compared to O(λ) as we take the Trotter limit N → ∞.Therefore we find the Trotter-Suzuki decomposition, where we have included n charges of the sl n algebra, here denoted by Nj , which satisfy the constraint n j=1 Nj = L.The coupling constants µ j play the role of generalized chemical potentials.Note that the product T (λ) T (λ) may be understood as the transfer matrix of a staggered lattice model.
Instead of reading the products that make up the partition function in the usual row-to-row direction one can read them from column to column.This gives rise to the quantum transfer matrix, (2.17) The new spectral parameter x is introduced such that commutativity holds T QT M (x), T QT M (y) = 0.With this, the partition function can be written as Now the single largest eigenvalue, Λ QT M max (x), is needed to evaluate the thermodynamic potential.Indeed,

Yangian analogue of Young tableau
The quantum transfer matrix inherits the integrability structure from its row-to-row counterpart.Both matrices possess the same intertwiner.Nevertheless, because of the alternation that defines Eq.(2.17), the reference states and the sectors associated to sl n charges are different.Here where nk,j = 0, 1 depending on whether site j is occupied by species k or not, i.e. nk,j are projectors onto orthogonal 1-dimensional spaces.On the account of n k=1 nk,j =id j we have the constraint Consequently, in terms of Bethe ansatz expressions, only vacuum expectation values and sectors are changed.Choosing the reference state |1, n, 1, n, . . ., 1, n⟩, the eigenvalue is given by where the eigenvalue functions λ j (x), j = 1, . . ., n are given by [30] where k , k = 1, . . ., m j are solutions of the Bethe ansatz equations, given in the form lim We are specially interested in the sector m 1 = m 2 = . . .= m n−1 = N 2 , which corresponds to Q k = 0 for all k, where the largest eigenvalue is found.
In Eq. (2.22) we made use of the Yangian version of the Young tableau [31][32][33].A rectangular tableau of width s and height a, where s, a are non-negative integers represents the product of functions (2.22) where the spectral parameter of the corresponding functions increases from top to bottom and decreases from left to right on the tableau in steps of i.It is important to observe the standard filling rules when building a tableau, namely j k,l ≤ j k,l+1 and j k,l < j k+1,l .
Finally, note that the sum over all possible fillings of Eq.(2.25) corresponds to the eigenvalue Λ (n) a,s (x) of the fused transfer matrix with anti-symmetric index a = 1, . . ., n − 1 and symmetric index s = 1, . . ., ∞.

Deductive proof of T-system and fusion hierarchy
Here we lay the groundwork to derive auxiliary functions in the forthcoming sections.The auxiliary functions can be written conveniently in terms of Young tableaux introduced above.They come in two sets, here denoted by uppercase {B j } and lowercase {b j } letters and are related to each other by B j (x) = 1+b j (x).Our goal is to find nontrivial multiplicative functional relations between the lowercase and uppercase functions.
Our definition is inspired by the well-known bilinear relations [27] among fused eigenvalues, which go by the name "T-system".
Although it is not our intention to discuss completeness of Bethe ansatz eigenvectors, we may assume that at least a few eigenvalues are given in terms of an analytical function in the manner we described above.
In this sense, if we find functional relations in terms of the Yangian analogue of the Young tableaux, it is reasonable that some sectors of the transfer matrix if not all do satisfy these functional relations, including the leading eigenvalue.This is certainly the case if the relations are derived by, for example, the fusion algebra for operators.
We want to show the relation As an example let us take two tableaux whose spectral parameters are shifted from each other by one in units of the crossing parameter.Let a be the number of rows and s = 1 the number of columns of these tableaux.They are meant to represent the product of T (a,1) x − i 2 T (a,1) x + i 2 .From Eq.(2.26) it is clear that in the case a = 4 we would like to find where an empty tableau is implicitly understood as the sum over all possible fillings.In the following we drop the explicit reference to the spectral parameter whenever no ambiguity is possible.One way to do so is to highlight at least the boxes that share one specific spectral parameter and then preferably to align their horizontal positions, as exemplified above.Before we jump into Eq.(2.27) per se let us look into the most fundamental fusion relation, which is actually the special case a = s = 1.We have (2.28) Indeed, say the white and colored boxes on the left hand side of Eq.(2.28) are filled with any non-negative integer indices j and j h respectively.Since either j < j h or j h ≤ j, which are precisely the filling rules for Young tableaux of width 1 and height 1, respectively, it is clear that the summation over all indices leads to the right hand side of Eq.(2.28).Now let us consider the left hand side of Eq.(2.27).There are several possibilities of rearranging the boxes according to their fillings, for example or which correspond, respectively, to the conditions: Note conditions i)-iv) satisfy all allowed filling rules for a combination of two column tableaux of sizes 3 and 5, i.e. the first term on the RHS of Eq.(2.27), while condition v) describes exactly the filling of a tableau of width 2 and height 4, i.e. the second term on the RHS of Eq. (2.27).What is left to show is that conversely any term appearing in the first term on the RHS of Eq.(2.27) appears in one of the cases i)-iv).
This can be done which we leave to the reader.This shows that the sum of the cases i)-iv) amounts to the first term on the RHS of Eq.(2.27)where each of the two columns can be summed over independently.This completes the proof of identity Eq.(2.29).One can apply the same reasoning recursively to show Eq.(2.26) for arbitrary values of a, s.
In the Yangian analogue of the Young tableau, boxes can only exchange positions when they have the same spectral parameter.This helps to pinpoint the possible outcomes in relations like the above.For instance, on the basis of the same analysis we may find It is then a pleasant game to use these "fusion moves" in order to derive functional relations among transfer matrices.

Adjacency matrices
There are circumstances where the set of functional relations provided by the T-system truncates by itself.
This can be viewed as a generalization of the inversion identity [27,34], which can either be used to derive finite sets of NLIE or to describe the spectrum in terms of the zeros of the eigenvalue functions.Nevertheless, in the general case such a "spontaneous" truncation does not take place, which makes it more natural to represent spectral values in terms of (fused) eigenvalues and Bethe ansatz, where extra polynomials like the q j (x) functions also enter explicitly.More generally, other types of functions may appear in different sets of functional equations for suitable auxiliary functions based on Young diagrams with restricted filling rules.In this section, we obtain these functions, henceforth called EAFs for elementary analytic functions or elementary analytic factors.

Analytical consequences of Bethe ansatz for fundamental representations
The eigenvalue Eq.(2.21) as function of the spectral parameter may have poles at points related to the zeros of the functions q j (x).Since the eigenvectors of the transfer matrix do not depend on the spectral parameter, the corresponding eigenvalues like all matrix elements must be analytical functions of x.This way one recovers the Bethe ansatz equations Eq.(2.24) whose zeros are those of the q j (x), j = 1, . . ., n − 1.
For the first fundamental representation (a = 1) the removability of these poles happens pairwise through the sum λ j (x) + λ j+1 (x).On the other hand, for the cases in which the sum of tableaux is truncated, due to restricted filling rules, not all poles are removed.To see this, one must first study case-by-case the relationship between any two tableaux that belong to the same Young diagram resp.representation and check whether they are connected via a Bethe ansatz equation.Then we encode this information into socalled adjacency matrices from which later we define the additional EAFs that appear in the definition of the auxiliary functions with the proper poles.
The adjacency matrices are represented by graphs where each vertex corresponds to a term of the eigenvalue Λ (n) a,1 and each edge to the shared pole between the vertices.In Fig. 1 we present the graphs of sl 4 with s = 1 and a = 1, 2 as examples.For the first fundamental representation the graph would simply amount to the A n−1 Dynkin diagram.See Fig. 1(a).The situation becomes more intricate as we move to fundamental representations with a > 1.
While a graphical approach consisting in identifying the proper EAFs from cuts of these graphs is possible, for representations with a > 3 this would become unfeasible on the account that these graphs spread to higher dimensions.To this end the adjacency matrices are more versatile.Nevertheless we refer to Appendix A in which we show how the graphs are used for the sl 4 case.
We enumerate the vertices by ordering them resp.the corresponding tableaux such that the tableau with smaller indices precedes the rest.For instance, in the case of the a = 2 representation of sl 4 we set Having established this ordering, set the entries of the adjacency matrix at position (i, j) to 0 if the vertices No. i and No. j are not connected or to the corresponding q k (x) of the removed set of poles if they are connected.These matrices have to be worked out individually.As an example, the matrices for all Figure 1: Pole removal graphs for the fundamental representations of sl 4 .The q-functions on the bonds refer to the shared pole between the vertices it connects.We use the notation: fundamental representations of the sl 4 invariant model with s = 1 are given by

Identification of EAFs and notation
As we can see from the graphs, the Bethe equations imply analyticity for any eigenvalue function, i.e. for the sum over all terms represented by a certain Young tableau.But also partial sums over such terms show the phenomenon of pole cancellation due to the Bethe equations, however with certain poles not being removed.The partial sums appear in factorizations of the later to be introduced auxiliary functions.It will be important to write the partial sums in terms of analytical functions divided by products of q-functions with certain shifts of the arguments.The use of the Fourier transform for the logarithmic derivatives of all functional equations will then allow us to derive non-linear integral equations for the auxiliary functions.
In general, one may define the analytical factor functions by summing over tableaux, checking from the adjacency matrix which poles are being removed and which are not and establishing the combination as a polynomial function p(x).Here, certain index, shift and sign functions ι : and s : β → ± appear.The label α indexes the set of poles of individual terms without cancellation in the sum, β indexes the factors of type Φ ± appearing jointly in all terms of the sum thanks to definition Eq. (2.22).The functions ι and s denote the type of poles and zeros, the function δ gives the shift of the arguments of the elementary factors.Then, the nonremoved poles -namely the ones appearing as elements of the adjacency matrices in rows of any term of the sum in question, but in a column of a term not appearing in that sum -result in the numerator of Eq.(3.3).The factors Φ − and Φ + appearing in the denominator originate from factoring out q 0 = Φ − or q n = Φ + jointly appearing in all terms of the sum.The eigenvalues a,1 (x) are also EAFs, obtained when we consider the full matrix.Of course, in this case there are no poles left behind.
For instance, take the first fundamental representation of sl 4 .The sum 1 + 2 removes the poles at q 1 (x) but keeps the poles at q 2 (x) (see 1,1 (x)) above.Moreover common trivial zeros may be factorized and one can show that is an entire function.More specifically, it is a polynomial of degree N 2 + m 2 .On the other hand, as an example for the second representation, we may use It is worth noting that representation matrices with larger a may give the already obtained EAFs as well as new ones.Also, the matrix A (n) n−1,1 (x) does not contribute with new EAFs for the sl n case except for Λ The polynomials p (n) a,k (x) together with q 1 (x), . . ., q n (x) and the eigenvalues Λ a,1 will be called the EAFs of the system.They are blocks in terms of which we will describe the (bi)linear relations underlying the definition of auxiliary functions.In other words, lowercase and uppercase auxiliary functions (see the definition in the first paragraph of Sect.2.3) are to be factorized solely in terms of these blocks.They are completely defined through the partial sum Eq.(3.3) over Young tableaux.
Before we continue, we set a convenient notation for the EAFs.Let {i 1 , i 2 , . . ., i a }, {j 1 , j 2 , . . ., j a } be the sets of indices corresponding to the initial and final tableaux, respectively, in the ordering explained in Eq.(3.1).i ℓ ≤ j ℓ for ℓ = 1, . . ., a. Then we set where the dependence on the spectral parameter has been omitted and the sum over k ℓ satisfies the standard filling rule.This is a good EAF for a Young tableau corresponding to the fundamental representation a if no factorization takes place, which happens for j ℓ ≥ i ℓ+1 with ℓ = 1, 2, . . ., n − 1.These are the column EAFs.It is not hard to see how this generalizes to other kinds of tableaux.

A canonical set of functions 4.1 Prescription from fusion
For the sl n invariant model in the a-th fundamental representation the T-system reads From Eq.( 4.1) the Y-system is defined as, see e.g.[30] thus Y (a,1) (x) = 1 + y (a,1) (x).As mentioned before we want to generate sets of auxiliary functions that are related to this very same expression.One way to achieve this is to modify the filling rules for the tableaux of Eq.(4.2) so that not all possibilities are allowed at once -as a consequence, instead of having solely the eigenvalues Λ (a,s) (x) in the definition of these auxiliary functions they will be constructed mainly with the EAFs introduced in the last section, which opens up a number of new combinations.
Let j 1 < . . .< j a , j k = 1, . . ., n be fixed numbers that act as boundaries for the indices in the boxes on the left hand side of Eq.(4.1).In order to fill up a column tableau from top to bottom we allow only a certain range of fillings determined according to its spectral parameter.More specifically, where we used the compact notation introduced in Eq.(3.6).As an example, take n ≥ 5 and a = 4 and apply these filling rules to the tableaux on the left hand side of Eq.(4.1).By means of the "fusion moves" explained in Sec.2.3 the rearranging of boxes yields Observe that the tableaux in Eq.(4.4) in general are not further reduced in the sense that there might be fillings of the individual boxes precluded by the filling rules of the Young tableau.For instance, if one takes j 1 = 1 then of course one cannot start the following box with the value 1.Instead, one should start with value 2, even though we indicate there that it may start from j 1 = 1.Therefore, reduction is circumstantial to the choice of the j k .On the other hand, the symmetrically fused tableau on the right hand side reduces completely regardless the choice of j k .Our construction transforms the sum in the symmetrically-fused tableau into a single term which can be factorized in terms of column tableaux.This truncation is important to avoid indefinite growth of EAFs (and therefore of auxiliary functions) thus making the NLIEs finite in number.In fact, for a sl n -symmetric model our approach provides 2 n − 2 = n−1 a=1 n a uppercase and lowercase auxiliary functions.This is also the number of EAFs that naturally appear in this construction.
Hence, as long as proper domains of the complex plane can be chosen for all EAFs, they can be eliminated in Fourier space thus providing a consistent set of NLIEs.

NLIE for sl 4 and sl 5
In this section we will apply the canonical construction in order to obtain auxiliary functions for sl 4 -and sl 5 -symmetric models.We separate the functions into sets corresponding to each fundamental representation of sl n , which is labeled by the index a = 1, . . ., n − 1.Each representation contains d a uppercase and lowercase auxiliary functions, which corresponds to its dimension.One can verify that all functions in a given representation a = k are conjugated to the ones in a = n − k.Also, each representation is invariant by species conjugation, i.e. the change in indices k ↔ n − k + 1, etc.It is worth noting that the auxiliary functions of sl 2 and sl 3 presented in Appendix B can be obtained through this approach.Later on, we derive sets of NLIEs involving these auxiliary functions.

sl 4
Using the canonical prescription, one finds 14 auxiliary functions for the sl 4 case.We observe that 10 of them coincide with the result of Damerau and Klümper [24].In the first representation, a = 1, one finds 1,3 (x) = 3 x−i/2 1,4 (x) = 4 x−i/2 1, 3 followed by the self-conjugated set of functions for a = 2, B  x−i/2 1, 2 2,6 (x) = 3 4 x−i/2 1, 2 and, finally, conjugated to the first representation, the functions corresponding to a = 3 are x−i/2 1 2, 3 3, 4 Now we want to replace the tableaux by the functions they represent.In the sl 4 case we have three eigenvalue expressions, namely although we should note that only Λ 1,1 (x) and Λ 3,1 (x) appear directly in our expressions.For convenience, in the calculations below we are going to consider the normalized quantities Also, we define functions to replace the subgraphs appearing in (4.6a)-(4.8d)as follows from Sec.3, 1,2 (x) q 1 (x)q 3 (x) , 2,1 (x) , . (4.13) By inspection, one finds that the zeros of all these EAFs are located close to certain lines of the complex plane.In the case µ j = 0, this is summarized in the table below.

EAF Im(x)
q j (x), j = 1, . . ., 3 0 1,j (x), j = 1, . . ., 5 ±1 2,4 (x) ±3/2 It is worth noting that the zeros of p 1,j (x) are known as the hole solutions to the Bethe ansatz equations λ j /λ j+1 = −1 for j = 1, 2, 3. We want to obtain nontrivial relations among {b (4) a,j (x)} and {B (4) a,j (x)}.To achieve this, the idea is to compute the logarithmic derivative of all functions and, in Fourier space, express the EAFs in terms of the {B (4) a,j (x)}.Inserting this into the expressions for the {b (4) a,j (x)} and taking the inverse Fourier transform yields the NLIEs.Thus we require these functions to be analytical, nonzero and with constant asymptotics (ANZ) in the strip − 1 2 < Im(x) < 1 2 .Taking into account the analytical structure of the EAFs according to the table above, we see that we must shift the argument of four functions as follows: where we relabeled these functions as denoted.Then, B 1,1 (x) = , (4.15b) , (4.15c) 3,1 (x) = 2,1 (x) e −βµ2 , (4.15m) e −βµ1 , ( while their lowercase counterparts are e β(−µ1+µ2+µ3+µ4) .(4.16n) It is important to notice that when taking the Fourier transform of the logarithmic derivative of the auxiliary functions one must treat the cases k < 0 and k > 0 separately.After eliminating the unwanted variables, we transform the resulting system back to real space and integrate with respect to x, which leads us to a set of NLIEs involving only b a,j (x) , given by log b (4) (x) = − c (4) + βJd (4) (x) − K (4)  * log B (4) (x), (4.17) where b (4) a,j (x) and B (4 a,j (x) , a = 1, . . ., 3, j = 1, . . ., d a .Here the logarithm is to be understood as applied to each function in the set.Also, c is the constant of integration c = c which is obtained from the knowledge of the auxiliary functions as x → ∞ together with the relation below for the convolution terms where F[f (x)] denotes the Fourier transform of the function f (x).The driving term d (4) (x) is where d (4,j) (x) = ∞ −∞ e ikx d(4,j) (k)dk, j = 1, 2, 3, which in turn are given by where v j are column vectors of dimension d (4) j whose entries are 1, and T j are diagonal matrices of this same dimension that encode the shifts of the auxiliary functions.As long as we do not violate ANZ strips, we can dispose of these shifts conveniently.For instance, for zero generalized chemical potentials we may assume that only T 2 is different from identity and reads diag(y −1 , 1, 1, 1, 1, y), with y = e k/2 .K (4) (x) is the kernel matrix, which can be written in terms of nine submatrices, They satisfy the following symmetries: where we have defined the reflection matrix [Π j ] mn = δ d (4) j +1−m,n .Therefore only four blocks of them need to be explicitly written.
In Fourier space, these submatrices and their entries are given by K( 4) where K( 4) K( 4) K( 4) K( 4)  4) K( 4) K( 4) K( 4) and finally  4) K( 4) K( 4) where the function K (a,b) [4] (k) appearing everywhere is defined by The NLIEs (4.17) allow to calculate the functions {B (4) a,j } by for instance numerical iterations.Now, computing the largest eigenvalue of the QTM Λ with d (4) (x) given by (4.20).It is noteworthy that the above results are much alike those of [24].We recall that the authors have found by trial and error the same number of auxiliary functions as one does by following our method and actually many of them coincide with ours.Furthermore the functions that are distinct are related to each other by means of the function f (4) (x), are given in Eqs.(B.3a)-(B.3e).Note that f (4) (x) is self-conjugated.When compared with the Y-system, we obtain Hence, the products of our uppercase functions within each representation do not give the Y-system per se but are related to them via f (4) (x) in a simple way.As a final remark, regarding the NLIEs we note that if one disregards the shifts (4.14b), whose effect is to introduce the factors y = e k/2 (or 1/y) our driving term and kernel matrices look the same and are given by the same functions as in [24] except for (4.32), (4.33) and (4.46) which differ by e −k , e k and 1, respectively.

Models and Numerics
For lower rank models, sl 2 and sl 3 , there are no differences between our proposed set of auxiliary functions and the successful approaches in the past [20].Nevertheless, with respect to the sl 4 we can already see that the canonical functions here proposed by us differ from [24].Both formulations allow for stable and fast algorithms, so in terms of numerical efficiency there is little difference.Using the same set of physical parameters (temperature and chemical potentials), as well as for the iterative solution of the NLIEs (discretization parameters in Fast Fourier Transform), we obtained the solution for the canonical functions after 42 iterations, whereas solutions for the old formulation were obtained after 38 iterations, Fig. 2. We see that the functions in question have the following advantages and disadvantages.The definition of the canonical set of functions is obtained in a systematic manner, however some functions take a slightly more involved trace accounting for slightly longer computation times in comparison to the functions used in [20,24] for which there is no systematic construction principle as to now.Nevertheless, this tradeoff reveals to be worthwhile, for now we may tackle the sl 5 case.Here we show that our approach works by displaying physical properties for some particular choices of physical parameters.When all chemical potentials are equal we have a balanced distribution n i = 1 5 , i = 1, . . ., 5. The physical significance of such an assertion depends on how we interpret the model.If we regard the permutation operator as the Hamiltonian describing exchange interaction in a spin chain with 2S + 1 = 5, then equal species densities result especially in zero magnetization.Nevertheless, here we propose the following Hamiltonian ) which describes hard core bosons of four types with nearest-neighbour hopping and interactions on the chain.The operators b α,i are bosonic in nature and may correspond to pairs of electrons that can be in triplet or singlet states, corresponding to the four possibilities of the indices α, β.Furthermore, one of the five local states plays the role of the vacant site.In this case, equal species density means that we have a fraction of holes of 1  5 which makes the density of the pairs to be 4 5 with no net magnetization.From the numerics, we can see that the system reaches the expected high and low temperature limits with respect to the entropy and specific heat, see Fig. 3.For high temperature all states become equally probable and we have S = log 5, whereas for the specific heat we find the Schottky anomaly, i.e. the specific heat approaches zero for low and high temperatures and thus has at least one temperature maximum.This is typical of lattice systems with a finite degree of freedom per site.
Although we keep the chemical potentials equal, we can evaluate derivatives of the thermodynamical potential with respect to them, leading to the compressibility χ i = ∂ni ∂µi T,µ\i and convertibility χ , where in the latter j ̸ = i and T, µ\i is the set of all chemical potentials except for µ i .Because of symmetry, we have χ i = 4χ i,j .In this sense, only one independent susceptibility exists from which the physical properties such as usual susceptibilities of the actual spin chain may be obtained.Moreover, at equal chemical potentials, we observe indeed that ∂ni ∂T µ = 0.In contrast, as we move to lower pair densities, there arises an asymmetry between hole and particle states.In view of this, two independent generalized susceptibilities exist, which we take to be χ 5 and χ α .
For concreteness, we let the local state No. 5 represent a vacancy, whereas α = 1, 2, 3, 4 represent local states populated by hard-core bosons in the corresponding particular configuration.Then χ 5 provides the charge susceptibility, while χ α gives specific susceptibility of particle species α.
Such a difference indicates a separation phenomenon between the charge degree of freedom and the remaining ones, like the spin and possibly others.This is also apparent when comparing specific heat profiles at high and low particle densities.One can observe that different types of elementary excitations take part in this model.Specially at pair density 0.1, the emergence of a second maximum seems to be imminent.Such structures were observed in the context of the t-J model, where Luttinger liquid properties were probed [35].In the present case, although our model is bosonic in character and it explicitly ignores pair breaking processes it still shows separation features.

Closing remarks
In our present understanding of the QTM approach to thermodynamics, NLIE have to be derived in order to work with the unusual discrete Bethe root patterns, namely discrete distributions with accumulation points, especially when considering the infinite Trotter limit.However, we bypass these problems by use of theorems for functions of complex variables applied to certain auxiliary functions with suitable analyticity properties.The actual set up of such auxiliary functions is a problem by itself and, so far, non-systematic procedures were used to define them.In this paper we have presented a systematic procedure.
The problem can be separated into different stages, for not only the auxiliary functions must have good analytical properties but also they must combine into sets that allow one to solve the resulting linear equations in Fourier space.Therefore the job is threefold: i) definition of auxiliary functions {B j = 1 + b j } with same number of EAFs that appear in their factorization, ii) analysis of their analytical properties, iii) expressing the EAFs in terms of the uppercase auxiliary functions, and in turn expressing the lowercase auxiliary functions in terms of the EAFs.
In this work we set out to generate functions i) systematically for sl n -invariant models as well as to derive the associated NLIEs.Our results are satisfactory in applications at least up to n = 5 which we demonstrated by explicit solutions to the NLIE.Our new route realizes a previous conjecture according to which the number of auxiliary functions should be 2 n − 2, while it abandons the heuristic guiding principle of factorization of the Y-system [36].Nevertheless, what is really essential is the factorization of Young tableaux with higher symmetric spin index.As a consequence, all EAFs are column-like tableaux and they come in finite number.There are differences between the previous approach and our actual proposition for NLIE in case of the sl 4 -invariant model, however the new approach in detail recovers the previously derived NLIE for the sl 2 and sl 3 -invariant cases [20].
We also elaborated on the determination of the EAFs.Using pole cancellation graphs, we extracted the fundamental column tableaux that do not factorize.This further stimulated us to introduce a compact notation for the EAFs.Naturally, many questions arise in this context.How many auxiliary functions are there that can be written in terms of column tableaux?Since there are many sets of equations for the same model, which one is the most suitable?Do they all lead to a complete set of NLIEs?Can we still find factorization for the Y-system?This addresses the question how to conform results with previous works.Factorization of the Y-system establishes the connection between the TBA and the QTM approach [37], and it would be desirable to carry out the program to higher-rank models.
Besides the lack of factorization of the Y-functions for n > 3, another difficulty of the new method in the elementary form presented here concerns the analyticity of the auxiliary functions.Some of these functions contain EAFs for which the distribution of zeros and poles requires a shift of the integration contour into the complex plane.It is not clear from the algebraic construction which ones or even how many functions require such a shift.In addition, these modifications introduce multiplicative factors that break some desired symmetries in the kernel matrix.We expect that working with functions based on Young tableaux with higher symmetric index s solves the "analyticity issue".
To sum up, we have presented a machinery to generate what we call "a canonical set of functions".However, as the questions posed by us indicate these are far from being unique or even completely understood.
x , B

15 K
(x) is relatively simple.Indeed, one of the expressions of the EAFs in terms of the {B

Figure 2 :
Figure 2: Traces of the curves of the canonical auxiliary functions that differ from the set of functions in [24].T = 0.1, µ = 0.

Figure 5 :
Figure 5: Entropy and specific-heat at pair density 0.1