Torus partition function of the six-vertex model from algebraic geometry

We develop an efficient method to compute the torus partition function of the six-vertex model exactly for finite lattice size. The method is based on the algebro-geometric approach to the resolution of Bethe ansatz equations initiated in a previous work, and on further ingredients introduced in the present paper. The latter include rational $Q$-system, primary decomposition, algebraic extension and Galois theory. Using this approach, we probe new structures in the solution space of the Bethe ansatz equations which enable us to boost the efficiency of the computation. As an application, we study the zeros of the partition function in a partial thermodynamic limit of $M \times N$ tori with $N \gg M$. We observe that for $N \to \infty$ the zeros accumulate on some curves and give a numerical method to generate the curves of accumulation points.

Computing partition functions of two-dimensional lattice models is one of the central problems in statistical mechanics. When the lattice model is integrable, one can often compute the partition function exactly. However, even for integrable models, exact results are typically only available in two limits, i.e., when the lattice size is very small or in the thermodynamic limit where the lattice size is infinite. For the intermediate case, obtaining exact results for the partition function is actually a hard task.
For definiteness, we consider the torus partition function of the six-vertex model at its isotropic point, with lattice size M × N . This is a well-known integrable model which is equivalent to the Heisenberg XXX spin chain [1]. The partition function can be computed by the transfer matrix method. Usual folklore of integrability tells us that we can diagonalize the transfer matrix by Bethe ansatz. The partition function can be written in terms of the eigenvalues of the transfer matrix. This works for any lattice size. However, the eigenvalues of the transfer matrix obtained in this way are only formal, since they are written in terms of Bethe roots, which are solutions of Bethe ansatz equations (BAE). To compute the partition function explicitly, we need to actually solve the BAE and find all physical solutions and then plug in the eigenvalues of the transfer matrix. What is usually not stressed in the literature is that it is in fact a highly non-trivial task to find all the solutions of the BAE, even numerically.
Morever, even if one finds all the solutions of BAE numerically, the result is not exact. Although numerical results are sufficient for many purposes, our goal here is to find exact results. In this work, we propose an efficient method to compute the partition function of integrable lattice model for finite-size system exactly and analytically without using any numerics. This method involves several recent developments in integrability, such as rational Q-systems and the algebro-geometric approach to Bethe ansatz equations. In this sense, the current work is a continuation of the work [2] initiated by two of us.
The general goal of the program started in [2] is to study the solution space of Bethe ansatz equations systematically by algebraic geometry and develop new analytical methods for physical applications in different contexts. In this paper, we extend the results of our previous work by introducing two new ingredients from algebraic geometry, namely primary decomposition and algebraic extension. Using these methods, we can probe structures in the solution space of the BAE. We find that the solution space can be decomposed naturally into non-intersecting subspaces upon primary decomposition on Q. The physical interpretation of such decomposition is related to the decomposition of the transfer matrix with respect to the total momentum. This decomposition can be performed more thoroughly The rest of this paper is organized as follows. We introduce the six-vertex model and its torus partition function in section 2. This part is standard and can be skipped by experienced readers. In section 3, we discuss the algebro-geometric approach to compute the torus partition function. We present explicit results for the partition function in section 4. For M ≤ 6, we have closed-form expressions for any N . For larger M , we give some partial results for fixed large N to convey an idea about the exact results. The full results can be found in the github repository. The partition function zeros in the partial thermodynamic limit are discussed in section 5, where we also compute the corresponding limiting curves. In section 6, we show how the algebro-geometric decomposition of the solution space of the BAE can be refined by working with an algebraic extension of Q. Using primary decomposition over this extension and some Galois theory provides us with a finer and computationally more efficient decomposition, which we can physically relate to the construction of momentum sectors for the transfer matrix. We state our conclusions and perspectives for further work in section 7. Appendices A-D contain details on more technical aspects, and the tables giving the number of physical solutions of the BAE for 6 ≤ M ≤ 18 are relegated to Appendix E.

