<figure><img src="image_1.png" /></figure>

HAL
open science

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

Jesper Lykke Jacobsen, Yunfeng Jiang, Yang Zhang

► To cite this version:

Jesper Lykke Jacobsen, Yunfeng Jiang, Yang Zhang. Torus partition function of the six-vertex model from algebraic geometry. Journal of High Energy Physics, 2019, 2019 (03), pp.152. 10.1007/JHEP03(2019)152. hal-01960385

HAL Id: hal-01960385

https://hal.science/hal-01960385v1

Submitted on 21 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Torus partition function of the six-vertex model from algebraic geometry

Jesper Lykke Jacobsen,ª,ᵇ,ᶜ Yunfeng Jiangª,ᵈ,ᵉ and Yang Zhangª,*f

E-mail: jesper.jacobsen@ens.fr, yunfeng.jiang@cern.ch, yzhang@mpp.mpg.de

**ABSTRACT:** 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 algebraic 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.

**KEYWORDS:** Bethe Ansatz, Differential and Algebraic Geometry, Lattice Integrable Models

ARXIV ePRINT: 1812.00447

<footnote>ªLaboratoire de Physique Théorique, Département de Physique de l'ENS,
École Normale Supérieure, Sorbonne Université, CNRS,
PSL Research University, 75005 Paris, France</footnote>

<footnote>ᵇSorbonne Université, École Normale Supérieure, CNRS,
Laboratoire de Physique Théorique (LPT ENS), 75005 Paris, France</footnote>

<footnote>ᶜInstitut de Physique Théorique, Paris Saclay, CEA,
CNRS, 91191 Gif-sur-Yvette, France</footnote>

<footnote>ᵈInstitut für Theoretische Physik, ETH Zürich,
Wolfgang Pauli Strasse 27, CH-8093 Zürich, Switzerland</footnote>

<footnote>ᵉCERN Theory Department, Geneva, Switzerland</footnote>

<footnote>fMax Planck Institut für Physik Werner Heisenberg Institut,
80805 München, Germany</footnote>