Torus partition function of the six-vertex model
In this section, we review some basic facts about the six-vertex model, which also serve to fix our notations. The six-vertex model is a prototype of integrable lattice models. It is well known that it can be mapped to the Heisenberg XXZ spin chain and solved by Bethe ansatz. Throughout this work, we consider the isotropic point of the six-vertex model, which can be mapped to the Heisenberg XXX spin chain. We leave the more general case for future investigation.
The six-vertex model is a two-dimensional lattice model. At each site there are six possible configurations obeying the so-called ice rule, namely the number of incoming arrows should equal the number of outgoing arrows. The six possible configurations are depicted in figure 2.1. Each configuration can be represented in two ways. One is by putting arrows on each edge and the other is by using thin and thick lines. Each configuration is associated with an interaction energy ε i (z, θ), (i = 1, · · · , 6) subject to the following constraints ε 1 = ε 2 , ε 3 = ε 4 , ε 5 = ε 6 . (2.1) Following Baxter [1], we denote the corresponding Boltzmann weights by ω j = exp(−βε j ) and define a = ω 1 = ω 2 , b = ω 3 = ω 4 , c = ω 5 = ω 6 , (2.2) where in the isotropic case Indeed, the anisotropy parameter is then ∆ = a 2 +b 2 −c 2 2ab = 1, so the corresponding spin chain is the XXX one. We consider a lattice of M columns and N rows with periodic boundary condition in both directions, as is shown in figure 2.2. The partition function is given by summing the product of local Boltzmann weights over all possible configurations. The partition function of the six-vertex model can be computed by the transfer matrix method (see for example [1]). To define the transfer matrix, we start from the R-matrix R an (z, θ) = Integrability of the six-vertex model is guaranteed by the fact that the R-matrix satisfies the Yang-Baxter equation The transfer matrix is defined as R an (z, θ n ) . (2.6) The partition function can be written in terms of the transfer matrix where the parameters z j and θ k characterize the Boltzman weights at each site of the lattice. We consider the homogeneous case θ k = 0 and z j = z. Then the partition function is simply given by It is clear from the above definitions that it is a polynomial of degree M N in the variable z with rational coefficients.
Transfer matrix. Our task is to compute the trace of T M (z) N . The transfer matrix T M (z) is a matrix of dimension 2 M . The most straightforward way to compute the trace in (2.8) is by first constructing the matrix explicitly, performing the matrix multiplication N times and then taking the trace. To simplify this task, it is useful to notice that the transfer matrix is block-diagonal and we can construct the transfer matrix in each spin sector separately. The spin sectors are labelled by K where K = 0, 1, . . . , M is the number of vertical thick black lines in the right panel of figure 2.2. In the spin chain language, K corresponds to the number of flipped spins, or the number of magnons with respect to the pseudo-vacuum state | ↑ M . The dimension of spin sector K is We denote the transfer matrix of the spin sector labelled by K as T M,K (z). The partition function is then given by In what follows, we shall call this approach the brute-force method. It serves as a check for other approaches. This approach becomes cumbersome very quickly since the dimension d M,K grows rapidly with M and K. For example, d 12,6 = 924 which implies that for M = 12 we already need to deal with matrices of dimension about 10 3 .
Bethe ansatz. A better method which makes use of the integrability of the model is the Bethe ansatz. One can diagonalize the transfer matrix by directly constructing its eigenvectors using Bethe ansatz. The construction can be done in each spin sector. Let us denote the eigenvectors by |u , parameterized by a set of parameters u = {u 1 , . . . , u K } called rapidities. We then have Notice that for the spin sector K, the number of rapidities that characterize the state is K.
The eigenvalue t u (z) is given by and is called the Baxter polynomial. The rapidities are constrained by the Bethe ansatz equations (BAE) Bethe ansatz is a very powerful analytical method and it leads to the solution of the six-vertex model in the thermodynamic limit, when M, N → ∞. However, for fixed and finite M and N , the expression for t u (z) is only formal since the parameters {u 1 , . . . , u K } are not known explicitly. To find them, we need to solve the system of algebraic equations (2.15).
This raises some serious problems if we want to compute the partition function exactly and explicitly as a polynomial in z. First of all, the BAE for generic M and K are impossible to solve analytically and can only be solved numerically. Even finding the numerical solutions turns out to be a highly non-trivial task. Worse, not all the solutions of the the BAE lead to true eigenvectors and eigenvalues of the transfer matrix. We know that within each spin sector labelled by K, the dimension of the transfer matrix is d M,K . One therefore needs to show that there are exactly d M,K physical solutions of the BAE for fixed M and K. This is intimately related to the completeness problem of Bethe ansatz and is quite subtle.
Physical solutions of the BAE. In general, the number of solutions to the BAE is more than the number of physical states. It is a non-trivial problem to characterize physical solutions of BAE among all the solutions. Fortunately, this problem has been studied systematically in [10]. The conclusion is that the physical solutions can be classified as pairwise distinct non-singular solutions and singular physical solutions. For a detailed discussion of these solutions, we also refer to [2]. To single out these solutions, one needs to impose further constraints in addition to the original set of BAE [2,10].
To find the physical solutions, it is actually more convenient to work with other formulations of the BAE, in particular Baxter's T Q-relation and the rational Q-system. Baxter's Q-operator provides a powerful method for solving integrable lattice models [1]. The central equation of this method is an operator equation called the T Q-relation. In terms of the eigenvalues of the T and Q operators, the T Q relation becomes the following functional equation (2. 16) Note that this equation is basically equivalent to (2.12) if we multiply both sides of the latter by Q u (z). For a state of length M and magnon number K, t(z) and Q(z) are polynomials in the variable z of degree M and K, respectively. To solve the T Q-relation (2.16), one first writes t(z) and Q(z) in the explicit polynomial form The unknown variables that we are solving for are {t 0 , . . . , t M , s 0 , . . . , s K−1 }. Plugging (2.17) and the explicit form (2.13) of a(z), d(z) into the equation (2.16) and demanding that it is satisfied for any value of z, one obtains a system of algebraic equations for the unknown variables. These equations are linear in both sets of variables {t 0 , . . . , t M } and {s 0 , . . . , s K−1 } and are easier to solve than the original set of BAE. Solving the T Q-relation gives at the same time both polynomials t(z) and Q(z). The zeros of Q(z) in turn provide the solution of the BAE. Another advantage of using the T Q-relation is that it automatically eliminates the solutions with coinciding rapidities, i.e., the solutions of the T Q-relation lead to pairwise distinct solutions of the BAE.
Working with the T Q-relation instead of the original BAE is already more efficient for our purpose. However, the use of the T Q-relation does not eliminate all the unphysical solutions. Very recently, an even more efficient method has been developed for solving the BAE which is the rational Q-system approach [11,12]. In this method, one defines a rational Q-system associated with a Young tableaux with specific boundary conditions. The configuration of the Young tableaux is related to the length and magnon number of the Bethe state. By requiring all the Q-functions at the vertices of the Young tableaux to be polynomials, one ends up with a set of algebraic equations called zero-remainder conditions (ZRC) where the unknown coefficients are the coefficients of the Q-functions, i.e., the analogues of the {s 0 , . . . , s K−1 } defined above. Solving the Q-system amounts to finding all the Q-polynomials on the Young tableaux. One specific Q-polynomial coincides with Q(z) defined above, so its zeros give the Bethe roots. For more details of this approach and some explicit examples, we refer to appendix A.
It turns out that solving the ZRC is even more efficient than solving the T Q-relation. Moreover, the solutions of the ZRC are in one-to-one correspondence with the physical solutions of the BAE, so no further constraints need to be imposed. One minor disadvantage of this method is that the equations themselves are not known explicitly for any length and magnon numbers and have to be derived case by case. This is not a big issue in practice since the equations can be generated rather efficiently for not too large M and K.
In what follows, to find all the physical solutions of the BAE, we will work with rational Q-systems. So far, the Q-system approach for BAE has only been developed for the isotropic limit (∆ = 1) corresponding to the XXX spin chain. To treat the XXZ spin chain at generic ∆, we still need to rely on the T Q-relation, together with additional constraints to select physical solutions. The strategy is to first find the Q-polynomial, and next use the T Qrelation to find the t-polynomials. After finding all the t-polynomials, it is straightforward to write down the torus partition function that we are after.
However, it is clear that if we want to find all the t-polynomials explicitly, no matter which approach (BAE, T Q-relation, Q-system) we are using, we have to solve the system of algebraic equations numerically. The results are thus approximate and not exact. To avoid solving equations and obtain instead exact and analytical results for the partition function, we can however apply methods in computational algebraic geometry. We discuss these approaches in the next section.
3 Algebro-geometric approach to the partition function In this section, we describe our method for computing the torus partition function using an algebro-geometric approach. The main tools that we are going to use are Gröbner basis, quotient ring and companion matrix. See [13] for a textbook reference to the corresponding mathematics. For a detailed introduction to these notions in the context of Bethe ansatz, we refer to [2].
As discussed in the previous section, in order to select all the physical solutions, we work with the rational Q-system. For the su(2) invariant XXX spin chain which is equivalent to the six-vertex model, the corresponding Young tableaux have two rows with the number of boxes being (M −K, K). The Q-polynomial that we are interested is Q 0,1 . The computation of the Q-polynomials relies on the definition of certain paths on the Young tableau (for more details, see appendix A). For the su(2) case, we can choose the path such that the unknown coefficients are precisely the coefficients of Q 0,1 (z), and we have Ideal. The zero-remainder conditions (ZRC) then give a set of algebraic equations for the K variables {s 0 , . . . , s K−1 }, where f k (s 0 , s 1 , · · · , s K−1 ) are polynomials in the variables {s 0 , · · · , s K−1 }. Here S is the number of equations and it depends on the path we choose. The polynomials f 1 , · · · , f S define an ideal in the polynomial ring C[s 0 , s 1 , · · · , s K−1 ], denoted A given ideal can be generated by different bases, among which the so-called Gröbner basis is particularly useful for us. We denote the Gröbner basis by G k where in general S and S are different. Notice that when computing the Gröbner basis we need to choose a partial ordering of the monomials formed of the variables {s 0 , s 1 , · · · , s K−1 }. For different orderings, the corresponding Gröbner basis can look quite different.
Quotient ring. The quotient ring is defined as and it is a finite-dimensional linear space. The dimension of this linear space equals the number of physical solutions of the BAE for given M and K, which is given by [10] Since Q M,K is a linear space, it can be spanned by a basis set. The standard basis of the quotient ring Q M,K is given by all the monomials of {s 0 , · · · , s K−1 } that cannot be divided by LT[G k ] (k = 1, · · · , S ) where "LT" stands for the leading term in some given partial ordering. In order to construct the standard basis, one needs to calculate the Gröbner basis of the ideal f 1 , · · · , f S . This is one of the main calculations of the current work which can be done by standard algorithms such as Buchberger's algorithm or the F4/F5 algorithm of Faugère [14,15]. These algorithms are implemented in several packages for algebraic geometry, such as Singular [16] .
Companion matrix. Let us denote the basis of the quotient ring by {e 1 , e 2 , · · · , e N M,K }. Any polynomial P (s 0 , s 1 , . . . , s K−1 ) can be mapped to a numerical matrix M P called the companion matrix of dimension N M,K . The algorithm for doing so is as follows. We multiply the polynomial P by one of the standard basis elements e j and then find the remainder of the polynomial reduction with respect to the Gröbner basis, The remainder r j (s 0 , · · · , s K−1 ) can be expanded in terms of the standard basis where the coefficients of the expansion are the elements of the companion matrix, i.e., The companion matrix satisfies the following homomorphism properties The main result from algebraic geometry which we are going to use is sol P (s 0 , · · · , s K−1 ) = Tr M P , (3.10) where the sum 'sol' is over all solutions of the system of equations (3.2). It is straightforward to see that if we start with equations whose coefficients are in Q, the right-hand side of (3.10) is rational (viz., it is a polynomial in z with rational coefficients), even though individual terms appearing on the left-hand side may be irrational.
Using the procedure described above, after the construction of quotient ring, we can map the function Q 0,1 (z) into a companion matrix by mapping each coefficient s k → S k , so that where I is the identity matrix of dimension N M,K . Our goal is to find the companion matrix of the transfer matrix t u (z), since this will permit us to access the partition function. This can be done by using the explicit expression of the transfer matrix (2.12) and the properties of the companion matrix (3.9). More explicitly, we have where T M,K (z) can be computed using the ideal generated by the T Q-relations 1 or from the companion matrix Q M,K (z) by the following relation Several comments are in order. The multiplicities M − 2K + 1, for K ≤ M/2, take into account the descendant states in the Bethe ansatz. These states are obtained by sending some of the Bethe roots to infinity. It is easy to see that adding a root at infinity does not change the eigenvalue of the transfer matrix. When solving the BAE or the ZRC, we only find regular solutions that do not have roots at infinity. Since the descendant states are indeed part of the spectrum, we need to take them into account in the computation of the partition function by the proper multiplicity.
The expression in (3.14) takes a very similar form to the one in (2.10). However, there are important differences. Firstly, the result in (3.14) makes use of the full su(2) symmetry, and hence the dimensions of the transfer matrices T M,K (z) are smaller than those of T M,K (z). The dimensions of T M,K (z) and T M,K (z) are N K,M and d K,M respectively, see (3.6) and (2.9). For example, for M = 14, K = 7, we have N 14,7 = 429 and d 14,7 = 3432. Therefore (3.14) is computationally more efficient than (2.10).
Secondly, one can impose further constraints on the solution space of the BAE and decompose the solution space into even smaller subspaces for (3.14). This means that we can make the matrices T M,K (z) in (3.14) block-diagonal and hence work with even smaller matrices, which improves the efficiency further. This point will be discussed in detail in section 6.

Explicit results
In this section, we give results of partition function for different values of M and N . We obtain closed-form results up to M = 6 for arbitrary N . The reason we stop at M = 6 is that the ZRC can be solved analytically up to this length in Q. For M = 7 and M = 8, analytical solutions of the ZRC can also be found by working in extended fields; see section 6 for more details. For M ≥ 9, we give an efficient algorithm for computing the partition function for fixed M and N (which can be large). The results are polynomials in z of high degrees with rational coefficients (typically large) and it does not make sense to write them down in this paper. Instead, interested readers can find the results on the repository which we mentioned in the introduction.

Closed-form expressions for M ≤ 6
In this section, we list the results for M ≤ 6. We denote the partition function of an M × N toroidal lattice by Z M,N . M = 1. This is the simplest case. The sum in (3.14) contains only K = 0, and the result is given by (4.1) M = 2. We have to sum over K = 0 and K = 1. The partition function is given by We again have to sum over K = 0 and K = 1. The result is We see that irrational numbers start to show up within some eigenvalues of the transfer matrix. However, for any N ∈ N the result for Z 3,N is a rational-coefficient polynomial in z after simplification. M = 4. We have to sum over K = 0, 1, 2. The result is We have to sum over K = 0, 1, 2. The final result reads We see that the eigenvalues of the transfer matrix now become more complicated, with double square roots showing up in the coefficients. M = 6. We have to sum over K = 0, 1, 2, 3. The final result reads Again, we see that some of the terms in (4.5)-(4.6) contain multiple square roots. However, once we sum over all terms for any N ∈ N, we obtain polynomials whose coefficients are rational numbers, as expected.

Partition functions for higher M and N
For M = 7 and M = 8, we can also work out the analytical results. However, the closedform results involve complicated multi-square roots and are not very illuminating to write down explicitly here. For M ≥ 9, we are not able to find analytical solutions anymore. For these cases, we give an efficient approach to compute the partition function for fixed M and N based on the companion matrix.
Although our method is much more efficient than the brute-force approach, the complexity still grows exponentially with M . On a laptop, we are able to compute the Gröbner basis and companion matrices of the Q-polynomial Q(z) and the transfer matrix Approximately this number is c 50 ≈ 6.75536 × 10 125 . Using numerical methods such as the function Rationalize in Mathematica to guess such a large number will be quite difficult in practice since it requires working with floating point numbers with very high accuracy.
Since the results for partition functions for large N are typically large, we find it more useful to give the companion matrices. These companion matrices contain all the information we need. To find the explicit eigenvalues of Q u (z) and t u (z), we can diagonalize the companion matrices. In general this diagonalization can only be done numerically, but the matrices that we need to diagonalize are much smaller and can be handled much more efficiently. The zeros of Q u (z) give the Bethe roots. That is to say, we can straightforwardly find all physical solutions of the BAE up to length M = 14 using our results. The eigenvalues t u (z) contain all the information about conserved charges of the system, i.e., momentum, energy and higher conserved charges of the Bethe states.
To find the exact partition function Z M,N , we need to take matrix powers of the companion matrices T M,K (z) N and then take the trace. We will be interested in the case of large N . Naively, we would need to perform N matrix multiplications involving the companion matrix before taking the trace. When the size of the matrix is large, the analytic computation of matrix multiplication become time consuming. To reach a high value of N , we actually need a better way to do the computation. This is described in appendix C.

Consistency check
Since our results are usually large polynomials, it is important to perform some checks for their correctness. One important consistency check is the following. We can compute the lattice partition function in two ways, corresponding to two different choices of the transfer direction, and the result should be the same, as is shown in figure 4.3. Specifically, we consider the transformation of the partition function where we rotate the lattice by π/2. To analyze the effect of this transformation, we consider the M × N lattice with M vertical lines and N horizontal lines. We associate to each vertical line the same spectral parameter θ, and to each horizontal line the same spectral parameter z. We denote this partition function by Z M,N (z − θ). There are six possible configurations at each site with three different Boltzmann weights given by The partition function can be written as where f m,n,k denotes the multiplicity of configurations with m, n, k vertices with Boltzmann weight a, b, c respectively. Rotating the lattice by π/2, it is clear, from the first line of figure 2.1, that we have the following transformation Therefore, under this transformation, which we denote by R, the partition function transforms as The number of type-c vertices k is even due to the ice rule. Therefore, if M N is even, the transformation leaves the partition function invariant, otherwise it gives an additional minus sign. This invariance is a strong consistency check.
We use this relation to check the correctness of the companion matrices as follows. For a given length M , we construct the corresponding companion matrices and compute the partition function Z M,N (z − θ) for all N ≤ M . These partition functions can also be computed as Z N,M (θ − z) using the companion matrices constructed for length N . If the results are correct, the two calculations should give the same result. We have checked the companion matrix in this way up to M = 14.