<footnote n="9">⁹PRISMA Cluster of Excellence, Johannes Gutenberg University,
55128 Mainz, Germany</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Torus partition function of the six-vertex model</td><td>4</td></tr><tr><td>3</td><td>Algebro-geometric approach to the partition function</td><td>8</td></tr><tr><td>4</td><td>Explicit results</td><td>11</td></tr><tr><td>4.1</td><td>Closed-form expressions for M ≤ 6</td><td>11</td></tr><tr><td>4.2</td><td>Partition functions for higher M and N</td><td>14</td></tr><tr><td>4.3</td><td>Consistency check</td><td>15</td></tr><tr><td>5</td><td>Zeros of partition functions</td><td>16</td></tr><tr><td>5.1</td><td>Partition functions for different M and N</td><td>17</td></tr><tr><td>5.2</td><td>Limiting curves</td><td>19</td></tr><tr><td>6</td><td>Primary decomposition</td><td>21</td></tr><tr><td>6.1</td><td>Primary decomposition over Q</td><td>23</td></tr><tr><td>6.2</td><td>Primary decomposition over F<sub>M</sub></td><td>25</td></tr><tr><td>6.3</td><td>Results on decomposition over F<sub>M</sub></td><td>29</td></tr><tr><td>6.4</td><td>Relation to naive momentum-sector diagonalization</td><td>30</td></tr><tr><td>7</td><td>Conclusions and discussions</td><td>33</td></tr><tr><td>A</td><td>Rational Q-system</td><td>35</td></tr><tr><td>A.1</td><td>A simple example</td><td>36</td></tr><tr><td>A.2</td><td>From Q-system to BAE: su(2) spin chain</td><td>37</td></tr><tr><td>B</td><td>More on algebraic geometry</td><td>39</td></tr><tr><td>C</td><td>Power of companion matrices</td><td>39</td></tr><tr><td>D</td><td>Details on finding limiting curves</td><td>40</td></tr><tr><td>E</td><td>More results of primary decomposition on F<sub>M</sub></td><td>42</td></tr></table>
# 1 Introduction

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 \times 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 $\mathbb{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 on an algebraically extended field. The structure we find is interesting on its own right. Furthermore, as an application, by using this decomposition and Galois theory, we gain a huge boost in the efficiency of the partition function computation.

Let us summarize the main results of this paper.

1. We perform the decomposition with respect to total momentum on an extended field $\mathbb{F}_M$ (defined in section 6). We denote each subspace by $I_{M,K,l}$ where $M$ is the length of the spin chain, or equivalently the size of one direction of the lattice; $K$ is
the number of magnons and $ℓ$ denotes the momentum sector with total momentum $2\piℓ/M$. We compute the Gröbner basis and construct the quotient ring for each $I_{K,M,ℓ}$. The dimensions of the quotient rings, or equivalently the number of physical solutions of the BAE are given in tables 3–15 for length $6 \le M \le 18$.

2. We compute the companion matrices of the transfer matrix $T_{M,K,ℓ}(z)$ and Baxter's $Q$-polynomials $Q_{M,K,ℓ}(z)$ within each subspace $I_{M,K,ℓ}$ up to $M=14$. Using this data, one can compute the exact torus partition function for any number of $N$. Of course, as $N$ grows, the results become quickly large (meaning here high-degree polynomials with huge rational coefficients). In practice, we can compute the partition function up to rather large $N \sim 100$ for all $M$ up to $M=14$. In addition, we can also diagonalize these matrices numerically and find Bethe roots and eigenvalues of the transfer matrix up to $M=14$ efficiently.

3. Our results for the partition function with $N \gg M$ are close to the partial thermodynamic limit of fixed $M$ and $N \to \infty$. Powerful techniques for studying this limit have been developed in the framework of the (antiferromagnetic) Potts model [3–9]. Building on elements of these works, we here study the partition function zeros for the isotropic six-vertex model in this limit. As $N$ grows, the zeros accumulate on curves in accordance with the Beraha-Kahane-Weiss theorem. We give a numerical method to determine these limiting curves of accumulation points and discuss the universal behaviors of the curves for different values of $M$.

All the results mentioned above can be downloaded from the github repository,

https://github.com/yzhphy/BAE_AG/tree/master/Results/Exact_Partition_Function/partition_function_M14

For the future reference, we also present the Gröbner basis results for the six-vertex model,

https://github.com/yzhphy/BAE_AG/tree/master/Results/XXX_Groebner_basis

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 \le 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 $\mathbb{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
<figure><img src="image_3.png" /><figcaption>Figure 1. Configurations of 6-vertex model.</figcaption></figure>

technical aspects, and the tables giving the number of physical solutions of the BAE for $6 \le M \le 18$ are relegated to appendix E.

## 2 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 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 $\epsilon_i(z, \theta)$, ($i = 1, \dots, 6$) subject to the following constraints

$$ \epsilon_1 = \epsilon_2, \quad \epsilon_3 = \epsilon_4, \quad \epsilon_5 = \epsilon_6. \qquad (2.1) $$

Following Baxter [1], we denote the corresponding Boltzmann weights by $\omega_j = \exp(-\beta\epsilon_j)$ and define

$$ a = \omega_1 = \omega_2, \quad b = \omega_3 = \omega_4, \quad c = \omega_5 = \omega_6, \qquad (2.2) $$

where in the isotropic case

$$ a(z, \theta) = z - \theta + i/2, \quad b(z, \theta) = z - \theta - i/2, \quad c(z, \theta) = i. \qquad (2.3) $$

Indeed, the anisotropy parameter is then $\Delta = \frac{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. The partition function is given by summing the product of local Boltzmann weights over all possible configurations.
<figure><img src="image_4.png" /><figcaption>Figure 2. One possible configuration for the six-vertex model of size $M \times N$.</figcaption></figure>

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, \theta) = \begin{pmatrix} a(z, \theta) & 0 & 0 & 0 \\ 0 & b(z, \theta) & c(z, \theta) & 0 \\ 0 & c(z, \theta) & b(z, \theta) & 0 \\ 0 & 0 & 0 & a(z, \theta) \end{pmatrix}. \qquad (2.4)$$

Integrability of the six-vertex model is guaranteed by the fact that the R-matrix satisfies the Yang-Baxter equation

$$R_{an}(u)R_{bn}(v)R_{ab}(u-v) = R_{ab}(u-v)R_{bn}(v)R_{an}(u). \qquad (2.5)$$

The transfer matrix is defined as

$$T_M(z) = \mathrm{Tr}_a \left( \prod_{n=1}^{M} R_{an}(z, \theta_n) \right). \qquad (2.6)$$

The partition function can be written in terms of the transfer matrix

$$Z_{M,N} = \mathrm{Tr} [T_M(z_1)T_M(z_2) \cdots T_M(z_N)], \qquad (2.7)$$

where the parameters $z_j$ and $\theta_k$ characterize the Boltzman weights at each site of the lattice. We consider the homogeneous case $\theta_k = 0$ and $z_j = z$. Then the partition function is simply given by

$$Z_{M,N} = \mathrm{Tr} [T_M(z)^N]. \qquad (2.8)$$

It is clear from the above definitions that it is a polynomial of degree $MN$ 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. 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

$$d_{M,K} = \binom{M}{K}. \qquad (2.9)$$

We denote the transfer matrix of the spin sector labelled by K as $T_{M,K}(z)$. The partition function is then given by

$$Z_{M,N} = \sum_{K=0}^{M} \mathrm{Tr} [T_{M,K}(z)^N]. \qquad (2.10)$$

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\rangle$, parameterized by a set of parameters $u = \{u_1, \dots, u_K\}$ called *rapidities*. We then have

$$T_{M,K}(z)|u\rangle = t_{u}(z)|u\rangle. \qquad (2.11)$$

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

$$t_{u}(z) = a(z) \frac{Q_{u}(z-i)}{Q_{u}(z)} + d(z) \frac{Q_{u}(z+i)}{Q_{u}(z)}, \qquad (2.12)$$

where

$$a(z) = \left(z + \frac{i}{2}\right)^M, \qquad d(z) = \left(z - \frac{i}{2}\right)^M, \qquad (2.13)$$

and

$$Q_{u}(z) = \prod_{j=1}^{K} (z - u_{j}) \qquad (2.14)$$

is called the *Baxter polynomial*. The rapidities are constrained by the Bethe ansatz equations (BAE)

$$\left(\frac{u_j + i/2}{u_j - i/2}\right)^M = \prod_{k \neq j}^{K} \frac{u_j - u_k + i}{u_j - u_k - i}, \qquad j = 1, \dots, K. \qquad (2.15)$$
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 \to \infty$. However, for fixed and finite $M$ and $N$, the expression for $t_u(z)$ is only formal since the parameters $\{u_1, \dots, 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 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 $TQ$-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 $TQ$-relation. In terms of the eigenvalues of the $T$ and $Q$ operators, the $TQ$ relation becomes the following functional equation

$$t(z)Q(z) = a(z)Q(z-i) + d(z)Q(z+i). \qquad (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 $TQ$-relation (2.16), one first writes $t(z)$ and $Q(z)$ in the explicit polynomial form

$$t(z) = \sum_{j=0}^{M} t_j z^j, \qquad Q(z) = z^K + \sum_{k=0}^{K-1} s_k z^k. \qquad (2.17)$$

The unknown variables that we are solving for are $\{t_0, \dots, t_M, s_0, \dots, 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, \dots, t_M\}$ and $\{s_0, \dots, s_{K-1}\}$ and are easier to solve than the original set of BAE. Solving the $TQ$-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 *TQ*-relation is that it automatically eliminates the solutions with coinciding rapidities, i.e., the solutions of the *TQ*-relation lead to pairwise distinct solutions of the BAE.

Working with the *TQ*-relation instead of the original BAE is already more efficient for our purpose. However, the use of the *TQ*-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, \dots, 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 *TQ*-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 ($\Delta = 1$) corresponding to the XXX spin chain. To treat the XXZ spin chain at generic $\Delta$, we still need to rely on the *TQ*-relation, together with additional constraints to select physical solutions. The strategy is to first find the *Q*-polynomial, and next use the *TQ*-relation 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, *TQ*-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

$$Q(z) = Q_{0,1}(z) = z^K + \sum_{k=0}^{K-1} s_k z^k. \qquad (3.1)$$

**Ideal.** The zero-remainder conditions (ZRC) then give a set of algebraic equations for the $K$ variables $\{s_0, \dots, s_{K-1}\}$,

$$f_1(s_0, s_1, \dots, s_{K-1}) = f_2(s_0, s_1, \dots, s_{K-1}) = f_S(s_0, s_1, \dots, s_{K-1}) = 0, \qquad (3.2)$$

where $f_k(s_0, s_1, \dots, s_{K-1})$ are polynomials in the variables $\{s_0, \dots, s_{K-1}\}$. Here $S$ is the number of equations and it depends on the path we choose. The polynomials $f_1, \dots, f_S$ define an ideal in the polynomial ring $\mathbb{C}[s_0, s_1, \dots, s_{K-1}]$, denoted

$$I_{M,K} = \langle f_1, \dots, f_S \rangle. \qquad (3.3)$$

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$

$$I_{M,K} = \langle f_1, \dots, f_S \rangle = \langle G_1, \dots, G_{S'} \rangle \qquad (3.4)$$

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, \dots, s_{K-1}\}$. For different orderings, the corresponding Gröbner basis can look quite different.

**Quotient ring.** The quotient ring is defined as

$$Q_{M,K} = \mathbb{C}[s_0, s_1, \dots, s_{K-1}]/I_{M,K}, \qquad (3.5)$$

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]

$$\mathcal{N}_{M,K} = \binom{M}{K} - \binom{M}{K-1}. \qquad (3.6)$$

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, \dots, s_{K-1}\}$ that cannot be divided by $LT[G_k]$ ($k = 1, \dots, 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 $\langle f_1, \dots, f_S \rangle$. 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, \dots, e_{N_{M,K}}$}. Any polynomial $P(s_0, s_1, \dots, 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,

$$P e_j = \sum_{k=1}^{S'} a_k G_k + r_j. \qquad (3.7)$$

The remainder $r_j(s_0, \dots, 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.,

$$r_j = \sum_{k=1}^{N_{M,K}} M_{jk} e_k, \quad M_P = (M_{ij}). \qquad (3.8)$$

The companion matrix satisfies the following homomorphism properties

$$\begin{aligned}
M_{P_1 \pm P_2} &= M_{P_1} \pm M_{P_2}, \\
M_{P_1 P_2} &= M_{P_1} \cdot M_{P_2}, \\
M_{P_1 / P_2} &= M_{P_1} \cdot M_{P_2}^{-1}.
\end{aligned} \qquad (3.9)$$

The main result from algebraic geometry which we are going to use is

$$\sum_{\text{sol}} P(s_0, \dots, s_{K-1}) = \text{Tr} M_P, \qquad (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 $\mathbb{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 \mapsto S_k$, so that

$$Q_{0,1}(z) \mapsto Q_{M,K}(z) = z^K \mathbf{I} + \sum_{k=0}^{K-1} S_k z^k. \qquad (3.11)$$

where $\mathbf{I}$ is the identity matrix of dimension $N_{M,K}$. Our goal is to find the companion matrix of the transfer matrix $t_{\mathbf{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

$$t_{\mathbf{u}}(z) \mapsto T_{M,K}(z), \qquad (3.12)$$

where $T_{M,K}(z)$ can be computed using the ideal generated by the TQ-relations<fnref n="1" /> or from the companion matrix $Q_{M,K}(z)$ by the following relation

$$T_{M,K}(z) = [a(z)Q_{M,K}(z-i) + d(z)Q_{M,K}(z+i)] \cdot Q_{M,K}(z)^{-1}, \qquad (3.13)$$

<footnote n="1">¹More precisely, to select physical solutions, we need to combine the equations coming from the TQ-relations and the rational Q-system and then eliminate the variables s.</footnote>
valid whenever $\mathbf{Q}_{M,K}(z)$ is non-singular. The partition function is given by

$$Z_{M,N} = \sum_{K=0}^{[M/2]} (M - 2K + 1) \mathrm{Tr} (\mathbf{T}_{M,K}(z)^N). \qquad (3.14)$$

Several comments are in order. The multiplicities $M-2K+1$, for $K \le 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 $\mathbf{T}_{M,K}(z)$ are smaller than those of $T_{M,K}(z)$. The dimensions of $\mathbf{T}_{M,K}(z)$ and $T_{M,K}(z)$ are $\mathcal{N}_{K,M}$ and $d_{K,M}$ respectively, see (3.6) and (2.9). For example, for $M=14$, $K=7$, we have $\mathcal{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 $\mathbf{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.

# 4 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 $\mathbb{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 \ge 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.

## 4.1 Closed-form expressions for $M \le 6$

In this section, we list the results for $M \le 6$. We denote the partition function of an $M \times 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

$$Z_{1,N} = 2 \times (2z)^N. \qquad (4.1)$$
M = 2. We have to sum over $K = 0$ and $K = 1$. The partition function is given by

$$Z_{2,N} = 3 \left( 2z^2 - \frac{1}{2} \right)^N + \left( 2z^2 + \frac{3}{2} \right)^N. \qquad (4.2)$$

M = 3. We again have to sum over $K = 0$ and $K = 1$. The result is

$$Z_{3,N} = 4 \left( 2z^3 - \frac{3}{2}z \right)^N + 2 \left[ \left( 2z^3 + \frac{3}{2}z + \frac{\sqrt{3}}{2} \right)^N + \left( 2z^3 + \frac{3}{2}z - \frac{\sqrt{3}}{2} \right)^N \right]. \quad (4.3)$$

We see that irrational numbers start to show up within some eigenvalues of the transfer matrix. However, for any $N \in \mathbb{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

$$Z_{4,N} = 5 \left( 2z^4 - 3z^2 + \frac{1}{8} \right)^N + 3 \left( 2z^4 + z^2 + 2z + \frac{1}{8} \right)^N + 3 \left( 2z^4 + z^2 - 2z + \frac{1}{8} \right)^N + 3 \left( 2z^4 + z^2 - \frac{7}{8} \right)^N + \left( 2z^4 + 3z^2 + \frac{13}{8} \right)^N + \left( 2z^4 + 3z^2 - \frac{3}{8} \right)^N. \qquad (4.4)$$

M = 5. We have to sum over $K = 0, 1, 2$. The final result reads

$$Z_{5,N} = 6 \left( 2z^5 - 5z^3 + \frac{5}{8}z \right)^N + 4 \left( 2z^5 - \frac{1}{2}\sqrt{25+10\sqrt{5}}z^2 - \frac{1}{8}(5-4\sqrt{5})z + \frac{1}{8}\sqrt{5-2\sqrt{5}} \right)^N + 4 \left( 2z^5 + \frac{1}{2}\sqrt{25+10\sqrt{5}}z^2 - \frac{1}{8}(5-4\sqrt{5})z - \frac{1}{8}\sqrt{5-2\sqrt{5}} \right)^N + 4 \left( 2z^5 - \frac{1}{2}\sqrt{25-10\sqrt{5}}z^2 - \frac{1}{8}(5+4\sqrt{5})z + \frac{1}{8}\sqrt{5+2\sqrt{5}} \right)^N + 4 \left( 2z^5 + \frac{1}{2}\sqrt{25-10\sqrt{5}}z^2 - \frac{1}{8}(5+4\sqrt{5})z - \frac{1}{8}\sqrt{5+2\sqrt{5}} \right)^N + 2 \left( 2z^5 + 3z^3 + \frac{21}{8}z \right)^N + 2 \left( 2z^5 + 3z^3 - \frac{1}{2}\sqrt{10-2\sqrt{5}}z^2 + \frac{1}{8}(1+4\sqrt{5})z - \frac{1}{8}\sqrt{50+22\sqrt{5}} \right)^N + 2 \left( 2z^5 + 3z^3 + \frac{1}{2}\sqrt{10-2\sqrt{5}}z^2 + \frac{1}{8}(1+4\sqrt{5})z + \frac{1}{8}\sqrt{50+22\sqrt{5}} \right)^N + 2 \left( 2z^5 + 3z^3 - \frac{1}{2}\sqrt{10+2\sqrt{5}}z^2 + \frac{1}{8}(1-4\sqrt{5})z + \frac{1}{8}\sqrt{50-22\sqrt{5}} \right)^N + 2 \left( 2z^5 + 3z^3 + \frac{1}{2}\sqrt{10+2\sqrt{5}}z^2 + \frac{1}{8}(1-4\sqrt{5})z - \frac{1}{8}\sqrt{50-22\sqrt{5}} \right)^N.$$

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

$$
\begin{align}
Z_{6,N} ={}& 7 \left( 2z^6 - \frac{15}{2}z^4 + \frac{15}{8}z^2 - \frac{1}{32} \right)^N \tag{4.6} \\
& + 5 \left( 2z^6 - \frac{3}{2}z^4 - \frac{25}{8}z^2 + \frac{11}{32} \right)^N \nonumber \\
& + 5 \left( 2z^6 - \frac{3}{2}z^4 + \sqrt{3}z^3 - \frac{21}{8}z^2 - \frac{3\sqrt{3}}{4}z - \frac{1}{32} \right)^N \nonumber \\
& + 5 \left( 2z^6 - \frac{3}{2}z^4 - \sqrt{3}z^3 - \frac{21}{8}z^2 + \frac{3\sqrt{3}}{4}z - \frac{1}{32} \right)^N \nonumber \\
& + 5 \left( 2z^6 - \frac{3}{2}z^4 + 3\sqrt{3}z^3 + \frac{11}{8}z^2 - \frac{\sqrt{3}}{4}z - \frac{1}{32} \right)^N \nonumber \\
& + 5 \left( 2z^6 - \frac{3}{2}z^4 - 3\sqrt{3}z^3 + \frac{11}{8}z^2 + \frac{\sqrt{3}}{4}z - \frac{1}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 - \frac{25}{8}z^2 + \frac{3}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 + \frac{1}{8}(15 - 8\sqrt{5})z^2 - \frac{1}{32}(21 + 8\sqrt{5}) \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 + \frac{1}{8}(15 + 8\sqrt{5})z^2 - \frac{1}{32}(21 - 8\sqrt{5}) \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 - \sqrt{3}z^3 + \frac{11}{8}z^2 - \frac{5\sqrt{3}}{4}z - \frac{9}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 + \sqrt{3}z^3 + \frac{11}{8}z^2 + \frac{5\sqrt{3}}{4}z - \frac{9}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 - \frac{\sqrt{54 - 6\sqrt{17}}}{2}z^3 + \frac{(2\sqrt{17} - 3)}{8}z^2 \right. \nonumber \\
& \qquad \left. - \frac{\sqrt{54 - 6\sqrt{17}}(3 + \sqrt{17})}{16}z + \frac{(9 + 2\sqrt{17})}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 + \frac{\sqrt{54 - 6\sqrt{17}}}{2}z^3 + \frac{(2\sqrt{17} - 3)}{8}z^2 \right. \nonumber \\
& \qquad \left. + \frac{\sqrt{54 - 6\sqrt{17}}(3 + \sqrt{17})}{16}z + \frac{(9 + 2\sqrt{17})}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 - \frac{\sqrt{54 + 6\sqrt{17}}}{2}z^3 - \frac{(2\sqrt{17} + 3)}{8}z^2 \right. \nonumber \\
& \qquad \left. - \frac{\sqrt{54 + 6\sqrt{17}}(\sqrt{17} - 3)}{16}z + \frac{(9 - 2\sqrt{17})}{32} \right)^N \nonumber \\
& + 3 \left( 2z^6 + \frac{5}{2}z^4 + \frac{\sqrt{54 + 6\sqrt{17}}}{2}z^3 - \frac{(2\sqrt{17} + 3)}{8}z^2 \right. \nonumber \\
& \qquad \left. - \frac{\sqrt{54 + 6\sqrt{17}}(\sqrt{17} - 3)}{16}z + \frac{(9 - 2\sqrt{17})}{32} \right)^N \nonumber
\end{align}
$$
$$ 
\begin{aligned}
& + \left( 2z^6 + \frac{9}{2}z^4 + \frac{23}{8}z^2 - \sqrt{3}z - \frac{1}{32} \right)^N + \left( 2z^6 + \frac{9}{2}z^4 + \frac{23}{8}z^2 + \sqrt{3}z - \frac{1}{32} \right)^N \\
& + \left( 2z^6 + \frac{9}{2}z^4 + \frac{7-8\sqrt{13}}{8} + \frac{31-8\sqrt{13}}{32} \right)^N \\
& + \left( 2z^6 + \frac{9}{2}z^4 + \frac{7+8\sqrt{13}}{8} + \frac{31+8\sqrt{13}}{32} \right)^N \\
& + \left( 2z^6 + \frac{9}{2}z^4 + \frac{15}{8}z^2 - \frac{25}{32} \right)^N.
\end{aligned}
$$

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 \in \mathbb{N}$, we obtain polynomials whose coefficients are rational numbers, as expected.

## 4.2 Partition functions for higher $M$ and $N$

For $M = 7$ and $M = 8$, we can also work out the analytical results. However, the closed-form results involve complicated multi-square roots and are not very illuminating to write down explicitly here. For $M \ge 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 $T_{M,K}(z)$ up to $M = 14$ after algebraic extension (see more details in section 6). The dimensions of the companion matrices with fixed $M, K, \ell$ are given in tables 3–15. In general, the partition function $Z_{M,N}(z)$ is a polynomial of order $MN$. Let us consider one example and take $M = 14, N = 100$. We have

$$ Z_{14,100}(z) = \sum_{k=0}^{700} c_n z^{2n} \qquad (4.7) $$

where $c_n$ are rational numbers. The full result is too large to show, we present simply one coefficient, say $c_{50}$, here just to give an idea about the result. It can be written as $c_{50} = N_{50}/D_{50}$ where $N_{50}$ and $D_{50}$ are integers and are given by

$$ 
\begin{aligned}
N_{50} = & 3549714199509718765414261648405948375346908444631449814070999015959834 \\
& 5566548305771114077497148041577082024237243782436360433278999500001467 \\
& 7415020115130416461041186374238411444900889750187515469354178173296042 \\
& 6263542784005444527370571255655399436082869031810512249627664031748092 \\
& 4072971596196276059412147685707041129531252668420023067545577282800096 \\
& 9295580666703059631065928833005677223301695252356346794040665489002919 \\
& 252969420827675164522054378930548720509758103112303817657569605
\end{aligned}
$$
and

$$D_{50} = 5254662920039596746236382353531688074822975081008511328653941899536358 \\ 8535283647520133924979358051065937411251928593712733718650179494204788 \\ 1824032065695719746688948719970424831685787857667697873122163831147965 \\ 2780562037740551931521573661427623891169657954033123606560887017391291 \\ 1075991202517005852196882800252419269626961951377918107052351183818088 \\ 1989632.$$

Approximately this number is $c_{50} \approx 6.75536 \times 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.

## 4.3 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 3. Specifically, we consider the transformation of the partition function where we rotate the lattice by $\pi/2$. To analyze the effect of this transformation, we consider the $M \times N$ lattice with $M$ vertical lines and $N$ horizontal lines. We associate to each vertical line the same spectral parameter $\theta$, and to each horizontal line the same spectral parameter $z$. We denote this partition function by $Z_{M,N}(z - \theta)$. There are six possible configurations at each site with three different Boltzmann weights given by

$$a(z, \theta) = z - \theta + i/2, \quad b(z, \theta) = u - \theta - i/2, \quad c(z, \theta) = i. \qquad (4.8)$$
<figure><img src="image_5.png" /><figcaption>Figure 3. Lattice modular invariance of the lattice partition function.</figcaption></figure>

The partition function can be written as

$$Z_{M,N}(z-\theta) = \sum_{\substack{m,n,k \ge 0 \\ m+n+k=MN}} f_{m,n,k} a(z, \theta)^m b(z, \theta)^n c(z, \theta)^k, \quad (4.9)$$

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 $\pi/2$, it is clear, from the first line of figure 1, that we have the following transformation

$$a(u, \theta) \mapsto b(\theta, u) = -a(u, \theta), \quad (4.10)$$

$$b(u, \theta) \mapsto a(\theta, u) = -b(u, \theta),$$

$$c(u, \theta) \mapsto c(\theta, u) = c(u, \theta).$$

Therefore, under this transformation, which we denote by $\mathbb{R}$, the partition function transforms as

$$Z_{M,N}^{\mathbb{R}}(u-\theta) = Z_{N,M}(\theta-u) = \sum_{\substack{m,n,k \ge 0 \\ m+n+k=MN}} (-1)^{m+n} f_{m,n,k} a(u, \theta)^m b(u, \theta)^n c(u, \theta)^k. \quad (4.11)$$

The number of type-c vertices $k$ is even due to the ice rule. Therefore, if $MN$ 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-\theta)$ for all $N \le M$. These partition functions can also be computed as $Z_{N,M}(\theta-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$.

## 5 Zeros of partition functions

The torus partition function $Z_{M,N}(z)$ is a polynomial of order $MN$ 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–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 \to \infty$, 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 \times N$ system is obtained by letting $M, N \to \infty$ with a fixed and finite ratio, $0 < M/N < \infty$. But another means of obtaining relevant information is to fix a finite value of $M$, and study the accumulation points of zeros as $N \to \infty$. 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 \to \infty$). 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 \to \infty$.

## 5.1 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 \to \infty$ (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.

In figure 4 we show the partition function zeros for $M = 4, 5, \dots, 8$ with various aspect ratios, namely $N = \rho M$ for $\rho = 10, 20, 40, 80$ (and in one case $\rho = 160$). The lower right corner shows the result of our largest computation with $M = 14$ and $N = 100$. The partition functions for $M = 4, 5, 6$ have been computed by the direct application of the explicit formulae (4.4)–(4.6), and the remaining results by the algebra-geometric method. These are polynomials of degree $MN$, 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.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Partition function zeros in the complex z-plane, shown in reading direction for $M = 4, 5, \dots, 8$ and $M = 14$.</figcaption></figure>

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 4 reveal several interesting features. As the aspect ratio $\rho$ 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-$\rho$ effects are small, but further away their importance increases, and the fine structure of the limiting curves are barely visible even at the largest $\rho$ 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 $\rho$ 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 $\rho \to \infty$.

## 5.2 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

$$Z_{M,N}(z) = \sum_{i} \alpha_{i}(z) \Lambda_{i}(z)^{N}, \qquad (5.1)$$

where we shall refer to the $\Lambda_i(z)$ as eigenvalues, and the $\alpha_i(z)$ as the corresponding multiplicities. For a given $z$, let us order the eigenvalues by norm, so that $|\Lambda_1(z)| \ge |\Lambda_2(z)| \ge \dots$, and we call an eigenvalue $\Lambda_i(z)$ dominant (at $z$) if its norm is maximal, $|\Lambda_i(z)| \le |\Lambda_1(z)|$. Supposing a mild non-degeneracy condition, the BKW theorem then states that the accumulation set of zeros, as $N \to \infty$, 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., $|\Lambda_1(z_0)| > |\Lambda_2(z_0)|$) and the corresponding multiplicity vanishes (i.e., $\alpha_1(z_0) = 0$). A curve of accumulation points occurs when there are at least two dominant eigenvalues (i.e., $|\Lambda_1(z)| = |\Lambda_2(z)|$), and the relative phase $\phi(z) \in \mathbb{R}$ defined by $\Lambda_2(z) = e^{i\phi(z)}\Lambda_1(z)$ varies along the curve. The speed of variation of $\phi(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, $\alpha_i(z)$ and $\Lambda_i(z)$ depend on $M$ and $K$, and moreover $\alpha_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 \le 8$.

The resulting limiting curves for $M = 4, 5, \dots, 8$ are shown in figure 5. A number of qualitative features can be read off from these examples. First, the curves are invariant under the independent sign changes of Re $z$ and Im $z$. Second, they all contain the point
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Limiting curves of accumulation points of partition function zeros in the complex z-plane, shown in reading direction for $M = 4, 5, \dots, 8$.</figcaption></figure>

$z = i/2$. Third, they contain a number of branches extending to infinity; the number of such branches within each quadrant appears to be 3, 3, 5, 4, 7 for the sizes considered. Fourth, for even $M$ the curves do not intersect the real axis, while for odd $M$ they contain an exact vertical ray $z \in [-i/2, i/2]$. For odd $M$, there are further intersections with the real axis, namely $z \approx \pm3.5970$ for $M = 5$, as well as $z \approx \pm3.0096$ and $z \approx \pm6.0139$ for $M = 7$. Fifth, we only find T-points and now higher-order bifurcations.

A close scrutiny of figure 5 reveals that the limiting curves some tantalizingly tiny features for $M \ge 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.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left panel: partial results for the inner parts of the curves of accumulation points for odd $M = 5, 7, \dots, 17$. Right panel: extrapolation of the first real intersection, using a polynomial fit in the variable $1/M$.</figcaption></figure>

From a numerical point of view, the diagonalization computations become increasingly difficult as we approach the highly degenerate points $z = \pm i/2$. Practical details about the computational approach to limiting curves can be found in appendix D.

Comparing figures 4-5 gives convincing evidence that the partition function zeros indeed accumulate on the limiting curves, in the limit of large aspect ratio $\rho \to \infty$. 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 $\rho$ shown in figure 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 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 = \pm 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 6. The first intersection with the real axis appears to extrapolate to $z_* = 2$ as $M \to \infty$ (see the right panel of figure 6). Based on this, we conjecture that the enclosed region tends to a rectangle, given by $0 < \text{Re } z < 2$ and $-1/2 < \text{Im } z < 1/2$, in the thermodynamic limit, for $M$ odd.

## 6 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
<figure><img src="image_9.png" /><figcaption>Figure 7. Comparison between the limiting curves and the partition function zeros for $M = 7$.</figcaption></figure>

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 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 $\mathbb{Q}$ from some interesting observations about the partition function. We will see that it is useful to perform the decomposition
on a larger field $\mathbb{Q}(i, \xi_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].

## 6.1 Primary decomposition over $\mathbb{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):

$$ \begin{aligned} \Lambda_1 &= 2z^6 + \frac{5}{2}z^4 - \frac{\sqrt{54 - 6\sqrt{17}}}{2}z^3 + \frac{(2\sqrt{17} - 3)}{8}z^2 - \frac{\sqrt{54 - 6\sqrt{17}(3 + \sqrt{17})}}{16}z + \frac{(9 + 2\sqrt{17})}{32}, \\ \Lambda_2 &= 2z^6 + \frac{5}{2}z^4 + \frac{\sqrt{54 - 6\sqrt{17}}}{2}z^3 + \frac{(2\sqrt{17} - 3)}{8}z^2 + \frac{\sqrt{54 - 6\sqrt{17}(3 + \sqrt{17})}}{16}z + \frac{(9 + 2\sqrt{17})}{32}, \\ \Lambda_3 &= 2z^6 + \frac{5}{2}z^4 - \frac{\sqrt{54 + 6\sqrt{17}}}{2}z^3 - \frac{(2\sqrt{17} + 3)}{8}z^2 - \frac{\sqrt{54 + 6\sqrt{17}(\sqrt{17} - 3)}}{16}z + \frac{(9 - 2\sqrt{17})}{32}, \\ \Lambda_4 &= 2z^6 + \frac{5}{2}z^4 + \frac{\sqrt{54 + 6\sqrt{17}}}{2}z^3 - \frac{(2\sqrt{17} + 3)}{8}z^2 - \frac{\sqrt{54 + 6\sqrt{17}(\sqrt{17} - 3)}}{16}z + \frac{(9 - 2\sqrt{17})}{32}. \end{aligned} \quad (6.1) $$

One can check that although each $\Lambda_1, \dots, \Lambda_4$ is complicated and has irrational coefficients for generic $z$, their symmetric power sums

$$ \Lambda_1^n + \Lambda_2^n + \Lambda_3^n + \Lambda_4^n \quad (6.2) $$

are always polynomials whose coefficients are rational numbers! Since each $\Lambda_i$ corresponds to a solution of the BAE or the $TQ$-relation, this implies that we can group the four corresponding solutions of the BAE. Notice that we cannot make the decomposition further on $\mathbb{Q}$. If we further divide the four solutions into two groups, say $\Lambda_1, \Lambda_2$ and $\Lambda_3, \Lambda_4$, then the coefficients of the symmetric power sums $\Lambda_1^n + \Lambda_2^n$ and $\Lambda_3^n + \Lambda_4^n$ are no longer rational. This implies that these four solutions form an *irreducible* or *primary* block on $\mathbb{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,
<table><caption>**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.</caption><thead><tr><th>K</th><th>M = 6</th></tr></thead><tbody><tr><td>1</td><td>5 = 1+2+2</td></tr><tr><td>2</td><td>9 = 1+2+2+4</td></tr><tr><td>3</td><td>5 = 1+2+2</td></tr></tbody></table>

<table><caption>**Table 2.** The values of $l$ and eigenvalues of $U^3$ for physical solutions of BAE with $M = 6, K = 2$ in the four subspaces under primary decomposition. All the solutions in the same subspace have the same eigenvalues of $U^3$.</caption><thead><tr><th></th><th>dimensions</th><th>values of <i>l</i></th><th>eigenvalues of <i>U</i><sup>3</sup></th></tr></thead><tbody><tr><td>A</td><td>1</td><td>{3}</td><td>-1</td></tr><tr><td>B</td><td>2</td><td>{6, 6}</td><td>+1</td></tr><tr><td>C</td><td>2</td><td>{1, 5}</td><td>-1</td></tr><tr><td>D</td><td>4</td><td>{2, 2, 4, 4}</td><td>+1</td></tr></tbody></table>

with dimensions 1,2,2,4. In particular, the dimension-4 subspace corresponds to the four eigenvalues given in (6.1).

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

$$U_M = (-i)^M T_M(i/2). \qquad (6.3)$$

For a closed spin chain of length $M$, the allowed eigenvalues of the shift operator are

$$\exp\left(\frac{2\pi i l}{M}\right), \quad l = 1, \dots, M. \qquad (6.4)$$

Let us denote the four subspaces as A, B, C, D; we can then compute the values of $l$ for $M = 6$ for each subspace. The result is shown in table 2. We see that the value of $l$ is 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 $\mathbb{Q}$. Let us denote $\xi_M = \exp(2\pi i/M)$. It is clear that $\xi_M^l$ is not always rational for all $l = 1, \dots, M$. Therefore one cannot perform the decomposition over momentum sector completely on $\mathbb{Q}$. For each $M$, we can find the smallest integer $1 \le m \le M$ such that 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

$$A: \quad U = (-i)^6 t(i/2) = -1, \qquad (6.5)$$

$$B: \quad U = (-i)^6 t(i/2) = +1,$$

$$C: \quad U^3 = [(-i)^6 t(i/2)]^3 = -1, \quad U = (-i)^6 t(i/2) \neq -1,$$

$$D: \quad U^3 = [(-i)^6 t(i/2)]^3 = +1, \quad U = (-i)^6 t(i/2) \neq +1.$$
Notice that we need to include the constraints $U \neq \pm 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 $\xi_M^\ell$ is not always a rational number, we cannot perform the decomposition completely. However, we are not constrained to work on the field $\mathbb{Q}$. If we extend the field slightly to include $\xi_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 $\mathbb{F}_M = \mathbb{Q}(i, \xi_M)$ where $i$ is the imaginary unit.

After the decomposition into momentum sectors, we have $M$ subspaces<fnref n="2" /> (corresponding to $\ell = 1, \dots, 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.

## 6.2 Primary decomposition over $\mathbb{F}_M$

We explain in detail how to implement the decomposition over $\mathbb{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

$$ \prod_{j=1}^{K} \frac{u_j + i/2}{u_j - i/2} = \frac{Q_s(-i/2)}{Q_s(+i/2)} = \xi_M^\ell, \quad \ell = 1, \dots, M. \qquad (6.6) $$

Here $s = \{s_0, \dots, s_{K-1}\}$ are the unknown coefficients that we solve for in the rational $Q$-system. Alternatively, we can write $Q_s(z)$ as

$$ Q_s(z) = \prod_{j=1}^{K} (z - u_j), \qquad (6.7) $$

where $\{u_1, \dots, u_j\}$ are the Bethe roots. When the solution of BAE is singular, i.e., two of the Bethe roots are $\pm i/2$, we have $Q_s(\pm i/2) = 0$, whence the left-hand side of (6.6) is singular. This singularity can be eliminated by using the $TQ$-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(\pm i/2)$ is non-singular, we can write the momentum condition (6.6) in the following polynomial form

$$ P_{M,K,\ell} = Q_s(-i/2) - \xi_M^\ell Q_s(i/2) = 0. \qquad (6.8) $$

Consider the polynomial ring $\mathbb{A}_{M,K} = \mathbb{F}_M[s_0, \dots, s_{K-1}]$, and by abuse of notation we denote by $I_{M,K}$ the ideal in $\mathbb{A}_{M,K}$ generated by the rational $Q$-system. We define $M$ ideals in $\mathbb{A}_{M,K}$ as

$$ I_{M,K,\ell} = (I_{M,K} + \langle P_{M,K,\ell}, wQ_s(i/2) - 1 \rangle) \cap \mathbb{F}_M[s_0, \dots, s_{K-1}], \qquad (6.9) $$

<footnote n="2">²Some of the subspaces might not exist for some values of $M$ and $K$, as witnessed by the tables in the next subsection.</footnote>
where $\ell = 1, \dots, M$, and “+” means the sum of two ideals. Here $w$ is an auxiliary variable to remove the singular Bethe roots.

We further define an additional ideal for the singular case,

$$I_{M,K,\infty} \equiv I_{M,K} + \langle Q_s(i/2) \rangle. \qquad (6.10)$$

Let $\mathcal{Z}(\dots)$ be the common solution of a set of equations, or equivalently the algebraic set of the corresponding ideal, in the algebraic closure $\tilde{\mathbb{Q}}$ of rational numbers. We claim that,

$$\mathcal{Z}(I_{M,K}) = \left( \bigcup_{\ell=1}^{M} \mathcal{Z}(I_{M,K,\ell}) \right) \bigcup \mathcal{Z}(I_{M,K,\infty}). \qquad (6.11)$$

and thus the BAE roots are classified into $M+1$ subsets.

From the construction of these ideals, we see that for any point $\mathbf{x} \in \mathcal{Z}(I_{M,K,\ell})$,

$$Q_{\mathbf{x}}(i/2) \neq 0, \qquad \frac{Q_{\mathbf{x}}(-i/2)}{Q_{\mathbf{x}}(+i/2)} = \xi_M^\ell \qquad (6.12)$$

and it is clear that for $\mathbf{x} \in \mathcal{Z}(I_{M,K,\infty})$,

$$Q_{\mathbf{x}}(i/2) = 0. \qquad (6.13)$$

Hence, for different $\ell \in \{1, \dots, M, \infty\}$, the algebraic sets $\mathcal{Z}(I_{M,K,\ell})$ have no intersection, so the union (6.11) is disjoint. Since the $Q$-system equation has no coinciding Bethe roots by construction, $I_{M,K,\ell}$ are all radical ideals. By Hilbert's Nullstellensatz,

$$I_{M,K} = \left( \bigcap_{\ell=1}^{M} I_{M,K,\ell} \right) \bigcap I_{M,K,\infty}. \qquad (6.14)$$

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 $\ell \in \{1, \dots, M, \infty\}$, each $I_{M,K,\ell}$ is primary over the field $\mathbb{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),

$$Z_{M,N}(z) = \sum_{K=0}^{[M/2]} (M - 2K + 1) \left( \sum_{\ell \in \{1, \dots, M, \infty\}} \mathrm{Tr} (\mathbf{T}_{M,K,\ell}(z)^N) \right), \qquad (6.15)$$

where the companion matrix $\mathbf{T}_{M,K,\ell}(z)$ is the companion matrix for

$$(a(z)Q_{M,K}(z-i) + d(z)Q_{M,K}(z+i))Q_{M,K}(z)^{-1} \qquad (6.16)$$
in the ideal $I_{M,K,l}$. Note that $T_{M,K,l}(z)$ contains polynomials in $i$ and $\xi_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,\infty}$, i.e., singular roots. The singularity of the eigenvalue of $U$ in terms of $Q_s(\pm i/2)$ is actually spurious. They can be eliminated by using the $TQ$-relation. We can combine the equations from the rational $Q$-system and the $TQ$-relation and then eliminate the variables $s$. The elimination procedure is actually quite simple, due to the structure of equations from $TQ$-relations. This gives us a set of equations that only involve the variables $t = \{t_0, t_1, \dots, t_M\}$. The momentum conditions in terms of $t$ variables are simply

$$(-i)^M t_t(i/2) = \xi_M^l, \quad l = 1, 2, \dots, M \qquad (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,\infty}$.

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 by-product 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,l}$ are related by Galois group actions; therefore, instead of exhaustively going through the sum over all $l$'s in (6.15), we just need to compute a few $l$'s, namely one for each orbit.

Note that for two different decomposed BAE with $1 \le l_1, l_2 \le M$ described in the previous subsection, if we replace $\xi_M^{l_1}$ by $\xi_M^{l_2}$ in the generators from $I_{M,K,l_1}$ (6.9), without changing the imaginary unit $i$ or any rational coefficient, then the ideal $I_{M,K,l_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 = \text{Gal}(F_M/\mathbb{Q}(i))$ where

$$\mathbb{Q}(i) = \{a + bi : a, b \in \mathbb{Q}\}. \qquad (6.18)$$

Specifically, if there exists an element $g \in G$ such that,

$$g(\xi_M^{l_1}) = \xi_M^{l_2}, \qquad (6.19)$$

then by the field automorphism of Gröbner basis computation,

$$g(\mathbf{T}_{M,K,l_1}(z)) = \mathbf{T}_{M,K,l_2}(z), \qquad (6.20)$$
so the computation for $\mathbf{T}_{M,K,\ell_2}$ is no longer needed. Here $g$ acts on each element of the companion matrix. Then, instead of taking the sum over $\ell = 1, \dots, M$, we just need to find the orbits, under the $G$-action, of the set $\{\xi_M, \xi_M^2, \dots, 1\}$ and compute only one companion matrix for each orbit.

Hence it is important to analyze the structure of the algebraic extension $[\mathbb{F}_M : \mathbb{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 $\mathbb{F}_M$ is the cyclotomic field,

$$ \mathbb{F}_M = \mathbb{Q}(\xi_M, i) = \mathbb{Q}(\xi_{4M}) = \mathbb{Q}\left(e^{\frac{2\pi i}{4M}}\right). \qquad (6.21) $$

Note that $i = \xi_{4M}^M$. The large Galois group, $\mathrm{Gal}(\mathbb{F}_M/\mathbb{Q})$, is the multiplication group $(\mathbb{Z}/(4M)\mathbb{Z})^\times$ with the size $\phi(4M) = \phi(4)\phi(M) = 2\phi(M)$. Here $\phi(\cdot)$ is the Euler totient function. We have

$$ |G| = |\mathrm{Gal}(\mathbb{F}_M/\mathbb{Q}(i))| = \phi(M) \qquad (6.22) $$

From elementary number theory, there exists a $g \in G$ such that,

$$ g(\xi_M^{\ell_1}) = \xi_M^{\ell_2}, \quad g(i) = i. \qquad (6.23) $$

Hence $I_{M,K,\ell_1}$ and $I_{M,K,\ell_2}$ are equivalent if and only if $\mathrm{gcd}(\ell_1, M) = \mathrm{gcd}(\ell_2, M)$. We conclude that in this case, the ideals (6.9) are classified by the greatest common divisors. Under the Galois group $G$, the number of orbits is $\sigma_0(M)$, where $\sigma_0(\dots)$ denotes the divisor function which counts the number of divisors of $M$. Furthermore, with the ideal $I_{M,K,\infty}$, we need to compute $\sigma_0(M) + 1$ companion matrices when $M$ is odd.

2. $M$ is even and $4 \nmid M$. In this case, the field $\mathbb{F}_M$ is the cyclotomic field

$$ \mathbb{F}_M = \mathbb{Q}(\xi_M, i) = \mathbb{Q}(\xi_{2M}) = \mathbb{Q}\left(e^{\frac{2\pi i}{2M}}\right). \qquad (6.24) $$

$\mathrm{Gal}(\mathbb{F}_M/\mathbb{Q})$ is the multiplication group $(\mathbb{Z}/(2M)\mathbb{Z})^\times$ with the size $\phi(2M) = \phi(4)\phi(M/2) = 2\phi(M)$. We have

$$ |\mathrm{Gal}(\mathbb{F}_M/\mathbb{Q}(i))| = \phi(M). \qquad (6.25) $$

As in the previous case, $I_{M,K,\ell_1}$ and $I_{M,K,\ell_2}$ are equivalent if and only if $\mathrm{gcd}(\ell_1, M) = \mathrm{gcd}(\ell_2, M)$. We need to compute $\sigma_0(M) + 1$ companion matrices, as in the preceding case.

3. $4|M$. This case is different from the previous ones. The field $\mathbb{F}_M$ is the cyclotomic field

$$ \mathbb{F}_M = \mathbb{Q}(\xi_M, i) = \mathbb{Q}(\xi_M) = \mathbb{Q}\left(e^{\frac{2\pi i}{M}}\right). \qquad (6.26) $$

Note that $i$ is a power of $\xi_M$. $G = \mathrm{Gal}(\mathbb{F}_M/\mathbb{Q}(i))$ is the subgroup of $\mathrm{Gal}(\mathbb{F}_M/\mathbb{Q})$ which keeps $i$ invariant.
The classification of $\ell$'s is more complicated in this case, since $i \in Q(\xi_M)$. From detailed Galois theory analysis, $I_{M,K,\ell_1}$ and $I_{M,K,\ell_2}$ are equivalent if and only if

$$ \gcd(M, \ell_1) = \gcd(M, \ell_2), \quad \frac{\ell_2}{\ell_1} = 1 \quad \mod \gcd\left(\frac{M}{\gcd(M, \ell_1)}, 4\right). \qquad (6.27) $$

Note that the denominator of the reduced fraction $\ell_2/\ell_1$ is relatively prime to $M/\gcd(M, \ell_1)$, so the congruence condition for $\ell_2/\ell_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,

$$ s_i \mapsto (-1)^{k-i}s_i, \quad i = 0, \dots, K. \qquad (6.28) $$

Under this transformation, we find that up to $M \le 16$ the ideal $I_{M,K}$ is invariant. Furthermore, the polynomial for the momentum condition transforms as

$$ P_{M,K,\ell}(i) \mapsto (-1)^{k}P_{M,K,\ell}(-i). \qquad (6.29) $$

This means that the imaginary unit $i$ in $P_{M,K,\ell}$ transforms to $-i$. Hence, for two integers $\ell_1$ and $\ell_2$, such that $\gcd(M, \ell_1) = \gcd(M, \ell_2)$ but which do not satisfy (6.27), when the enhanced symmetry (6.28) exists, $I_{M,K,\ell_1}$ and $I_{M,K,\ell_2}$ are still equivalent.

In summary, with the enhanced symmetry (6.28), for any positive integer $M$, $I_{M,K,\ell_1}$ and $I_{M,K,\ell_2}$ are equivalent if and only if $\gcd(M, \ell_1) = \gcd(M, \ell_2)$. Hence there are always $\sigma_0(M) + 1$ classes.

## 6.3 Results on decomposition over $\mathbb{F}_M$

In this subsection, we give some results of BAE decomposition over the extended field $\mathbb{F}_M$.

**M = 6.** In this case $\mathbb{F}_6 = \mathbb{Q}(e^{\frac{2\pi i}{12}})$. From the discussion of Galois theory in the previous section, the decomposed BAE sub-systems are classified by the value of $\ell$:

$$ \{1, 5\}, \{2, 4\}, \{3\}, \{6\}, \{\infty\}. \qquad (6.30) $$

Two sub-systems, whose $\ell$-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,\ell}$ with $\ell = 1, 2, 3, 4, 5, 6, \infty$, and get the Bethe root counting for the decomposed BAE with $M = 6$ in table 3. For the singular Bethe roots, via the $TQ$-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 $\ell = 3$ means that $I_{6,2,3}$ has no solution but there is one singular Bethe root whose regularized momentum value is $e^{3 \times 2\pi i/6} = e^{\pi i}$.

**M = 7.** In this case $\mathbb{F}_7 = \mathbb{Q}(e^{\frac{2\pi i}{28}})$. The decomposed BAE sub-systems are classified by the value of $\ell$:

$$ \{1, 2, 3, 4, 5, 6\}, \{7\}, \{\infty\}. \qquad (6.31) $$

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.
<table><caption>**Table 3.** Number of Bethe roots for decomposed BAE with $M = 6$. There are $\sigma_0(6) + 1 = 5$ classes, which correspond to $\{1, 5\}, \{2, 4\}, \{3\}, \{6\}, \{\infty\}$. The numbers in brackets are the singular Bethe roots.</caption><thead><tr><th></th><th>{1, 5}</th><th>{2, 4}</th><th>{3}</th><th>{6}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>1</td><td>2</td><td>0 (1)</td><td>2</td><td>1</td></tr><tr><th>K = 3</th><td>1</td><td>0</td><td>2</td><td>0 (1)</td><td>1</td></tr></tbody></table>

<table><caption>**Table 4.** Number of Bethe roots for decomposed BAE with $M = 7$.</caption><thead><tr><th></th><th>{1, 2, 3, 4, 5, 6}</th><th>{7}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>2</td><td>2</td><td>0</td></tr><tr><th>K = 3</th><td>2</td><td>2</td><td>0</td></tr></tbody></table>

<table><caption>**Table 5.** Number of Bethe roots for decomposed BAE with $M = 8$.</caption><thead><tr><th></th><th>{1, 3, 5, 7}</th><th>{2, 6}</th><th>{4}</th><th>{8}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>2</td><td>3</td><td>2 (1)</td><td>3</td><td>1</td></tr><tr><th>K = 3</th><td>4</td><td>3</td><td>3</td><td>2 (1)</td><td>1</td></tr><tr><th>K = 4</th><td>1</td><td>2</td><td>0 (3)</td><td>3</td><td>3</td></tr></tbody></table>

**M = 8.** In this case $\mathbb{F}_8 = \mathbb{Q}(e^{\frac{2\pi i}{8}})$. The decomposed BAE sub-systems are classified by the value of $\ell$, via the condition (6.27):

$$ \{1, 5\}, \{2\}, \{3, 7\}, \{4\}, \{6\}, \{8\}, \{\infty\}. \qquad (6.32) $$

However, due to the enhanced symmetry (6.28), classes like $\{1, 5\}$ and $\{3, 7\}$ are equivalent, and similarly $\{2\}$ and $\{6\}$ are also equivalent. Therefore we can use the new classification,

$$ \{1, 3, 5, 7\}, \{2, 6\}, \{4\}, \{8\}, \{\infty\} \qquad (6.33) $$

The root counting for the decomposed BAE with $M = 8$ is given in table 5. Note that for the cases $M \le 8$, each BAE subsystem contains at most four roots. Hence the Bethe roots form a *solvable algebraic extension* over $\mathbb{F}_M$ by Galois theory, which means that each Bethe root can be expressed as radicals of rational numbers for $M \le 8$. The closed-form expressions for the partition function with $M = 7, 8$ will be presented in the future work.

Similar results for $9 \le M \le 18$ are presented in appendix E (see tables 6–15). We use the Gröbner basis method with coefficients in the algebraic extension $\mathbb{F}_M$ to get the Bethe roots classification. Our code is powered by SINGULAR [16].

## 6.4 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 $\ell$. 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 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$.

With $K = 0$ there is just one state, $|000000\rangle$, so we set $g_0 = 1$. With $K = 1$ there is one orbit with $g_1 = 6$ having the representative state $|100000\rangle$. With $K = 2$ there are three orbits with representative states

$$|110000\rangle (g_2^A = 6), \quad |101000\rangle (g_2^B = 6), \quad |100100\rangle (g_2^C = 3). \qquad (6.34)$$

And finally, with $M = 3$ there are four orbits with representative states

$$|111000\rangle (g_3^A = 6), \quad |110100\rangle (g_3^B = 6), \quad |101100\rangle (g_3^C = 6), \quad |101010\rangle (g_3^D = 2). \qquad (6.35)$$

Let the momentum label $\ell = 1, 2, \dots, M$ be defined as before. In general, a given orbit $s$ of length $g_s$ is compatible with $\ell$ if and only if

$$\ell g_s = 0 \mod M. \qquad (6.36)$$

Since obviously $g_s|M$, it is not hard to see from (6.36) that the number of compatible orbits must be constant on each class identified above on by Galois theory, except that we do not have the $\ell = \infty$ class in the present case. The classes for $M = 6$ then follow from (6.30), namely $\{1, 5\}, \{2, 4\}, \{3\}, \{6\}$. In table form the numbers of compatible orbits now read:

<table><thead><tr><th></th><th>{1, 5}</th><th>{2, 4}</th><th>{3}</th><th>{6}</th></tr></thead><tbody><tr><th>K = 0</th><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>K = 2</th><td>2</td><td>3</td><td>2</td><td>3</td></tr><tr><th>K = 3</th><td>3</td><td>3</td><td>4</td><td>4</td></tr></tbody></table>

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} = \binom{M}{K}$. One may now further block-diagonalize the $T_{M,K}(z)$ into momentum sectors $T_{M,K,\ell}(z)$, having the dimensions $d_{M,K,\ell}$ given by the above table, by using a procedure similar to appendix A.4 of [24]. To this end we write

$$T_{M,K,\ell}(z) = S_{\text{out}} T_{M,K}(z) S_{\text{in}}. \qquad (6.37)$$

Here $S_{\text{in}}$ is a $d_{M,K,\ell} \times d_{M,K}$ matrix that maps each compatible orbit into its representative state, with weight $g_s$. And $S_{\text{out}}$ is a $d_{M,K} \times d_{M,K,\ell}$ matrix that maps each state into a representative (and hence into an orbit), and attributes a weight $(\xi_M)^{\ell \cdot 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 $\xi_M = \exp(2\pi i/M)$, as before. One may now verify by explicit construction of
the matrices $T_{M,K,l}(z)$ that the spectrum of $T_{M,K}(z)$ is indeed the union of the spectra of the momentum blocks $T_{M,K,l}(z)$.

As it stands, this method does not yet take into account the $\mathfrak{su}(2)$ symmetry of the XXX spin chain. This means that each $T_{M,K,l}(z)$ contains all the highest-weight states with equal or higher spin (i.e., $K' \le 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:

<table><thead><tr><th></th><th>{1, 5}</th><th>{2, 4}</th><th>{3}</th><th>{6}</th></tr></thead><tbody><tr><th>$K = 0$</th><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><th>$K = 1$</th><td>1</td><td>1</td><td>1</td><td>0</td></tr><tr><th>$K = 2$</th><td>1</td><td>2</td><td>1</td><td>2</td></tr><tr><th>$K = 3$</th><td>1</td><td>0</td><td>2</td><td>1</td></tr></tbody></table>

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, \ell)$ sectors. Notice that the present approach does not particularize the singular case denoted $\{\infty\}$ in table 3, but assigns to it straight away the correct regularized momentum, namely $\{\infty\} \to \{3\}$ for $K = 2$, and $\{\infty\} \to \{6\}$ for $K = 3$.<fnref n="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, \ell)$. 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,\ell}$ is the sum over those $N(g_s; M, K)$ that respect the criterion (6.36), and the highest-weight combinations $d_{M,K,\ell}-d_{M,K-1,\ell}$ provide precisely the numbers of tables 3–15, up to the assignment of a definite momentum to the $\{\infty\}$ classes.

We have written a simple algorithm that carries out this computation. It produces the tables for $M \le 24$ in less than one minute. For $M \le 18$ these are in full agreement with tables 3–15, after assigning to each $\{\infty\}$ case the corresponding regularized momentum. We note that for $M$ even, this assignment appears to obey a simple rule: $\{\infty\} \to \{M/2\}$ when $K$ is even, and $\{\infty\} \to \{M\}$ when $K$ is odd. When $M$ is prime, the decomposition of the BAE is very simple, and there is no $\{\infty\}$. It remains to discuss the cases of odd non-prime $M \le 18$, namely $M = 9$ and $M = 15$. For $M = 9$ we find that $\{\infty\} \to \{3, 6\}$ when $K = 3$ (see table 6), and for $M = 15$ we have $\{\infty\} \to \{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 $TQ$-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 $\mathfrak{su}(2)$ symmetry of the XXX

<footnote n="3">³The first of these identifications agrees as well with table 2 for the subspace called A there.</footnote>
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.

# 7 Conclusions and discussions

In this paper we developed a method to compute the torus partition function of the six-vertex 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, algebraic extension and Galois theory.

Using this approach, we probed new structures in the solution space of the BAE. We found that the solution space can be decomposed into subspaces $I_{M,K,l}$ on an algebraically extended field $\mathbb{F}_M = \mathbb{Q}(i, \xi_M)$, where $M, K$ are the length, magnon number, and $l = 1, \dots, M$ is related to the total momentum $2\pi l/M$ of the solutions of the BAE. We classified the subspaces that are related by the action of the Galois group and form an orbit. For each orbit, we computed the Gröbner basis, quotient ring and the companion matrices of the transfer matrix $T_{M,K,l}(z)$ and Baxter's $Q$-operator $Q_{M,K,l}(z)$ up to $M=14$. This decomposition can also be understood in terms of a naive momentum-sector decomposition described in section 6.4. The decomposition of the solution space into $I_{M,K,l}$ and the use of Galois theory led to a huge boost in the computation of the partition function. For $M \le 6$, we have closed-form expressions for the partition function for any $N$. For larger $M$ up to $M=14$, the partition functions for fixed $N$ can be computed straightforwardly from the companion matrices $T_{M,K,l}(z)$.

The exact partition functions are polynomials in the spectral parameter $z$ of order $MN$ 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 \gg 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 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,<fnref n="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 \neq 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

<footnote n="4">⁴Notice that these have previously been studied by more combinatorial methods in the special case of the chromatic polynomial [7].</footnote>
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 algebro-geometric 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(\mathfrak{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.

## Acknowledgments

Y. Jiang is partially supported by the Swiss National Science Foundation through the NCCR SwissMap. J.L. Jacobsen acknowledges support from the European Research Council through the advanced grant NuQFT. Y. Zhang is supported by funding from Swiss National Science Foundation (Ambizione grant PZ00P2 161341) and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 725110). We thank Zoltan Bajnok, Janko Böhm, Matthias Staudacher and Yuwei Zhang for helpful discussions.

## A Rational $Q$-system

In this appendix, we review the rational $Q$-system method for solving Bethe ansatz equations (BAE) proposed by Marboe and Volin [11]. One of the main advantages of this method is that unlike the BAE, the solutions of the rational $Q$-system are all physical. In addition, it is by far the most efficient way of finding Bethe roots for fixed length $M$ and particle number $K$.

Let us briefly review how the rational $Q$-system works. For more details, we refer to the original papers [11, 12]. The method works for a large class of rational spin chains with
su $(m, n|k)$ symmetry. For our purpose, we restrict to the su (2) case. The main procedure is as follows:

1. For a BAE of length $M$ and magnon number $K$, we draw a Young diagram of two rows $(M - K, K)$ with $M - K \ge K$.<fnref n="5" />

2. To each node (i.e., a corner of a box of the Young diagram) we associate a Q-function and require that all the Q-functions are polynomials. The order of a given Q-polynomial is given by the number of boxes on the upper-right part of the corresponding node (see figure 8 for an illustration).

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

$$Q_{a+1,s}(u)Q_{a,s+1}(u) = Q_{a+1,s+1}^{+}(u)Q_{a,s}^{-}(u) - Q_{a+1,s+1}^{-}(u)Q_{a,s}^{+}(u), \quad (A.1)$$

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

$$Q_{a,s}^{\pm}(u) = Q_{a,s}(u \pm i/2), \quad Q_{a,s}^{++}(u) = Q_{a,s}(u + i), \quad Q_{a,s}^{--}(u) = Q_{a,s}(u - i). \quad (A.2)$$

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.

5. 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 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

<footnote n="5">⁵Since the solution of BAE gives the same eigenstate as the corresponding dual solution, this restriction already covers the whole Hilbert space.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The Young diagram for $M = 6$, $K = 2$. The number on each node denotes the degree of the $Q$-polynomial.</figcaption></figure>

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. \qquad (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 8. Apart from the two boundary $Q$-polynomials, there are two other unknown $Q$-polynomials which are parameterized by

$$Q_{0,1}(u) = u^2 + c_{0,1}^{(1)} u + c_{0,1}^{(0)}, \quad Q_{1,1} = u + c_{1,1}^{(0)}. \qquad (A.4)$$

The remaining task is to use the $QQ$-relations (A.1) to determine the remaining unknown $Q$-polynomials as well as the coefficients $c_{0,1}^{(1)}$, $c_{0,1}^{(0)}$ and $c_{1,1}^{(0)}$. In fact, in our case there are only three remaining non-trivial functions to be determined, namely $Q_{1,0}$, $Q_{2,0}$ and $Q_{3,0}$.

Let us first determine $Q_{1,0}(u)$. Taking $a = s = 0$, the $QQ$-relation leads to

$$Q_{1,0}(u) = c \frac{Q_{1,1}^{+}(u)Q_{0,0}^{-}(u) - Q_{1,1}^{-}(u)Q_{0,0}^{+}(u)}{Q_{0,1}(u)} \qquad (A.5)$$

where $c$ is some normalization constant to make the polynomial monic. The quotient and the remainder can be computed straightforwardly. Plugging (A.3) and (A.4) into (A.5), we find a non-trivial remainder in terms of the unknown coefficients $c_{0,1}^{(1)}$, $c_{0,1}^{(0)}$ and $c_{1,1}^{(0)}$. 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.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The Young tableaux with two rows. The first row has $K$ boxes and the second row has $M-K$ boxes.</figcaption></figure>

To this end, let us consider a generic Young tableaux with two rows of $(M-K, K)$ boxes, as shown in figure 9. We have $K \le 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, \dots, u_K$ and write

$$Q_{0,1}(u) = \prod_{j=1}^{K} (u - u_j). \qquad (A.6)$$

We consider the $QQ$-relations in the two shaded boxes, which read

$$Q_{1,1}Q_{0,2} = Q_{1,2}^{+}Q_{0,1}^{-} - Q_{1,2}^{-}Q_{0,1}^{+}, \qquad (A.7)$$

$$Q_{1,0}Q_{0,1} = Q_{1,1}^{+}Q_{0,0}^{-} - Q_{1,1}^{-}Q_{0,0}^{+}. \qquad (A.8)$$

At $u = u_k$, $Q_{0,1}(u_k) = 0$, and by shifting (A.7) by $\pm i/2$ we obtain the following relations:

$$Q_{1,1}^{+}(u_k)Q_{0,2}^{+}(u_k) = -Q_{1,2}(u_k)Q_{0,1}^{++}(u_k), \quad Q_{1,1}^{-}(u_k)Q_{0,2}^{-}(u_k) = Q_{1,2}(u_k)Q_{0,1}^{--}(u_k). \qquad (A.9)$$

Using the boundary condition $Q_{0,2} = Q_{1,2} = 1$, we obtain

$$Q_{1,1}^{+}(u_k) = -Q_{0,1}^{++}(u_k), \quad Q_{1,1}^{-}(u_k) = Q_{0,1}^{--}(u_k). \qquad (A.10)$$

Meanwhile, evaluating (A.8) at $u = u_k$ produces

$$Q_{1,1}^{+}(u_k)Q_{0,0}^{-}(u_k) - Q_{1,1}^{-}(u_k)Q_{0,0}^{+}(u_k) = 0, \qquad (A.11)$$

and plugging (A.10) into (A.11) we obtain

$$Q_{0,1}^{++}(u_k)Q_{0,0}^{-}(u_k) + Q_{0,1}^{--}(u_k)Q_{0,0}^{+}(u_k) = 0. \qquad (A.12)$$

Use the boundary condition and the expression for $Q_{0,1}(u)$,

$$Q_{0,0}(u) = u^M, \quad Q_{0,1}(u) = \prod_{j=1}^{K} (u - u_j), \qquad (A.13)$$

we finally arrive at

$$\left( \frac{u_k + i/2}{u_k - i/2} \right)^M = - \prod_{j=1}^{K} \frac{u_k - u_j + i}{u_j - u_k - i}, \qquad (A.14)$$
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 $\mathbb{F}$ be a field and $A = \mathbb{F}[z_1, \dots, z_n]$ be a polynomial. Any ideal $I$ in $A$ decomposes into the intersection of several primary ideals.

$$I = I_1 \cap I_2 \cap \dots \cap I_k, \qquad (B.1)$$

where each $I_j$ is primary. A primary ideal $J$ is an ideal such that if $a, b \in A$, $ab \in J$, then either $a \in J$ or $b^n \in 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

$$\mathcal{Z}(I) = \mathcal{Z}(I_1) \cup \mathcal{Z}(I_2) \cup \dots \cup \mathcal{Z}(I_k). \qquad (B.2)$$

This property is useful for the study of the complicated zero set $\mathcal{Z}(I)$. Especially, when $\mathcal{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 $\mathbb{F}$. For example, within $\mathbb{Q}[x, y]$,

$$\langle x^2 - 2y^2 \rangle \qquad (B.3)$$

is already primary and hence cannot be decomposed further. However, within the algebraic closure $\overline{\mathbb{Q}}[x, y]$, we have the primary decomposition:

$$\langle x^2 - 2y^2 \rangle = \langle x - y\sqrt{2} \rangle \cap \langle x + y\sqrt{2} \rangle. \qquad (B.4)$$

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

$$Z_{M,N} = \sum_{K=0}^{[M/2]} (M - 2K + 1) \mathrm{Tr} (\mathbf{T}_{M,K}(z)^N), \qquad (C.1)$$
where $\mathbf{T}_{M,K}(z)$ is the companion matrix of the polynomial $t_{M,K}(\mathbf{s}, z)$. Recall also that $\mathbf{T}_{M,K}(z)$, calculated from the Gröbner basis, contains only rational numbers. Therefore the whole computation is manifestly analytic.

In practice, although $\mathbf{T}_{M,K}(z)$ can be calculated from the straightforward Gröbner basis and polynomial division procedure, the matrix power $(\mathbf{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})$.

2. Divide $t_{M,K}(\mathbf{s}, z)$ towards $G(I_{M,K})$; the remainder $\overline{t_{M,K}(\mathbf{s}, z)}$ is called the normal form.

3. Assume $N$ is a power of 2, i.e., $N = 2^c$. Recursively compute $t_{M,K}^{(i)}(\mathbf{s}, z) \equiv (\overline{t_{M,K}^{(i-1)}(\mathbf{s}, z)})^2$, for $i = 2, \dots, c+1$, where $t_{M,K}^{(1)}(\mathbf{s}, z) = \overline{t_{M,K}^{(1)}(\mathbf{s}, z)}$. Divide $t_{M,K}^{(i)}(\mathbf{s}, z)$ towards $G(I_{M,K})$ and define $\overline{t_{M,K}^{(i)}(\mathbf{s}, z)}$ to be the remainder.

4. Calculate the companion matrix of $\overline{t_{M,K}^{(c+1)}(\mathbf{s}, z)}$. This is the demanded companion matrix $\mathbf{T}_{M,K}(z)^N$.

This strategy avoids the storage and multiplication of dense matrices. Here we used the algebraic geometry property (3.9) that the companion matrix of a product of two polynomials, equals the product of the corresponding two companion matrices. If $N$ is not a power of 2, we find the binary representation of $N$ which is $N = 2^{c_1} + 2^{c_2} + \dots$ and repeat the computation several times.

In this algorithm, we trim each intermediate polynomial $\overline{t_{M,K}^{(i)}(\mathbf{s}, z)}$ towards $G(I_{M,K})$ via the polynomial division, to get a much short polynomial $\overline{t_{M,K}^{(i)}(\mathbf{s}, z)}$. This significantly saves the RAM usage.

This algorithm can also be combined with the decomposition described in section 6. Instead of calculating $G(I_{M,K})$, we calculate $G(I_{M,K,l})$ over the algebraic extension $\mathbb{F}_M$. Since $\mathcal{Z}(I_{M,K,l})$ contains much fewer points than those in $\mathcal{Z}(I_{M,K})$, the quotient ring $\mathbb{F}_M[s_0, \dots, s_{K-1}]/I_{M,K,l}$ has much lower dimension. Hence the polynomial division step trims the polynomial size more dramatically and further speeds up the computation. By (6.20), we just need to compute $\sigma_0(M) + 1$ distinct $\mathbf{T}_{M,K,l}(z)^N$, and get the rest by the Galois group action.

Our algorithm is powered by the SINGULAR [16] code. For the application with decomposed BAE, we introduce $\xi_M$ through the `minpoly` command in SINGULAR with the explicit minimal polynomial for $\xi_M$ (cyclotomic polynomial).

## D Details on finding limiting curves

We here describe the numerical procedure employed in section 5.2 to find the limiting curves along which the partition function zeros $Z_{M,N}(z) = 0$ accumulate, in the limit $N \to \infty$, for a fixed value of $M$.
The key ingredient is obviously to be able to efficiently diagonalize the transfer matrix, so that the loci of equimodularity ($|\Lambda_1(z)| = |\Lambda_2(z)|$) can be identified. To this end, we do not need the entire spectrum of $T_M(z)$ but only the first few eigenvalues (in principle, just the first two, but see the remarks below). The most efficient means of finding those are iterative Krylov-subspace methods that depend only on implementing the multiplication of the transfer matrix with a vector, i.e., to compute $v' = T_M(z)v$ when given some vector $v$. Our method of choice is the Arnoldi algorithm, for which we use the ARPACK implementation of Arnoldi's algorithm for complex matrices [32].<fnref n="6" /> This implementation is 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 $2Md$ operations, where $d$ is the dimension of the matrix. The trace $\text{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, \ell)$ for the equimodular eigenvalues $\Lambda_1(z)$ and $\Lambda_2(z)$ were known beforehand, it would obviously be most efficient to diagonalize the smaller matrices $T_{M,K,\ell}(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, \dots, \lfloor M/2 \rfloor$ simultaneously and with decomposing the momentum with respect to the $\ell$ label. Note that the $su(2)$ highest-weight constraint is also not enforced in this approach, so the eigenvalues with $K < \lfloor M/2 \rfloor$ present degeneracies. We deal with this in practice by imposing the equimodularity criterion $|\Lambda_1(z)| = |\Lambda_\ell(z)|$, where $\ell \ge 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 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.<fnref n="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

<footnote n="6">⁶We thank C.R. Scullard for technical discussions about the use of ARPACK.</footnote>

<footnote n="7">⁷We chose $|\Delta z| = 10^{-2}$ in most cases, but decreased to $10^{-3}$ close to T-points and other fine details.</footnote>
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 present only a single connected component gives appealing evidence that they are actually complete.

## E More results of primary decomposition on $\mathbb{F}_M$

In this appendix, we list the results of primary decomposition on $\mathbb{F}_M$ for $9 \le M \le 18$.

<table><caption>**Table 6.** Number of Bethe roots for decomposed BAE with $M = 9$.</caption><thead><tr><th></th><th>{1, 2, 4, 5, 7, 8}</th><th>{3, 6}</th><th>{9}</th><th>{∞}</th></tr></thead><tbody><tr><th>$K = 1$</th><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>$K = 2$</th><td>3</td><td>3</td><td>3</td><td>0</td></tr><tr><th>$K = 3$</th><td>5</td><td>5 (1)</td><td>6</td><td>2</td></tr><tr><th>$K = 4$</th><td>5</td><td>4</td><td>4</td><td>0</td></tr></tbody></table>

<table><caption>**Table 7.** Number of Bethe roots for decomposed BAE with $M = 10$.</caption><thead><tr><th></th><th>{1, 3, 7, 9}</th><th>{2, 4, 6, 8}</th><th>{5}</th><th>{10}</th><th>{∞}</th></tr></thead><tbody><tr><th>$K = 1$</th><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>$K = 2$</th><td>3</td><td>4</td><td>2 (1)</td><td>4</td><td>1</td></tr><tr><th>$K = 3$</th><td>8</td><td>7</td><td>8</td><td>6 (1)</td><td>1</td></tr><tr><th>$K = 4$</th><td>8</td><td>10</td><td>4 (4)</td><td>10</td><td>4</td></tr><tr><th>$K = 5$</th><td>5</td><td>3</td><td>6</td><td>0 (4)</td><td>4</td></tr></tbody></table>

<table><caption>**Table 8.** Number of Bethe roots for decomposed BAE with $M = 11$.</caption><thead><tr><th></th><th>{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}</th><th>{11}</th><th>{∞}</th></tr></thead><tbody><tr><th>$K = 1$</th><td>1</td><td>0</td><td>0</td></tr><tr><th>$K = 2$</th><td>4</td><td>4</td><td>0</td></tr><tr><th>$K = 3$</th><td>10</td><td>10</td><td>0</td></tr><tr><th>$K = 4$</th><td>15</td><td>15</td><td>0</td></tr><tr><th>$K = 5$</th><td>12</td><td>12</td><td>0</td></tr></tbody></table>

<table><caption>**Table 9.** Number of Bethe roots for decomposed BAE with $M = 12$.</caption><thead><tr><th></th><th>{1, 5, 7, 11}</th><th>{2, 10}</th><th>{3, 9}</th><th>{4, 8}</th><th>{6}</th><th>{12}</th><th>{∞}</th></tr></thead><tbody><tr><th>$K = 1$</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>$K = 2$</th><td>4</td><td>5</td><td>4</td><td>5</td><td>4 (1)</td><td>5</td><td>1</td></tr><tr><th>$K = 3$</th><td>13</td><td>12</td><td>14</td><td>12</td><td>13</td><td>12 (1)</td><td>1</td></tr><tr><th>$K = 4$</th><td>22</td><td>24</td><td>21</td><td>25</td><td>18 (5)</td><td>24</td><td>5</td></tr><tr><th>$K = 5$</th><td>26</td><td>24</td><td>26</td><td>23</td><td>24</td><td>18 (5)</td><td>5</td></tr><tr><th>$K = 6$</th><td>9</td><td>12</td><td>10</td><td>12</td><td>4 (10)</td><td>14</td><td>10</td></tr></tbody></table>
<table><caption>**Table 10.** Number of Bethe roots for decomposed BAE with $M = 13$.</caption><thead><tr><th></th><th>{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}</th><th>{13}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>5</td><td>5</td><td>0</td></tr><tr><th>K = 3</th><td>16</td><td>16</td><td>0</td></tr><tr><th>K = 4</th><td>33</td><td>33</td><td>0</td></tr><tr><th>K = 5</th><td>44</td><td>44</td><td>0</td></tr><tr><th>K = 6</th><td>33</td><td>33</td><td>0</td></tr></tbody></table>

<table><caption>**Table 11.** Number of Bethe roots for decomposed BAE with $M = 14$.</caption><thead><tr><th></th><th>{1, 3, 5, 9, 11, 13}</th><th>{2, 4, 6, 8, 10, 12}</th><th>{7}</th><th>{14}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>5</td><td>6</td><td>4 (1)</td><td>6</td><td>1</td></tr><tr><th>K = 3</th><td>20</td><td>19</td><td>20</td><td>18 (1)</td><td>1</td></tr><tr><th>K = 4</th><td>44</td><td>47</td><td>38 (6)</td><td>47</td><td>6</td></tr><tr><th>K = 5</th><td>73</td><td>70</td><td>73</td><td>64 (6)</td><td>6</td></tr><tr><th>K = 6</th><td>69</td><td>74</td><td>54 (15)</td><td>74</td><td>15</td></tr><tr><th>K = 7</th><td>33</td><td>28</td><td>34</td><td>14 (15)</td><td>15</td></tr></tbody></table>

<table><caption>**Table 12.** Number of Bethe roots for decomposed BAE with $M = 15$.</caption><thead><tr><th></th><th>{1, 2, 4, 7, 8, 11, 13, 14}</th><th>{3, 6, 9, 12}</th><th>{5, 10}</th><th>{15}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>6</td><td>6</td><td>6</td><td>6</td><td>0</td></tr><tr><th>K = 3</th><td>23</td><td>24</td><td>22 (1)</td><td>24</td><td>2</td></tr><tr><th>K = 4</th><td>61</td><td>60</td><td>61</td><td>60</td><td>0</td></tr><tr><th>K = 5</th><td>109</td><td>109</td><td>110</td><td>110</td><td>0</td></tr><tr><th>K = 6</th><td>133</td><td>135</td><td>132</td><td>134</td><td>0</td></tr><tr><th>K = 7</th><td>96</td><td>94</td><td>96</td><td>94</td><td>0</td></tr></tbody></table>

<table><caption>**Table 13.** Number of Bethe roots for decomposed BAE with $M = 16$.</caption><thead><tr><th></th><th>{1, 3, 5, 7, 9, 11, 13, 15}</th><th>{2, 6, 10, 14}</th><th>{4, 12}</th><th>{8}</th><th>{16}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>6</td><td>7</td><td>7</td><td>6 (1)</td><td>7</td><td>1</td></tr><tr><th>K = 3</th><td>28</td><td>27</td><td>27</td><td>27</td><td>26 (1)</td><td>1</td></tr><tr><th>K = 4</th><td>77</td><td>80</td><td>81</td><td>74 (7)</td><td>81</td><td>7</td></tr><tr><th>K = 5</th><td>161</td><td>158</td><td>157</td><td>157</td><td>150 (7)</td><td>7</td></tr><tr><th>K = 6</th><td>224</td><td>231</td><td>231</td><td>210 (21)</td><td>231</td><td>21</td></tr><tr><th>K = 7</th><td>218</td><td>211</td><td>211</td><td>211</td><td>190 (21)</td><td>21</td></tr><tr><th>K = 8</th><td>85</td><td>93</td><td>94</td><td>60 (35)</td><td>95</td><td>35</td></tr></tbody></table>
<table><caption>**Table 14.** Number of Bethe roots for decomposed BAE with $M = 17$.</caption><thead><tr><th></th><th>{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}</th><th>{17}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>7</td><td>7</td><td>0</td></tr><tr><th>K = 3</th><td>32</td><td>32</td><td>0</td></tr><tr><th>K = 4</th><td>100</td><td>100</td><td>0</td></tr><tr><th>K = 5</th><td>224</td><td>224</td><td>0</td></tr><tr><th>K = 6</th><td>364</td><td>364</td><td>0</td></tr><tr><th>K = 7</th><td>416</td><td>416</td><td>0</td></tr><tr><th>K = 8</th><td>286</td><td>286</td><td>0</td></tr></tbody></table>

<table><caption>**Table 15.** Number of Bethe roots for decomposed BAE with $M = 18$.</caption><thead><tr><th></th><th>{1, 5, 7, 11, 13, 17}</th><th>{2, 4, 8, 10, 14, 16}</th><th>{3, 15}</th><th>{6, 12}</th><th>{9}</th><th>{18}</th><th>{∞}</th></tr></thead><tbody><tr><th>K = 1</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><th>K = 2</th><td>7</td><td>8</td><td>7</td><td>8</td><td>6 (1)</td><td>8</td><td>1</td></tr><tr><th>K = 3</th><td>37</td><td>36</td><td>38</td><td>37</td><td>38</td><td>36 (1)</td><td>1</td></tr><tr><th>K = 4</th><td>123</td><td>127</td><td>122</td><td>126</td><td>114 (8)</td><td>126</td><td>8</td></tr><tr><th>K = 5</th><td>308</td><td>304</td><td>308</td><td>304</td><td>308</td><td>296 (8)</td><td>8</td></tr><tr><th>K = 6</th><td>550</td><td>559</td><td>552</td><td>562</td><td>524 (28)</td><td>562</td><td>28</td></tr><tr><th>K = 7</th><td>742</td><td>733</td><td>740</td><td>730</td><td>740</td><td>702 (28)</td><td>28</td></tr><tr><th>K = 8</th><td>656</td><td>670</td><td>656</td><td>670</td><td>600 (56)</td><td>670</td><td>56</td></tr><tr><th>K = 9</th><td>276</td><td>262</td><td>279</td><td>265</td><td>280</td><td>210 (56)</td><td>56</td></tr></tbody></table>

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

**References**

[1] R. Baxter, *Exactly solved models in statistical mechanics*, Dover Publications, U.S.A. (2007).

[2] Y. Jiang and Y. Zhang, *Algebraic geometry and Bethe ansatz. Part I. The quotient ring for BAE*, JHEP 03 (2018) 087 [arXiv:1710.04693] [INSPIRE].

[3] J. Salas and A.D. Sokal, *Transfer matrices and partition function zeros for antiferromagnetic Potts models. 1. General theory and square lattice chromatic polynomial*, J. Statist. Phys. 104 (2001) 609 [cond-mat/0004330] [INSPIRE].

[4] J.L. Jacobsen and J. Salas, *Transfer matrices and partition function zeros for antiferromagnetic Potts models. 2. Extended results for square lattice chromatic polynomial*, J. Statist. Phys. 104 (2001) 701 [cond-mat/0011456] [INSPIRE].

[5] J.L. Jacobsen, J. Salas and A.D. Sokal, *Transfer matrices and partition function zeros for antiferromagnetic Potts models. 3. Triangular lattice chromatic polynomial*, J. Statist. Phys. 112 (2003) 921 [cond-mat/0204587] [INSPIRE].

[6] J.L. Jacobsen and J. Salas, *Transfer matrices and partition-function zeros for antiferromagnetic Potts models. IV. Chromatic polynomial with cyclic boundary conditions*, J. Statist. Phys. 122 (2006) 705 [cond-mat/0407444] [INSPIRE].
[7] J.L. Jacobsen and J. Salas, *Phase diagram of the chromatic polynomial on a torus*, Nucl. Phys. B 783 (2007) 238 [cond-mat/0703228] [INSPIRE].

[8] J.L. Jacobsen and J. Salas, *A generalized Beraha conjecture for non-planar graphs*, Nucl. Phys. B 875 (2013) 678 [arXiv:1303.5210] [INSPIRE].

[9] J.L. Jacobsen, J. Salas and C.R. Scullard, *Phase diagram of the triangular-lattice Potts antiferromagnet*, J. Phys. A 50 (2017) 345002 [arXiv:1702.02006].

[10] W. Hao, R.I. Nepomechie and A.J. Sommese, *Completeness of solutions of Bethe's equations*, Phys. Rev. E 88 (2013) 052113 [arXiv:1308.4645] [INSPIRE].

[11] C. Marboe and D. Volin, *Fast analytic solver of rational Bethe equations*, J. Phys. A 50 (2017) 204002 [arXiv:1608.06504] [INSPIRE].

[12] C. Marboe and D. Volin, *The full spectrum of AdS5/CFT4 I: representation theory and one-loop Q-system*, J. Phys. A 51 (2018) 165401 [arXiv:1701.03704] [INSPIRE].

[13] D.A. Cox, J. Little and D. O'Shea, *Ideals, varieties, and algorithms*, 4th edition, Undergraduate Texts in Mathematics, Springer, Germany (2015).

[14] J.C. Faugère, *A new efficient algorithm for computing Gröbner bases (F4)*, J. Pure Appl. Alg. 139 (1999) 61.

[15] J.C. Faugère, *A new efficient algorithm for computing Gröbner bases without reduction to zero (F5)*, in the proceedings of the 2002 International Symposium on Symbolic and Algebraic Computation (ISSAC'02), July 7–10, Lille, France (2002).

[16] W. Decker, G.M. Greuel, G. Pfister and H. Schönenmann, *SINGULAR 4-1-1 — A computer algebra system for polynomial computations*, http://www.singular.uni-kl.de (2018).

[17] C.N. Yang and T.D. Lee, *Statistical theory of equations of state and phase transitions*. I. Theory of condensation, Phys. Rev. 87 (1952) 404.

[18] M. Fisher, *The nature of critical points*, in *Lecture notes in theoretical physics*. vol. 7c, W. Brittin ed., University of Colorado Press, U.S.A. (1965).

[19] S. Beraha, J. Kahane and N.J. Weiss, *Limits of zeroes of recursively defined polynomials*, Proc. Natl. Acad. Sci. 72 (1975) 4209.

[20] D.A. Bini and L. Robol, *Solving secular and polynomial equations: a multiprecision algorithm*, J. Comp. Appl. Math. 272 (2015) 276.

[21] O. Aberth, *Iteration methods for finding all zeros of a polynomial simultaneously*, Math. Comp. 27 (1973) 339.

[22] Y. Jiang and Y. Zhang, *Algebraic Geometry and Bethe Ansatz. Part II. Primary decomposition and algebraic extensions*, to appear.

[23] R.I. Nepomechie and C. Wang, *Algebraic Bethe ansatz for singular solutions*, J. Phys. A 46 (2013) 325002 [arXiv:1304.7978] [INSPIRE].

[24] J.L. Jacobsen and H. Saleur, *Bootstrap approach to geometrical four-point functions in the two-dimensional critical Q-state Potts model: a study of the s-channel spectra*, JHEP 01 (2019) 084 [arXiv:1809.02191] [INSPIRE].

[25] Z. Bajnok et al., *QQ-relations for $U_q(sl_2)$ symmetric spin chains*, to appear.

[26] J. Belletête, A.M. Gainutdinov, J.L. Jacobsen, H. Saleur and R. Vasseur, *On the correspondence between boundary and bulk lattice models and (logarithmic) conformal field theories*, J. Phys. A 50 (2017) 484002 [arXiv:1705.07769] [INSPIRE].
[27] N. Read and H. Saleur, *Enlarged symmetry algebras of spin chains, loop models and S-matrices*, Nucl. Phys. **B** **777** (2007) 263 [cond-mat/0701259] [INSPIRE].

[28] A.M. Gainutdinov, J.L. Jacobsen, N. Read, H. Saleur and R. Vasseur, *Logarithmic conformal field theory: a lattice approach*, J. Phys. **A** **46** (2013) 494012 [arXiv:1303.2082] [INSPIRE].

[29] V. Pasquier and H. Saleur, *Common structures between finite systems and conformal field theories through quantum groups*, Nucl. Phys. **B** **330** (1990) 523 [INSPIRE].

[30] A.M. Gainutdinov, W. Hao, R.I. Nepomechie and A.J. Sommese, *Counting solutions of the Bethe equations of the quantum group invariant open XXZ chain at roots of unity*, J. Phys. **A** **48** (2015) 494003 [arXiv:1505.02104] [INSPIRE].

[31] S. Lang, *Algebra*, Graduate Texts in Mathematics volume 211, Springer, Germany (2002).

[32] F.M. Gomes and D.C. Sorensen, ARPACK++ — An object-oriented version of ARPACK eigenvalue package, https://www.caam.rice.edu/software/ARPACK/arpack++.html (2000).