Zeros of partition functions
The torus partition function Z M,N (z) is a polynomial of order M N in z. It is instructive to find the zeros of this partition function in the complex z-plane.
Partition function zeros for statistical models with one complex parameter have been studied in a variety of contexts-including Lee-Yang zeros (complex magnetic field) [17], Fisher zeros (complex temperature) [18], and graph polynomials such as the Q-colour chromatic polynomial [3][4][5][6][7][8][9]-and has given rise to an immense literature (see, e.g., references in [3]). The zeros of partition functions in the thermodynamic limit contain information about the phases and critical behavior of the model at hand. In many cases the zeros will accumulate on curves, for M, N → ∞, which "pinch" the real axis at one or more critical points. Isolated accumulation points provide another possible scenario. In the simplest cases-such as the Ising model with suitable boundary conditions-the curves of accumulation points can be proved to form circles, giving rise to so-called circle theorems. More generally, the density of zeros near a critical point obey scaling laws that can be related to critical exponents.
The true thermodynamic limit of an M ×N system is obtained by letting M, N → ∞ with a fixed and finite ratio, 0 < M/N < ∞. But another means of obtaining relevant information is to fix a finite value of M , and study the accumulation points of zeros as N → ∞. For a partition function of the form (3.14), and supposing a mild non-degeneracy condition, the Beraha-Kahane-Weiss theorem [19] states that the accumulation points will form curves. By standard analyticity theorems, any closed region delimited by such curves will constitute a thermodynamical phase (for N → ∞). Under reasonable (but not entirely innocuous) assumptions about the commutativity of limits, the phase diagram in the thermodynamic limit can then be inferred by studying the convergence of these curves upon taking M → ∞.

Partition functions for different M and N
In this subsection, we give the partition function zeros for different M and N . We fix the value of M and increase N to see how the distribution of the zeros change. In this way, we try to extrapolate the behaviors of the zeros to the (partial) thermodynamic limit N → ∞ (with M being fixed and finite). One might wonder what is the benefit of knowing the partition function exactly for finding the zeros. Naively one might expect that numerical approximations will be sufficient for finding the zeros of the partition function. However, it is known that the locations of zeros of a polynomial can be very sensitive to perturbations of coefficients, especially when the degree of the polynomial is large. One famous example is the so-called Wilkinson's polynomial where a change of one of the coefficients by 10 −7 leads to significant changes in the locations of the zeros. Having exact results eliminates this potential subtlety. These are polynomials of degree M N , and the coefficients can be normalized to be integers by multiplying Z M,N (z) by the overall factor 2 (M −1)N . Given the very large degree and the size of their coefficients, it is actually a non-trivial task to compute the zeros of these polynomials. These difficulties are however efficiently overcome by the application of the software MPSolve [20], which is a multiprecision implementation of the Aberth method [21]. The main advantage of the latter method is that it approximates all the roots of a univariate polynomial simultaneously.
The zero plots in figure 5.4 reveal several interesting features. As the aspect ratio ρ grows, the zeros tend to settle on certain curves-the limiting curves of accumulation points to be discussed further in the next subsection. In the regions close to the origin the finite-ρ effects are small, but further away their importance increases, and the fine structure of the limiting curves are barely visible even at the largest ρ shown. Moreover, far from the origin the density of zeros is very scarce. Regarding the case M = 14, it seems likely that it would develop rich details as those seen in the other plots, provided large ρ could be accessed. In particular, there are "stray" zeros around the central almost-horizontal branches that appear as precursors of multiple branches and T-points. While all these features could certainly be analyzed at length, we instead move on to the direct determination of the limiting curves as ρ → ∞.

Limiting curves
In the previous subsection, we have seen that as N increases, the zeros of the partition function accumulate on some curves. Following [3] we shall refer to these as limiting curves.
By the Beraha-Kahane-Weiss (BKW) theorem [19], this is a consequence of the form (2.10), or equivalently of (3.14), that relates the partition function Z M,N to sum over traces of the N 'th power of the transfer matrix T M,K (z), or of the corresponding companion matrix T M,K (z) given by (3.12).
More precisely, the BKW theorem applies to an expression of the form where we shall refer to the Λ i (z) as eigenvalues, and the α i (z) as the corresponding multiplicities. For a given z, let us order the eigenvalues by norm, so that |Λ 1 (z)| ≥ |Λ 2 (z)| ≥ · · · , and we call an Supposing a mild non-degeneracy condition, the BKW theorem then states that the accumulation set of zeros, as N → ∞, will form either isolated points or curves. An isolated accumulation point occurs for z = z 0 , when there is a unique dominant eigenvalue (i.e., |Λ 1 (z 0 )| > |Λ 2 (z 0 )|) and the corresponding multiplicity vanishes (i.e., α 1 (z 0 ) = 0). A curve of accumulation points occurs when there are at least two dominant eigenvalues (i.e., |Λ 1 (z)| = |Λ 2 (z)|), and the relative phase φ(z) ∈ R defined by Λ 2 (z) = e iφ(z) Λ 1 (z) varies along the curve. The speed of variation of φ(z) along the curves can be related to the density of partition function zeros [3]. Note also that the limiting curves may have T-points or higher-order bifurcations at a point z 0 where more than two eigenvalues are equimodular. We refer to [3] for more details on the BKW theorem and the detailed analysis of the generic setup.
In our context, α i (z) and Λ i (z) depend on M and K, and moreover α i (z) = M − 2K + 1 are simply constants. Therefore all accumulation points form curves, and not isolated points, in agreement with the observations of the preceding subsection. To trace these curves for a given M , we use an approach for identifying the loci of equimodularity that is described in appendix D. This consists in two steps: first we identify some points of equimodularity by a direct search (e.g., along suitably chosen straight lines), and second we trace the equimodular curves starting from each of those points, using a procedure explained in the appendix. While this approach may fail to detect very small curves of accumulation points, we believe to have obtained complete results for M ≤ 8.   A close scrutiny of figure 5.5 reveals that the limiting curves some tantalizingly tiny features for M ≥ 6, including almost-parallel curves and short stems linking the various branches. We have taken great case to represent (what we believe to be) all such features.
From a numerical point of view, the diagonalization computations become increasingly difficult as we approach the highly degenerate points z = ±i/2. Practical details about the computational approach to limiting curves can be found in appendix D.
Comparing figures 5.4-5.5 gives convincing evidence that the partition function zeros indeed accumulate on the limiting curves, in the limit of large aspect ratio ρ → ∞. That this is indeed the case is proved by the BKW theorem. However, it is also clear that some parts of the limiting curves are very scarcely populated by the zeros, even for the large values of ρ shown in figure 5.4. Moreover, some of the fine details of the limiting curves are hardly discernable on the plots of zeros, such as the short stem-like pieces connecting the almost-parallel branches for M = 8 or the (barely visible) sliver-shaped enclosed region for M = 6. Figure 5.7 shows a comparison between the limiting curves and the M = 7 partition function zeros for various aspect ratios. We have examined the curves delimiting the region close to the origin in some more detail. For even M , it has a corridor-like aspect, with the branches containing z = ±i/2 appearing to become more horizontal as M increases. For odd M , it forms an elongated bubble, with the above-mentioned vertical ray in the middle, whose upper and lower boundaries tend as well to become more horizontal as M increases. To illustrate the size-dependence of this bubble, we have produced partial results for this part of the curves for higher, odd values of M (up to M = 17), as shown in figure 5.6. The first intersection with the real axis appears to extrapolate to z = 2 as M → ∞ (see the right panel of figure 5.6). Based on this, we conjecture that the enclosed region tends to a rectangle, given by 0 < Re z < 2 and −1/2 < Im z < 1/2, in the thermodynamic limit, for M odd.

Primary decomposition
Let us summarize what we have achieved so far in computing the torus partition function of the six-vertex model. Computing the partition function by brute force, we need to work with matrices of dimension d M,K within the spin sector of K magnons. Using Bethe ansatz and the algebro-geometric method, we are able to reduce the problem to the computation of companion matrices of dimension N M,d = d M,K − d M,K−1 . This reduction in the dimensions of the matrices makes use of the full su(2) symmetry of the theory. Recall that we classify the Bethe states as primary states and their descendants with respect to the su(2) algebra. Since all the descendant states of a given primary state have the same eigenvalue of the transfer matrix, we can focus on the primary states only. The number of primary states is much less than the total number of states in a given spin sector.
Can we do better? Notice that we have not yet exploited all the symmetries of the model. For example, the model is also invariant under a lattice translation. This symmetry leads to the the total momentum of the Bethe state being quantized, taking only a finite number of possible values. Therefore, apart from decomposing the Hilbert space according to spin sectors, we can also decompose the Hilbert space according to momentum sectors, i.e., states with different values of the lattice momentum. These two decompositions can be performed simultaneously and leads to even smaller companion matrices. This will greatly enhance the efficiency of our approach.
Mathematically, the decomposition with respect to momentum sectors is intimately related to primary decomposition and algebraic extension in algebraic geometry. Physically, this decomposition also allows us to probe much deeper into the solution space of BAE and find new structures that have not been studied in the literature. In this section, we discuss the decomposition of the solution space with respect to momentum sectors. We first introduce the notion of primary decomposition on Q from some interesting observations about the partition function. We will see that it is useful to perform the decomposition on a larger field Q(i, ξ M ) obtained by an algebraic extension. In addition, exploiting Galois theory in the current context, we will show that many of the subspaces after the decomposition are actually related by the Galois group, and it is thus sufficient to perform the computation for a representative. The decomposition together with Galois theory lead to a huge boost in the efficiency of our computation. More details are given in appendix B and an upcoming publication [22].

Primary decomposition over Q
To see that the solution space of the BAE has more structure, we take a careful look at the closed-form results of the partition function in (4.6). We can see that it is natural to group some of the terms together since they take very similar forms. For example, we can group the following four eigenvalues in (4.6): , (6.1) , , One can check that although each Λ 1 , · · · , Λ 4 is complicated and has irrational coefficients for generic z, their symmetric power sums are always polynomials whose coefficients are rational numbers! Since each Λ i corresponds to a solution of the BAE or the T Q-relation, this implies that we can group the four corresponding solutions of the BAE. Notice that we cannot make the decomposition further on Q. If we further divide the four solutions into two groups, say Λ 1 , Λ 2 and Λ 3 , Λ 4 , then the coefficients of the symmetric power sums Λ n 1 + Λ n 2 and Λ n 3 + Λ n 4 are no longer rational. This implies that these four solutions form an irreducible or primary block on Q. Similarly, the remaining terms in (4.6) can be divided into such primary blocks. In geometrical terms, this grouping is equivalent to decomposing an affine variety into independent components. Such an operation is called primary decomposition in algebraic geometry. We refer to appendix B for more details.
Given an ideal, it is straightforward to compute the primary decomposition using standard algorithms. To understand the physical meaning of primary decomposition, we now analyze the example M = 6 carefully.
An example: M = 6. The result of the primary decomposition is given in table 1. Let us consider the spin sector K = 2. From table 1 we see that there are 9 physical solutions for M = 6, K = 2, and that these solutions can be divided into four groups, with dimensions 1,2,2,4. In particular, the dimension-4 subspace corresponds to the four eigenvalues given in (6.1). K M = 6 1 5=1+2+2 2 9=1+2+2+4 3 5=1+2+2 Table 1: Primary decomposition of solution space of BAE with M = 6. The numbers on the right-hand sides of each line represent the dimension of each subspace.
As we alluded to before, this decomposition is related to the lattice translational invariance which is generated by the shift operator U = e iP . As an operator, it is related to the transfer matrix as  not the same within each subspace. However, if we compute the eigenvalues of the operator U 3 , we find that they are the same within each subspace, as is shown in the last column of table 2. This is due to the fact that we work on the field Q. Let us denote ξ M = exp(2πi/M ). It is clear that ξ M is not always rational for all = 1, · · · , M . Therefore one cannot perform the decomposition over momentum sector completely on Q. For each M , we can find the smallest integer 1 ≤ m ≤ M such that the all the eigenvalues of U m are rational. Then we can perform the decomposition with respect to the eigenvalues of U m . As a result, we can restrict ourselves to each subspace by imposing an additional constraint on the original BAE. In our example, for M = 6, K = 2, the additional constraints for the four subspaces are Notice that we need to include the constraints U = ±1 in the cases C and D because otherwise they will include cases A and B.
Algebraic extension. As we see in the previous discussion, the primary decomposition is related to the decomposition with respect to the lattice momentum. Due to the fact that ξ M is not always a rational number, we cannot perform the decomposition completely. However, we are not constrained to work on the field Q. If we extend the field slightly to include ξ M and perform the primary decomposition on the extended field, then the decomposition with respect to the lattice momentum can be performed completely. More precisely, the extended field will turn out to be F M = Q(i, ξ M ) where i is the imaginary unit.
After the decomposition into momentum sectors, we have M subspaces 2 (corresponding to = 1, · · · , M ) in each spin sector K. In principle, we need to compute the Gröbner basis and companion matrices of all subsystems. However, we will show that by making use of the Galois group of the algebraic extension, we just need to calculate a very few BAE subsystems. We get the contribution from all subsystems by the Galois group actions.

Primary decomposition over F M
We explain in detail how to implement the decomposition over F M in practice. Since we work with the rational Q-system, it is most convenient to express the momentum condition in terms of Baxter polynomials Q(z) as Here s = {s 0 , · · · , s K−1 } are the unknown coefficients that we solve for in the rational Q-system. Alternatively, we can write Q s (z) as where {u 1 , · · · , u j } are the Bethe roots. When the solution of BAE is singular, i.e., two of the Bethe roots are ±i/2, we have Q s (±i/2) = 0, whence the left-hand side of (6.6) is singular. This singularity can be eliminated by using the T Q-relations, as we will comment on below.
Let I M,K be the ideal of the Q-system, for a spin-chain state of length M and magnon number K, in the variables s. When Q s (±i/2) is non-singular, we can write the momentum condition (6.6) in the following polynomial form From the construction of these ideals, we see that for any point x ∈ Z(I M,K, ), and it is clear that for x ∈ Z(I M,K,∞ ), This is the ideal decomposition which is crucial for the efficient computation of exact partition function via Gröbner bases. Note that in this paper, we do not prove that for ∈ {1, . . . , M, ∞}, each I M,K, is primary over the field F M , i.e., that there exists no further decomposition beyond the computation in this paper. This discussion is left for future work. Another comment is that one may well expect that in addition to the lattice translation symmetry, there can be other discrete symmetries such as reflection symmetry that may play a similar role. Namely, we can further decompose the solution space with respect to these symmetries. This interesting possibility is also left for future work.
With the decomposition (6.11), the exact partition function is presented as a sum over the contributions from the ideals in (6.14), in the ideal I M,K, . Note that T M,K, (z) contains polynomials in i and ξ M but no other algebraic numbers. As we will see, each T M,K,l (z) has a much smaller size than the original companion matrix T M,K (z), hence (6.15) provides a highly efficient way of computing the partition function.
Finally, we comment on the Bethe roots in I M,K,∞ , i.e., singular roots. The singularity of the eigenvalue of U in terms of Q s (±i/2) is actually spurious. They can be eliminated by using the T Q-relation. We can combine the equations from the rational Q-system and the T Q-relation and then eliminate the variables s. The elimination procedure is actually quite simple, due to the structure of equations from T Q-relations. This gives us a set of equations that only involve the variables t = {t 0 , t 1 , · · · , t M }. The momentum conditions in terms of t variables are simply (−i) M t t (i/2) = ξ M , = 1, 2, · · · , M (6. 17) and are free of singularities. We can of course directly work with the equations involving only t variables and there will not be the spurious class I M,K,∞ .
The reason that we also work with equations involving the s variables is that we can separate the regular and singular solutions in this case. Singular physical solutions are special among the solutions of the BAE, since naively plugging them into the eigenvalues and eigenstates lead to divergences and one needs to perform judicious regularizations [23]. According to the conjecture of [10], all the physical solutions of BAE can be divided into regular and singular physical solutions. These authors worked out the number of these two kinds of solutions up to M = 14, K = 7 and checked the validity of the conjecture. As a byproduct of the current paper, we can actually provide more data points up to M = 18, K = 9 and find that the conjecture still holds up to these values.
Galois theory. The decomposition (6.15) is a very convenient expression. Moreover, there is a further short-cut for the computation. The equations in different I M,K, are related by Galois group actions; therefore, instead of exhaustively going through the sum over all 's in (6.15), we just need to compute a few 's, namely one for each orbit.
Note that for two different decomposed BAE with 1 ≤ 1 , 2 ≤ M described in the previous subsection, if we replace ξ 1 M by ξ 2 M in the generators from I M,K, 1 (6.9), without changing the imaginary unit i or any rational coefficient, then the ideal I M,K, 2 is obtained. This implies that we need to consider the field automorphisms of F M which keeps i and the rational numbers invariant, or the Galois group G ≡ Gal(F M /Q(i)) where Specifically, if there exists an element g ∈ G such that, (6.19) then by the field automorphism of Gröbner basis computation, Hence it is important to analyze the structure of the algebraic extension [F M : Q(i)]. We consider three different cases for M . The analysis is based on elementary Galois theory. Here we just list the classification results, and the proof will be presented in the future work.
1. M is odd. In this case, the field F M is the cyclotomic field, (6.21) Note that the denominator of the reduced fraction 2 / 1 is relatively prime to M/ gcd(M, 1 ), so the congruence condition for 2 / 1 is meaningful.
The condition (6.27) for 4|M is complicated. However, it is possible to simplify it and get a similar condition as in the previous two cases. We notice that there is an enhanced symmetry for the Q-system, Two sub-systems, whose -values live in the same class, are equivalent by a Galois group action. Hence the number of solutions to the two sub-systems must be equal.
We compute the Gröbner basis of I 6,K, with = 1, 2, 3, 4, 5, 6, ∞, and get the Bethe root counting for the decomposed BAE with M = 6 in Table 3. For the singular Bethe roots, via the T Q-relation, we can find the values of their regularized momenta. These regularized values are indicated by the numbers between brackets. For example, in Table 3, the entry "0 (1)" for K = 2 and = 3 means that I 6,2,3 has no solution but there is one singular Bethe root whose regularized momentum value is e 3×2πi/6 = e πi . Note that since 7 is a prime number, the classification of BAE sub-systems is very simple. The root counting for the decomposed BAE with M = 7 is given in Table 4.  Bethe root can be expressed as radicals of rational numbers for M ≤ 8. The closed-form expressions for the partition function with M = 7, 8 will be presented in the future work.
Similar results for 9 ≤ M ≤ 18 are presented in appendix E (see Tables 6-15). We use the Gröbner basis method with coefficients in the algebraic extension F M to get the Bethe roots classification. Our code is powered by Singular [16].

Relation to naive momentum-sector diagonalization
We can relate the above results for the counting of Bethe roots for the decomposed BAE to a more naive diagonalization of the transfer matrix in sectors with specified magnon number K and momentum . To this end, we start from an example to parallel the discussion in section 6.1.
An example: M = 6. We first classify the states of the six-vertex model transfer matrix T M (z) given by (2.6). We work in the particle picture corresponding to the right panel of Figure 2.2 and we let 1 (resp. 0) denote the presence (resp. the absence) of a particle at a given lattice site. To access the momentum information, we pick a single representative state for each orbit of the cyclic group C M , and we denote its orbit length by g K .
By particle conservation, the full transfer matrix T M (z), of dimension 2 M , is a direct sum of blocks T M,K (z), each having dimension d M,K = M K . One may now further blockdiagonalize the T M,K (z) into momentum sectors T M,K, (z), having the dimensions d M,K, given by the above table, by using a procedure similar to Appendix A.4 of [24]. To this end we write T M,K, (z) = S out T M,K (z)S in . (6.37) Here S in is a d M,K, × d M,K matrix that maps each compatible orbit into its representative state, with weight g s . And S out is a d M,K × d M,K matrix that maps each state into a representative (and hence into an orbit), and attributes a weight (ξ M ) ·k /g s if a state from orbit s (not necessarily its representative) has to be shifted cyclically through k lattice steps (say, towards the right) in order to make it coincide with the representative state of s. We recall that ξ M = exp(2πi/M ), as before. One may now verify by explicit construction of the matrices T M,K, (z) that the spectrum of T M,K (z) is indeed the union of the spectra of the momentum blocks T M,K, (z).
As it stands, this method does not yet take into account the su(2) symmetry of the XXX spin chain. This means that each T M,K, (z) contains all the highest-weight states with equal or higher spin (i.e., K ≤ K) in its spectrum. To correct this, on the level on the counting, it suffices to subtract from each row in the table the one just above it, and we arrive at: This can finally be compared with Table 3. It is seen that the two tables are identical, in so far as they assign the same dimensions to the same (K, ) sectors. Notice that the present approach does not particularize the singular case denoted {∞} in Table 3, but assigns to it straight away the correct regularized momentum, namely {∞} → {3} for K = 2, and {∞} → {6} for K = 3. 3 The counting is consistent with the sum of numbers outside and inside the brackets in Table 3.
General case. The case of general M can be treated in the same way. To recover the results corresponding to Tables 3-15, we only need to know the number of compatible orbits under C M for each set of (K, ). In even simpler terms, pick a divisor g s |M , and let N (g s ; M, K) be the number of orbits of length g s with K magnons. For instance, we have N (g s ; 6, 2) = 0, 0, 1, 2 and N (g s ; 6, 3) = 0, 1, 0, 3 for g s = 1, 2, 3, 6, respectively. From this data, the dimensions d M,K, is the sum over those N (g s ; M, K) that respect the criterion (6.36), and the highest-weight combinations d M,K, − d M,K−1, provide precisely the numbers of Tables 3-15, up to the assignment of a definite momentum to the {∞} classes.
We have written a simple algorithm that carries out this computation. It produces the tables for M ≤ 24 in less than one minute. For M ≤ 18 these are in full agreement with Tables 3-15, after assigning to each {∞} case the corresponding regularized momentum. We note that for M even, this assignment appears to obey a simple rule: {∞} → {M/2} when K is even, and {∞} → {M } when K is odd. When M is prime, the decomposition of the BAE is very simple, and there is no {∞}. It remains to discuss the cases of odd nonprime M ≤ 18, namely M = 9 and M = 15. For M = 9 we find that {∞} → {3, 6} when K = 3 (see Table 6), and for M = 15 we have {∞} → {5, 10} when K = 3 (see Table 12). But we do not presently know how to establish such assignments for general M, K, without going through the actual computations of regularization via the T Q-relations.
To summarize, the computations described in this subsection appear to be an efficient short-cut for obtaining the decomposition dimension counting of Tables 3-15, without ever actually using the integrability of the XXX chain, analysing the BAE, or doing any algebraic geometry. This suggests that the solutions of the BAE simply decompose in a way that respects the conservation of spin and momentum, and the su(2) symmetry of the XXX chain, with no extra hidden structure. But obviously the decomposition of the BAE goes much further than the mere counting of dimensions; in particular the explicit results for the Gröbner bases make possible the efficient computations of the partition functions, as we have seen.

Conclusions and discussions
In this paper we developed a method to compute the torus partition function of the sixvertex model exactly and analytically. The method is based on an algebro-geometrical approach to the BAE, together with new ingredients that include the rational Q-system, primary decomposition, algbraic extension and Galois theory. The exact partition functions are polynomials in the spectral parameter z of order M N with rational coefficients. When M and N become large, we obtain polynomials with high orders and large coefficients. Since polynomials are essentially specified by their zeros, we solved for the zeros of the partition functions numerically to high accuracy and studied their behavior in the partial thermodynamic limit where N M and M is fixed. We observed that the zeros accumulate on some curves in this limit and gave a numerical method to generate the limiting curves of accumulation points. These curves exhibit some universal behaviors for different values of even and odd M which led us to formulate several observations and conjectures.
There are many open questions and new directions that one can pursue in the near future. We discuss some of them in what follows.
An immediate interesting direction is to generalize the current work to the the quantum deformed case. In this paper, we focussed on the six-vertex model at the isotropic point where the model is equivalent to the Heisenberg XXX spin chain. Away from the isotropic point, the six-vertex model is still integrable and is equivalent to the XXZ spin chain. In the XXZ spin chain, we have a new parameter q which is related to the anisotropy. The isotropic point corresponds to q = 1. Usually the BAE of the XXZ spin chain are written in terms of hyperbolic or trigonometric functions, and one might wonder how our approach, which seems to be restricted to rational functions, can be applied to this case. It is actually quite simple to perform a change of variables to recast the BAE in terms of rational functions.
To study the solution space of the BAE and the torus partition function as in this paper, we will however have to deal with several very interesting new features.
• First of all, one needs to distinguish between the cases where q takes a generic complex value and the cases where q is a root of unity. It is well known that the latter case is much more subtle than the former, in terms of solutions of BAE. The completeness problem for the generic q case is a straightforward generalization of the XXX case, namely the physical solutions consist of regular and physical singular solutions [25]. On the other hand, when q is a root of unity, due to the presence of the so-called exact K-strings, there are seemingly infinitely many solutions and the situation for the completeness problem is less clear. It is therefore not clear what are the physical solutions. Before we can compute the torus partition function, it seems that we need to sort out clearly the completeness problem first, which is an interesting question in its own right. Some preliminary calculations show that algebro-geometric methods in these cases are again very useful. For example, we observe that the Gröbner bases exhibit singularities when q is a root of unity and the quotient rings become affine varieties with positive dimensions (instead of a collection of points).
• The six-vertex model is closely related to another famous model, namely the Potts model. The latter can be represented in terms of the affine Temperley-Lieb (TL) algebra (see, e.g., [24,26] for a recent overview). The dimensions N M,K given in (3.6) appear naturally as the dimension of standard modules of the affine TL algebra, which we denote by W j,z (for the isotropic case, we take z = 1). The representations of the affine TL algebra take the graphical form of "link patterns" with pairwise connections (arcs) and defect lines (through-lines) where the total number of lines and the number of defect lines j play the role of length M and magnon number K, respectively, in our context. To compute the torus partition function of such models, 4 one can perform a decomposition within the standard module with respect to the lattice momentum [24], which is essentially the same as what we did in section 6.4. Let us recall that in our case the decomposition comes completely from studying the solution space of the BAE (or the rational Q-system) using algebraic geometry and the dimensions of the subspaces come from counting the number of solutions; while in the Potts model case, these come from studying the representation theory of the affine TL algebra. These similarities are quite remarkable and imply that the physical solutions of the BAE, studied here using the algebro-geometric approach, actually know a lot about the representation theory of affine TL algebra. It will be interesting to see to which extent these connections carry over to the q-deformed case.
• For generic q, we need to consider the standard module W j,z with non-trivial z, namely z = 1. In this case, there is another quantum number that appears which is associated with z. This is related to the momentum with which the defects spiral around the periodic direction. It will be intriguing to see how such a new quantum number can appear in our context by studying the solutions space of the BAE.
• The case when q is a root of unity is even more interesting. In that case, representations of affine TL should be reducible, but indecomposable. In practice, this will mean that the W j,z will have to be "glued" in various ways [27,28]. The complexity of these gluings and the expected appearance of Jordan cells will challenge the algebrogeometric approach. It will be very exciting to see how these structures carry over to the solution space of the BAE.
In the current work, we considered periodic boundary condition in both directions for the lattice. This corresponds to the torus partition function. It is also interesting to consider the partition function on other topologies, such as an annulus [6]. For this topology, we need to compute the transfer matrix of the spin chain with open boundary conditions. One nice starting point for this case is the quantum group invariant XXZ spin chain [29]. This spin chain is invariant under the quantum group U q (sl (2)). It has several nice properties. In particular, the completeness problem has been studied systematically in [30] both for generic q, and q at roots of unity. The relation with TL algebra has also been established. Working out this simpler example should also shed light on the more challenging periodic boundary conditions mentioned above.
For the zeros of partition function, it will be desirable to find an analytic approach to understand or even predict the condensation curves in the partial thermodynamic limit. It is also interesting to see how the q-deformation affects the distribution of the partition function zeros.
3. Some of the Q-functions at the boundary are completely fixed and do not need to be solved. The Q-functions at the upper right boundary are completely fixed to be 1 since there are no Bethe roots. In addition, the Q-polynomial at node (0, 0) is given by u M .
4. The rest of the Q-functions are determined by the QQ-relations where Q a,s denotes the Q-function associated with the node at position (a, s). The QQ-relation (A.1) is a relation between the four nodes around a box.
Here and in what follows, we introduce the shorthand notations To determine these polynomials, one makes an ansatz for the unknown Q-polynomials along certain path (see the example in the next subsection) and require that all the Q-functions on the Young tableaux are polynomials. This leads to a set of algebraic equations for the unknown coefficients. Solving this set of algebraic equations gives the Q-polynomials.

The
Bethe roots of length M and magnon number K are given by the zeros of Q 0,1 (u). This fact will be shown below.
In order to explain the above procedure, we give an explicit example in the next subsection.

A.1 A simple example
In order to explain the method, we consider a simple example with M = 6 and K = 2, where M is the length of the spin chain and K is the number of magnons. The corresponding Young diagram is given by (M − K, K) = (4, 2), and is shown in figure A.8. The number on each node denotes the degree of the corresponding Q-polynomial. There are thirteen Q-functions on the Young tableaux, each Q-polynomial being associated with a node of the Young tableaux. The boundary nodes (the ones labeled by 0) are simply taken to be 1. In addition, the Q-polynomial at the origin is taken to be Q 0,0 (u) = u 6 . (A.3) One needs to choose a path from the origin to the upper right boundary. We have chosen the one with yellow color in figure A.8. Apart from the two boundary Q-polynomials, there are two other unknown Q-polynomials which are parameterized by The remaining task is to use the QQ-relations (A.1) to determine the remaining unknown Q-polynomials as well as the coefficients c Let us first determine Q 1,0 (u). Taking a = s = 0, the QQ-relation leads to where c is some normalization constant to make the polynomial monic.
1,1 . Since by definition Q 1,0 (u) is a polynomial, the remainder should be zero. This leads to a set of equations for the unknown coefficients. These relations are called zero remainder conditions (ZRC). Repeating this analysis for all the other non-trivial Q-functions, we obtain the full set of ZRC, which are the systems of equations that we need to solve.
The paths can be chosen in different ways, which result in different forms of ZRC, but finally they lead to the same solution of rational Q-systems. The simplest choice of the path is the one that goes from (0, 0) to (0, 2) and then from (0, 2) to the rightmost node (K, 2). In this way, we only have one unknown Q-function, namely Q 0,1 , to determine.
After solving the ZRC, we find all the Q-polynomials. The Bethe roots are simply given by the zeros of Q 0,1 (u). This is proved in the next subsection.
A.2 From Q-system to BAE: su(2) spin chain In the previous section, we discussed how to find the solutions of all the Q-functions once the boundary conditions are fixed. In this section, we derive the BAE from the Q-system, which will demonstrate why the zeros of Q 0,1 are identified with the Bethe roots.
To this end, let us consider a generic Young tableaux with two rows of (M − K, K) boxes, as shown in figure A.9. We have K ≤ M − K. Let us consider the Q-function Q 0,1 (u) (the one at the node with the orange circle) and the related QQ-relations. From the power counting we know that Q 0,1 is a monic polynomial of order K. We assume that its roots are u 1 , · · · , u K and write We consider the QQ-relations in the two shaded boxes, which read At u = u k , Q 0,1 (u k ) = 0, and by shifting (A.7) by ±i/2 we obtain the following relations: Using the boundary condition Q 0,2 = Q 1,2 = 1, we obtain Meanwhile, evaluating (A.8) at u = u k produces and plugging (A.10) into (A.11) we obtain Use the boundary condition and the expression for Q 0,1 (u), we finally arrive at which is nothing but the BAE of SU (2) spin chain with length M and magnon number K.
It follows in particular that the roots u j of Q 0,1 (u) are precisely the Bethe roots, so we can identify the latter with the Baxter polynomial (2.14), viz. Q 0,1 (u) = Q(u), as previously claimed.

B More on algebraic geometry
In this appendix, we briefly review the important algebraic geometry technique used in our paper, primary decomposition. The mathematics reference is [31]. For the concept of Gröbner basis and companion matrix, we refer to the introduction in [2].
Let F be a field and A = F[z 1 , . . . z n ] be a polynomial. Any ideal I in A decomposes into the intersection of several primary ideals.
where each I j is primary. A primary ideal J is an ideal such that if a, b ∈ A, ab ∈ J, then either a ∈ J or b n ∈ J. The decomposition (B.1) is an analogy for the factorization of integers.
Note that with the decomposition (B.1), geometrically, the zero sets of I decompose into the union of several algebraic sets Z(I) = Z(I 1 ) ∪ Z(I 2 ) ∪ . . . ∪ Z(I k ) . (

B.2)
This property is useful for the study of the complicated zero set Z(I). Especially, when Z(I) is zero-dimensional, like the set of Bethe roots, this decomposition classifies the Bethe roots and we shall call each I j decomposed BAE.
Note that the definition of primary decomposition depends on the coefficient field F. For example, within Q[x, y], is already primary and hence cannot be decomposed further. However, within the algebraic closureQ[x, y], we have the primary decomposition: We refer to [31] for the introduction of algebraic extension and Galois theory. The computation of primary decomposition can be carried out by the computer algebra system Singular [16].

C Power of companion matrices
The main method we present in this paper to compute partition function is to calculate the companion matrices for BAE. Recall that, according to (3.14), the six-vertex model partition function is given by where T M,K (z) is the companion matrix of the polynomial t M,K (s, z). Recall also that T M,K (z), calculated from the Gröbner basis, contains only rational numbers. Therefore the whole computation is manifestly analytic.
In practice, although T M,K (z) can be calculated from the straightforward Gröbner basis and polynomial division procedure, the matrix power T M,K (z) N computation can be difficult.
We here present an alternative algorithm which speeds up the computation and saves RAM usage. The algorithm can be sketched as follows: 1. For given M and K, calculate the Gröbner basis G(I M,K ). capable of sorting the eigenvalues in order of decreasing norm.
A huge advantage of such iterative methods is that T M (z) is given as a product of sparse matrices via (2.6), where each R-matrix contains at most two non-zero entries per column. Therefore the computation of v in each iteration requires at most 2M d operations, where d is the dimension of the matrix. The trace Tr a is performed by going from M to M + 2 sites when the auxiliary space is inserted, and back to M sites after a row of the lattice has been completed and the trace operation performed.
If the magnon and momentum labels (K, ) for the equimodular eigenvalues Λ 1 (z) and Λ 2 (z) were known beforehand, it would obviously be most efficient to diagonalize the smaller matrices T M,K, (z) constructed in (6.37). But since the accumulation curves in practice turn out to have multiple branches and T-points where the sector labels may change, this approach would necessitate a considerable amount of manual intervention. We have thus chosen a more brute-force approach in which the entire matrix T M (z) is diagonalized, for all K = 0, 1, . . . , M/2 simultaneously and with decomposing the momentum with respect to the label. Note that the su(2) highest-weight constraint is also not enforced in this approach, so the eigenvalues with K < M/2 present degeneracies. We deal with this in practice by imposing the equimodularity criterion |Λ 1 (z)| = |Λ r (z)|, where r ≥ 2 is a suitably chosen (small) integer, which may need some adjustment as we run through the various branches of the equimodular curves.
The first step in our procedure is to acquire some approximate knowledge about where to start the search for the equimodular curves. In the case of Figure 5.5 we first made a rough plot of the norms of the first few eigenvalues along a few straight lines with constant Re z, or along the real axis, to get a finite list of points close to the equimodular curves. In a second step, we then launched a direct-search algorithm, taking each of these points as the initial point. The direct-search method is carefully described in [3]; it has the property of first locking onto any close-by equimodular curve and then following it in small steps. 7 The search is stopped whenever a branch of the equimodular wanders off to infinity, or if it starts overlapping with a part of the curve which is already known. One all starting points have been exploited, we have completed the second step.
The third and final step consists in making sure that the set of equimodular curves is complete. This requires in particular examining carefully the surroundings of any point where the curves present a discontinuous tangent vector, since this is the sign of a T-point or a higher-order bifurcation point. As in the first step, we make a rough plot of the norms along a small circle surrounding any potential bifurcation point. If a branch is identified that has not yet been traced out, we go back to the second step as many times as necessary.
Obviously, if the whole set of equimodular curves has some very small disconnected pieces, they may be missed by this procedure. However, the fact that the curves in Figure 5.5 present only a single connected component gives appealing evidence that they are actually complete.

E More results of primary decomposition on F M
In this appendix, we list the results of primary decomposition on F M for 9 ≤ M ≤ 18.