Algebraic geometry and Bethe ansatz (I) the quotient ring for BAE

In this paper and upcoming ones, we initiate a systematic study of Bethe ansatz equations for integrable models by modern computational algebraic geometry. We show that algebraic geometry provides a natural mathematical language and powerful tools for understanding the structure of solution space of Bethe ansatz equations. In particular, we find novel efficient methods to count the number of solutions of Bethe ansatz equations based on Gr\"obner basis and quotient ring. We also develop analytical approach based on companion matrix to perform the sum of on-shell quantities over all physical solutions without solving Bethe ansatz equations explicitly. To demonstrate the power of our method, we revisit the completeness problem of Bethe ansatz of Heisenberg spin chain, and calculate the sum rules of OPE coefficients in planar $\mathcal{N}=4$ super-Yang-Mills theory.


Introduction
Bethe ansatz is a powerful tool to find exact solutions of integrable models. Ever since the seminal work of Hans Bethe [1], the original method has been developed largely and the term 'Bethe anatz' now refers to a whole family of methods with different adjectives such as coordinate Bethe ansatz, algebraic Bethe ansatz [2,3], analytic Bethe ansatz [4] and off-diagonal Bethe ansatz [5]. A crucial step in the Bethe ansatz methods is to write down a set of algebraic equations called the Bethe ansatz equations (BAE). These equations can be derived from different point of views such as periodicity of the wavefunction, cancelation of 'unwanted terms' and analyticity of the transfer matrix.
The BAE is a set of quantization conditions for the rapidities (or momenta) of excitations 1 of the model, the solutions of which are called Bethe roots. Physical quantities such as momentum and energy of the system are functions of the rapidities. Once the BAE is known, one can solve it to find the Bethe roots and plug into the physical quantities. Therefore, in many cases solving an integrable model basically means writing down a set of BAE for the model. However, in many applications, simply writing down the BAE is not the end of the story. In fact, solving BAE is by no means a trivial task ! Due to the complexity of BAE, it can only be studied analytically in certain limits such as the thermodynamic limit [6] and the Sutherland limit [7] (or semi-classical limit [8,9]). In both cases the size of the system and the number of excitations are large or infinite. For finite system size and number of excitations, typically the BAE can only be solved numerically. While numerical methods are adequate for many applications in physics, they have their limits and shortcomings. Firstly, numerical solutions cannot give exact answers and one needs to find the solutions with high precisions to obtain reliable results. Also, numerical methods might suffer from additional subtleties such as numerical instabilities. Finally and most importantly, the algebraic structure and beauty of BAE can hardly be seen by solving the equations numerically.
From the mathematical point of view, BAE is a set of algebraic equations whose solutions are a collection of points in certain affine space and form a zero dimensional affine variety. It is therefore expected that algebraic geometry may play a useful role in studying the BAE. The first work in this direction was done by Langlands and Saint-Aubin who studied the BAE of six vertex model (or XXZ spin chain) using algebraic geometry [10]. Here we take a slightly different point of view and study BAE from the perspective of modern computational 1 It might also involve some auxiliary variables as in the case of integrable models with higher rank symmetry algebras. algebraic geometry. In particular, we propose that Gröbner basis and quotient ring are the proper language to describe BAE. The aim of our current work is to initialize a more systematic study of the structure of BAE using the powerful tool of algebraic geometry and at the same time developing efficient methods to derive exact results which previously relies on solving BAE numerically.
To demonstrate our points, we study two types of problems with algebro-geometric methods. The first type of problem is a revisit of the completeness problem of Bethe ansatz. This is a longstanding problem for Bethe ansatz which will be discussed in more detail in section 3 and appendix A. Despite the general belief that the Bethe ansatz is complete and many non-trivial progress, this problem does not have a complete and satisfactory solution. In terms of BAE, the completeness problem amounts to counting the number of physical solutions of BAE. Analytical formula for the number of solutions of BAE with various additional constraints in terms of quantum numbers 2 are still unknown 3 even for the simplest Heisenberg XXX spin chain. In order to find the number of solutions for fixed quantum numbers, one needs to solve BAE numerically and find all the solutions explicitly (see for example [11]). From the algebraic geometry point of view, the number of solutions is nothing but the dimension of the quotient ring of BAE which will be defined in section 2. The quotient ring of BAE is a finite dimensional linear space whose dimension can be found without solving any equations ! We propose a method based on Gröbner basis to find the dimension of the quotient ring efficiently.
The second type of problem appears more recently in the context of integrability in AdS/CFT [12,13]. We will give a more detailed introduction to the background of this problem in section 4 and appendix B. The problem can be formulated as the follows. Let us consider a set of BAE with fixed quantum numbers and some additional constraints 4 on rapidities. Typically the number of physical solutions is not unique. Consider a rational function F (u 1 , · · · , u N ) of the rapidities. The problem is to compute the sum of the function F (u 1 , · · · , u N ) evaluated at all physical solutions. The usual way to proceed is first solving BAE numerically and then plugging the solutions in F (u 1 , · · · , u N ) and finally performing the sum. We propose a different approach which avoids solving BAE. The main point is that the function F (u 1 , · · · , u N ) evaluated at the solutions of BAE can be mapped to a finite dimensional matrix called the companion matrix in the quotient ring. The summation over all physical solutions corresponds to taking the trace of this matrix. Importantly, the 2 Such as the length of the spin chain and number of excitations. 3 By this we mean the number of all solutions with pairwise distinct Bethe roots, the number of singular and physical singular solutions. The expected number of physical solutions is of course known from simple representation theory of the symmetry algebra. 4 Such as the condition that the total momentum of the state should be zero.
companion matrix can be constructed in purely algebraic way.
We would also like to mention that similar computational algebro-geometric methods for summing over solutions have been applied to a rather different field, which is the scattering amplitudes [14][15][16]. In the framework of Cachazo-He-Yuan formalism [17][18][19][20][21], the scattering amplitudes can be written as a sum of a given function over all possible solutions of the scattering equations. The scattering equations are also a set of algebraic equations like BAE 5 which can be studied by algebraic geometry. Compared to our case, the scattering equations are much simpler and the structure of the solutions are easier to study. For example, the number of physical solutions can be determined readily and an analytic formula is known.
The rest of this paper is structured as follows. In section 2, we review some basic algebraic geometry that is necessary to understand our methods. In section 3, we study the completeness problem of Bethe ansatz by algebro-geometric methods. We first give a detailed discussion of the physical problem and then provide the method to count the solutions of BAE under additional constraints. In section 4, we propose an analytical method to compute the sum of a given function evaluated at all physical solutions of BAE with fixed quantum numbers. We conclude in section 5 and give a list of open problems and future directions. More backgrounds and technical details are presented in the appendices.

Basics of algebraic geometry
In this section, we briefly review some rudiments of algebraic geometry. We refer to [22][23][24][25] for the mathematical details. See also the lecture notes [26] for the application of computational algebraic geometry for polynomial reductions in scattering amplitudes.

Polynomial ring, ideal and affine variety
Consider a polynomial ring A K = K[z 1 , . . . z n ]. An ideal I of A is a subset of A such that, 2. gf ∈ I, for f ∈ I and g ∈ A.
A polynomial ring is a Noether ring, which means any ideal I of A is finitely generated: for an ideal I, there exist a finite number of polynomials f i ∈ I such that any polynomial F ∈ I can be expressed as We may write I = f 1 , . . . , f k . Given an ideal I, we define quotient ring A/I as the quotient set specified by the equivalence relation: f ∼ g if and only if f − g ∈ I.
We are interested in the common solutions of polynomial equations, or in algebraic geometry language, the algebraic set. The algebraic set Z(S) of a subset S of A is the set in affine spaceK n , HereK is a field extension of the original field K, since frequently we need a field extension to get all the solutions.
It is clear that the algebraic set of polynomials is the same as the algebraic set of the ideal generated by these polynomials, Therefore, we usually only consider the algebraic set of an ideal.

Gröbner basis and quotient ring
An ideal I of A can be generated by different generating sets, or basis. In many cases, a "convenient" basis is needed. For polynomial equation solving and polynomial reduction problems, the convenient basis is the so-called Gröbner basis . A Gröbner basis is an analog of the row echelon form in linear algebra, because it makes the reduction in a polynomial ring possible. (Schematically, the polynomial reduction towards an arbitrary generating set is ill-defined since the result is non-unique, while the polynomial reduction towards a Gröbner basis provides the unique result.) To define a Gröbner basis, we first need to define monomial orders in a polynomial ring. A monomial orders ≺ is a total order for all monomials in A such that, • if u ≺ v then for any monomial w, uw ≺ vw.
• if u is non-constant monomial, then 1 ≺ u.
Some common monomial orders are lex (Lexicographic), deglex (DegreeLexicographic), and degrevlex (DegreeReversedLexicographic). Given a monomial order ≺, for any polynomial f ∈ A there is a unique leading term, LT(f ) which is the highest monomial of f in the order ≺.
A Göbner basis G(I) of an ideal I with respect to a monomial order ≺ is a generating set of I such that for any f ∈ I, (2.4) (Here a|b means a monomial b is divisible by another monomial a). Given a monomial order ≺, the corresponding Göbner basis can be computed by the Buchberger algorithm [27] or more recent F4/F5 [28,29] algorithms. Furthermore, for an ideal I, give a monomial order ≺, the so-called minimal reduced Göbner basis is unique. We give more details on the computation of Gröbner basis in appendix D.
The property (2.4) ensures that the polynomial division of a polynomial F ∈ A towards an ideal I in the order ≺, is well-defined: where g i 's are the elements of the Gröbner basis. r is called the remainder, which contains monomials not divisible by any LT(g i ). Given the monomial order ≺, the remainder r for F is unique.
Therefore, the polynomial division and Gröbner basis method provide the canonical representation of elements in the quotient ring A/I. For two polynomials F 1 and In particular, f ∈ I if and only if its remainder of the polynomial division is zero. This is a very useful application of Gröbner basis since it efficiently determines if a polynomial is inside the ideal or not.

Zero dimensional ideal
A zero dimensional ideal is a special case of ideals such that its algebraic set in an algebraic closed field is a finite set, i.e., |ZK(I)| < ∞. The study of zero dimensional ideals are crucial for our Bethe Ansatz computations.
One of the important properties of a zero dimensional ideal I define over K is that the number of solutions (in an algebraically closed field) equals the linear dimension of the quotient ring Note that the field K need not be algebraically closed, but the field extensionK must be algebraically closed for this formula. Let G(I) be the Gröbner basis of I in any monomial ordering. Since (A K /I) is linearly spanned by monomials which are not divisible by any elements in LT(G(I)), the number of solutions, |ZK(I)| equals the number of monomials which are not divisible by LT(G (I)). This statement provides a valuable method of determining the number of solutions. In practice, we can use the lattice algorithm [24] to list these monomials. If we only need the dimension dim K (A K /I), we can use the command 'syz' in Singular [30].
Let (m 1 , . . . , m k ) be the monomial basis of A K /I determined from the above Gröbner basis G(I). We can reformulate the algebraic structure of (A K /I) as matrix operations. For The k × k matrix c ji is called the companion matrix.
Furthermore, if a polynomial f is in the ideal g + I, we say the fraction f /g is a "polynomial" in the quotient ring A/I by the abuse of terminologies. The reason is that, in this case, f = gq + s, s ∈ I . . For a point ξ ∈ Z(I), if g(ξ) = 0, then f (ξ)/g(ξ) = q(ξ). In this sense, the computation of a fraction over the solution set is converted to the computation of a polynomial over the solutions.
Furthermore, we define M f /g ≡ M q . It is clear that when M g is an invertible matrix, Companion matrix is a powerful tool for computing the sum of values of f evaluated at the algebraic set (solutions) of I over the algebraically closed field extensionK. Let (ξ 1 , . . . , ξ k ) be the elements of |ZK(I)|, Hence this sum over solutions overK can be evaluated directly from the Gröbner basis over the field K. It also proves that this sum must be inside K, even though individual terms may not be.

Application I. Completeness of Bethe ansatz
As a first application of algebro-geometric approach, we revisit the completeness problem of Bethe ansatz in this section. The main calculation is to count the number of solutions of BAE under additional constraints. The usual way of finding the number of solutions is by solving the equations numerically and finding all the solutions explicitly [11,31]. However, if our aim is simply counting the number of solutions, this approach is overkilling. Using algebro-geometric approaches, we can avoid solving BAE and reduce the computation to simple algebraic manipulations.
We start by a detailed discussion on the completeness of Bethe ansatz, using the Heisenberg XXX spin chain as our example. Our goal is to explain why certain kinds of solutions of BAE are 'non-physical' and should be discarded. After that, we present a methods based on Gröbner basis and the quotient ring to count the number of solutions.

Completeness of Bethe ansatz for XXX spin chain
Many integrable models can be solved by Bethe ansatz [1]. In practice this means one has a systematic method to construct the eigenstates of the Hamiltonian and compute the corresponding eigenvalues. The completeness problem of Bethe ansatz is whether all the eigenstates of the Hamiltonian can be constructed by Bethe ansatz. This question turns out to be quite subtle and there is no general answer to it.
In this subsection, we consider the completeness of Bethe ansatz for SU(2) invariant Heisenberg XXX spin chain in the spin- 1 2 representation. There has been arguments for the completeness of Bethe ansatz in the thermodynamic limit where the length of the spin chain is infinite [1,3,32,33]. These arguments are based on the string hypothesis, which needs justification itself. The arguments lead to the correct number of states in the thermodynamic limit but were challenged in the more recent work [11], it is thus still unclear how to justify this kind of arguments in a more rigorous way. When the length of the spin chain is finite, the problem is more difficult and has been investigated in [34][35][36] (see also [37][38][39][40]). In [11] a conjecture for the number of solutions with pairwise distinct roots in terms of the number of singular solutions is proposed. This conjecture has been checked by solving BAE numerically up to L = 14 (see [31] for a generalization to higher spin representations and [41,42] for relations with rigged configurations). We will review this conjecture below. Following this approach, the statement of completeness of Bethe ansatz can be formulated in terms of numbers of solutions of BAE with various additional constraints.
The Heisenberg XXX spin chain is a one-dimensional quantum lattice model with the following Hamiltonian where L is the length of the spin chain and we have imposed periodic boundary condition.
Here σ = (σ 1 , σ 2 , σ 3 ) are the 2×2 Pauli matrices and σ k denotes the spin operator at position k. At each site, the spin can point either up or down, so the Hilbert space has dimension 2 L . The Heisenberg spin chain can be solved by Bethe ansatz [1,3]. In this approach, each eigenstate is labeled by a set of variables {u 1 , · · · , u N } called the rapidities where N is the number of flipped spins. The rapidities satisfy the following BAE The corresponding eigenvalue is given by Naively, one might expect that each solution of BAE corresponds to an eigenstate. However, this is not true and there are solutions of BAE which one should discard. In particular, the following four kinds of solutions need special care 1. Coinciding rapidities. The BAE allows solutions where two of the rapidities coincide, namely u i = u j for some u i , u j ∈ {u 1 , · · · , u N }. For Heisenberg spin chain (3.1), these solutions are not physical and should be discarded. However, we want to mention that whether this kind of solutions are allowed or not in fact depends on the model under consideration [31].

2.
Solutions with N > L/2. The BAE (3.2) can be solved for any N ≤ L. However, when we count the number of physical solutions, we do not consider the cases with magnon number N > L/2. This is because the eigenvectors corresponding to these solutions are not independent from the ones with N ≤ L/2.

3.
Solutions at infinity. The BAE also allows solutions at infinity, namely we can take some u i → ∞. This case corresponds to the descendant states which are necessary for the completeness of Bethe ansatz. However, when we consider the solutions of BAE, we usually count the number of primary states, i.e. no roots at infinity. The number of descendant states of a given primary state can be counted straightforwardly.

Singular solutions.
There are also solutions of BAE at which the eigenvalues diverge (3.3) and the eigenstates are also singular. These solutions are called singular solutions. To determine whether a singular solution is physical or not, one needs to perform a careful regularization. As it turns out, some of the singular solutions are physical and the others are not. The conditions for physical singular solutions are given in [43], which we quote in (A. 19).
For the readers' convenience, we give more detailed discussions on the above points in appendix A.
For the algebro-geometric approach, there's an additional subtlety which is the nontrivial multiplicities of certain solutions. While it is quite normal for algebraic equations to have solutions with multiplicities greater than one, physically we count them as one solution. The number of solutions is counted with multiplicity in algebro-geometric methods and we need to get rid of the multiplicities when counting the number of physical solutions.
By solving BAE for a few cases, we find that the multiple solutions are the ones contain u j = ±i/2, which are the singular solutions. In order to obtain the correct counting, our strategy is to consider separately the singular solutions and the rest ones. To obtain nonsingular solutions, we introduce an auxiliary variable w and add the constraint to the original set of BAE. We see that whenever u j = ±i/2, (3.4) cannot be satisfied.
To obtain the singular solutions, we put u 1 = i/2 and u 2 = −i/2 and solve for the rest variables.
Finally, the completeness of Bethe ansatz can be formulated as a statement of the numbers of solutions of BAE under various constraints. Let us denote the number of pairwise distinct (Pauli principle) finite solutions (primary state) for N ≤ L/2 by N L,N . Among these solutions, we denote the number of singular solutions by N s L,N and the singular physical solutions by N sphy L,N . The number of solutions are counted without multiplicities. The statement of completeness of Bethe ansatz is [11] This is the alluded conjecture in [11]. It has been confirmed by numerics up to L = 14.
The goal of algebro-geometric approach is twofold. The first goal is to provide more efficient and stable methods to find the number of solutions N L,N , N s L,N and N sphys L,N for given L and N and test the conjectures further. The second and more ambitious goal is to find analytical expressions for these numbers in terms of L and N. This requires a careful use of some powerful theorems in algebraic geometry such as the BKK theorem [44][45][46]. While the second goal is not yet achieved in the current work and is still under investigation, we provide an efficient method for the first goal in what follows.

Counting the number of solutions
In this section, we explain how to apply the method of Gröbner basis to compute the numbers N L,N , N s L,N and N sphys L,N for given L and N. The basic idea is that the number of solutions for a given set of polynomial equations is the dimension of the corresponding quotient ring. Instead of solving equations, we construct the quotient rings and compute their dimensions.
For a given L and N, let us define the following polynomials.
where Q u (u) is the Baxter polynomial defined by To have pairwise distinct roots, we define the polynomials This is a classical trick of getting distinct roots in algebraic geometry. For singular and singular physical solutions, we define the following polynomials Using these polynomials, we define the following ideals where the subscribes denote 'Non-Singular', 'Singular' and 'Singular Physical'. The corresponding quotient rings are defined as All We divide the dimensions by factorials to get rid of the permutation redundancy. Any permutation of the set of Bethe roots is considered to be the same solution, yet they correspond to different points in the affine variety. From the definitions of the ideals (3.12), it is straightforward to compute the corresponding Gröbner basis. Then we can construct the standard basis for the quotient rings and the dimensions of the quotient rings follows.

A symmetrization trick
Note that BAE (for non-singular and singular solutions) is totally symmetric in u 1 , . . . u n , i.e., the ideal for BAE is symmetric under the full permutation group of u i 's. We can take advantage of this feature to speed up the Gröbner basis computation. One immediate choice is to apply the symmetric ideal Gröbner algorithm, "symodstd.lib" in Singular. However, this approach is still not fast enough for our propose. Instead, we discovered the following trick: For a totally symmetric ideal I in variables u 1 , . . . u n , we add n auxiliary variables s 1 , . . . s n and n auxiliary equations to make a new idealĨ, (3.15) Therefore we define s k as the k-th elementary polynomials in u 1 , . . . u n . We find that with auxiliary variables and equations, and a block order [u 1 , . . . u n ] ≻ [s n , . . . s 1 ], the Gröbner basis computation is much faster. Furthermore the resulting Gröbner basis forĨ is much shorter comparing with that for I. We believe that the improvement comes from the fact that BAE is much simpler in terms of the symmetric variables s 1 , . . . , s n . The solutions of I' are in one-to-one correspondence to the solution of I, so this method is sufficient.
As a very interesting byproduct, this trick provides a new representation of BAE: The Gröbner basis G(Ĩ), in the block order mentioned above, eliminates the original variables u 1 , . . . u n and gives a set of equations only in s 1 , . . . , s n .   with Table 2. of [43] except for the case L = 12 and N = 5. 6 On a laptop with 16GB RAM and one processor Intel Core i7 without parallelization, we can perform the calculation up to L = 12, N = 6. We use both the software Singular [30] and FGb [47] for this computation.
We comment that this method is very efficient: for example, it only takes about 124 seconds to get the Gröbner basis for the BAE with L = 12 and N = 6, on the laptop mentioned above with the software FGb. Notice that the authors of [43] used clusters to compute these numbers while we are simply using laptops.
Finally, we would like to mention that in parallel with the Gröbner basis method, it is also possible to count the number of Bethe root with the so-called resultant method. The details of this direction are beyond the scope of this paper and we sketch it in the appendix C.

Application II. Sum over solutions of BAE
In this section, we study another kind of problem in integrable systems using algebrogeometric methods. Oftentimes, one encounters the problem of computing the following where the summation runs over all physical solutions 7 of BAE with fixed quantum numbers. For XXX spin chain, the quantum numbers are the length of the spin chain L and the number of particles N. Here F (u 1 , · · · , u N ) is a rational function of the rapidities and might also depend on other parameters. One example of such function is the (square of) OPE coefficient in the planar N = 4 SYM theory which we will discuss below.
The usual way to proceed is first finding all the physical solutions of BAE numerically to very high precisions, plugging into the function F (u 1 , · · · , u N ) and then computing the sum numerically. An interesting observation in [12] is that although each solution of BAE is a complicated irrational numbers and so is the resulting F (u 1 , · · · , u N ), when one sums over all the solutions, the final result gives a simple rational number ! This observation was made by carefully looking at the numerical patterns in the final result.
The numerical approach has certain disadvantages. To start with, finding all solutions of BAE is a highly non-trivial task even for simple models. Secondly, due to numerical instabilities, it is not always easy to estimate to which precision should one be working with in order to find the pattern of rational numbers mentioned above. Finally, it is not clear whether the final result should be a rational number or not.
We propose an alternative method based on algebraic geometry to perform the sum (4.1). Using this approach, there's no need to solve BAE and the computation is reduced to taking traces of numerical matrices whose matrix elements are rational numbers if the coefficients of F (u 1 , · · · , u N ) are rational numbers 8 , which is the case for OPE coefficients. It is then obvious that the final result should be a rational number.
In what follows, we first describe the general method with the help of a simple toy problem. Then we demonstrate how our method works in the context of [12] and how to generalize it to higher loop orders in this case. 7 For a solution to be physical, one usually needs to impose extra selection rules, as was discussed in the completeness problem of BAE. Sometimes, when the quantity under consideration has more symmetry, one can restrict to even smaller subsects of solutions. 8 Our notion of rational numbers also includes complex numbers whose real and imaginary parts are both rational.

Description of the method
In this section, we present more details for the brief discussions on companion matrix in section 2.3 for our current problem. We start with a set of polynomial equations (which we can think as BAE and additional constraints) Let us denote the ideal generated by F 1 , · · · , F n by I F = F 1 , · · · , F n . The corresponding quotient ring is defined by Because the quotient ring Q F is a finite dimensional linear space, it can be spanned by a set basis {e 1 , · · · , e s }. Consider any polynomial P(u 1 , · · · , u N ) in C[u 1 , · · · , u N ]. After imposing the 'on-shell conditions' (4.2), P(u 1 , · · · , u N ) becomes a function in the quotient ring Q F and can be represented in terms of a matrix called the companion matrix.
To be more precise, we have where M P is a numerical matrix of dimension s × s. Here [e j ] denotes the conjugacy class of the basis e j under the identification Our method is based on the following crucial result Two comments are in order. Firstly, the companion matrix M P contains all the information about the on-shell quantity P(u 1 , · · · , u N ). If one diagonalizes the s × s matrix M P , each eigenvalue correspond to P(u 1 , · · · , u N ) with u 1 , · · · , u N at one of the physical solutions of (4.2). Secondly, if the equations (4.2) are symmetric with respect to some of the variables, we should divide the number by proper symmetric factors when performing the sum (4.6) to get rid of permutation redundancy. An alternative way to get ride of the permutation redundancy is to rewrite the polynomial P(u 1 , · · · , u N ) in terms of elementary symmetric polynomials P s (s 1 , · · · , s N ) and perform the calculation in the quotient ring of symmetrized BAE.
The main task is to construct the basis {e j } for Q F and find M P . This can be done using the Gröbner basis. Let us denote the Gröbner basis of I F to be G 1 , · · · , G n , which can be computed from F 1 , · · · , F n and Then the standard basis of quotient ring can be constructed by the method given in section 2.3. The companion matrix can be constructed as follows. First multiply the polynomial P(u 1 , · · · , u N ) with one of the basis e j and then divide the result by the Gröbner basis, where a n are polynomials in u 1 , · · · , u N and P j is the remainder. Since {G 1 , · · · , G n } are Gröbner basis, the remainder P j is well-defined. Now that P j is a polynomial defined in the quotient ring Q F , it can be expanded in terms of the standard basis as This gives the j-th row of the matrix M P . Repeating this process for j = 1, · · · , s, we obtain the matrix M P . In this way, we can construct the companion matrix of any polynomials. To find the companion matrices for rational functions, we can make use the properties of the companion matrices (2.8) and (2.10).

A simple example
To illustrate our approach, we consider a simple example. We first solve the problem by a numerical approach and then by our algebro-geometric approach in order to make a comparison. Let us take F 1 = x 4 y 2 + 3xy + 1, F 2 = y 3 + y 2 − 2 (4.10) and P(x, y) = x 3 3 + y 3 7 + 4xy(x + y) + 2x + 1. It is easy to find that the equations F 1 = F 2 = 0 have 12 solutions. Setting the working precision to 22 digits, we can find the numerical solutions quite easily Finally we take the sum of the 12 values in (4.13) and obtain (4.14) where we see a clear pattern. After rationalization, we obtain simply The result is the following Gröbner basis G 1 = 3xy 2 + 3xy + y + 2x 4 + 1, Now we can construct the standard basis for the quotient ring Q eg = C[x, y]/ G 1 , G 2 . Using the lexicographical ordering for monomials x ≻ y, the standard basis of Q eg is given by Notice that the dimension of Q eg equals the number of solutions of F 1 = F 2 = 0. The next step is to construct the companion matrix M P . Let us first consider e 1 . It is straightforward to calculate that 9 P(x, y)e 1 = a 1 G 1 + a 2 G 2 + P 1 (4.19) where It can be expanded in terms of the basis (4.18) as Working out the other rows (M P ) ij in the same way, we obtain We notice immediately that from the second approach, we directly manipulate the polynomials in a purely algebraic way and there is no need to solve any equations. Therefore we completely avoid all the subtleties of numerical approach. As a bonus, it is clear that the final result should be a rational number since all the manipulations, including the computation of Gröbner basis and companion matrix, involve only simple addition, substraction, multiplication and division of rational numbers and there is no room to create irrational numbers from these operations.

Sum rule of OPE coefficients
In this section, we revisit the calculation of [12] for OPE coefficients in planar N = 4 Super-Yang-Mills theory (N = 4 SYM) using algebro-geometric approach. Let us first give the minimal background of this calculation. It is now well accepted that N = 4 SYM theory is integrable in the planar limit [48]. In practice this means one can use integrability-based methods to compute physically interesting quantities of the theory. For a conformal field theory like N = 4 SYM theory, the most fundamental quantities of interest are the so-called conformal data which consists of the scaling dimensions of all primary operators and the OPE coefficients among these operators.
In order to check the predictions of integrability-based methods, one needs to compare with results from other approaches, such as direct field theoretical calculations based on Feynmann diagrams. The most convenient source of data for OPE coefficients are the fourpoint functions of BPS operators, which are known up to three loops in perturbation theory (see [49] and references therein). By performing operator product expansions of the fourpoint functions, one has access to the information of OPE coefficients. However, it is usually hard to extract a single OPE coefficient from four-point functions. The best one can do is to give predictions for the so-called sum rules defined in (4.25). We give more details of the OPE coefficients and sum rules in appendix B. To summarize, one needs to compute the following quantity is the OPE coefficient of two BPS operators and one non-BPS operator and γ u is the anomalous dimension of the non-BPS operator. They are both functions of the rapidities u ≡ {u 1 , · · · , u S }. The structure constant is given by and (4.28) The quantities in the above expressions such as the momentum p(u), the S-matrix S(u, v) and f (u, v) are known functions of the coupling constant g, where We expand these quantities at weak coupling when g → 0 and consider the result up to 1-loop, namely O(g 2 ) order. We consider the leading order in this subsection and discuss the one-loop result in the next subsection. At the leading order, the various quantities are given by The anomalous dimension γ u only starts to contribute at one-loop order and is given by Let us now consider the sum rule in (4.25). The OPE coefficients depend on four integers L, S, l, N and a set of rapidities {u 1 , · · · , u S }. For fixed L and S, these rapidities satisfy the BAE of SL(2) spin chain In addition, we also need to impose the zero momentum condition The summation in (4.25) runs over all possible solutions of (4.32) and (4.33) for fixed L and S. For generic values of L and S, the solutions of (4.32) and (4.33) are not unique. This is precisely the same type of problem which we discussed in the previous subsection. We can apply our method to perform this sum. Since the coefficients that appear in the sum rule are all rational numbers, it is guaranteed from our approach that the final result will be a rational number as well. We give more details on the implementation of our method in what follows.
We first write down a basis that generate the ideal I S corresponding to (4.32) and (4.33). In order to obtain a polynomial basis, we can write BAE as F 1 = · · · = F S = 0 where F j = (u j + i/2) L Q u (u j + i) + (u j − i/2) L Q u (u j − i), j = 1, · · · , S. Solving these constraints naively, there are solutions with coinciding roots. These solutions are not allowed since they are not physical. To eliminate these solutions, we need to impose extra constraints. These constraints can be imposed in various ways. For example, we can define the following polynomials and impose K ij = 0. The ideal I S is then given by The computations of the Gröbner basis of I S and the basis of the quotient ring Q S = C[u 1 , · · · , u S ]/I S are standard. Once the basis for the quotient ring has been constructed, we can follow the same method described in the previous subsection to construct the companion matrix for the summand F (u 1 , · · · , u S ) = (C ••• u ) 2 e γu . (4.38) As an example, we can consider the case with L = 4, S = 4, l = 2, N = 1. In this case there are 5 allowed solutions and the sum rule (4.25) at the leading order is F = 16/63. We find the dimension of the quotient ring is dim Q S = 120 = 5 × 4!. We use the lattice algorithm [24], implemented in our Mathematica code to determine the 120 monomials in the basis of Q S . As we explained before, the S! permutation redundancy is due to the fact that the BAE and zero momentum condition are completely symmetric with respect to all the rapidities. For this example, our method leads to a matrix M F of 120 × 120 which we will not write down explicitly.
The function F is a rational function and can be written as the ratio of two polynomials F = P/Q. Let us denote their corresponding multiplication matrices as M P and M Q . We then have M F = M P · M −1 Q . Taking the trace of the matrix, we confirm that We checked several other examples and in all the cases, we reproduce the same results as in [12].
To improve the efficiency, we can also use the symmetrization trick in (3.15). Define s i to be the i-th elementary symmetric polynomials in u 1 , . . . , u 4 , i = 1, . .  Here the factor 4! is no longer needed.
Finally, let us comment on the efficiency of our method. For the current calculation, we focus on the SL(2) sector. The BAE of SL(2) spin chain is actually quite easy to solve numerically and the solutions have very nice behaviors, such as all the roots are real and there are no solutions with higher multiplicities.
The analytical counterpart of solving BAE and additional constraints is the construction of the quotient ring. In this specific case, numerical solution is actually faster than constructing the quotient ring. One of the main reasons for this is the permutation redundancy which grows factorially with the number of magnons. This difficulty can be overcome partially by the symmetrization trick described before, but we suspect that further improvements should be possible. Since the elementary symmetric polynomials of rapidities are nothing but the coefficients of Baxter polynomials. It is thus very natural to work with Q-systems instead of BAE. Very recently, Marboe and Volin [50,51] proposed an efficient method to find Bethe roots based on the QQ-relations. Instead of solving BAE, they solve a system of zero remainder conditions (ZRC) where the unknowns are exactly the coefficients of Baxter polynomials. This approach also has other merits such as it automatically selects physical solutions. The method turns out to be quite efficient and works for a large family of spin chains.
Combining our approach and the methods in [50,51], we can in fact construct the quotient ring for ZRC of the corresponding Q-system. In this case, constructing the quotient ring is even faster than solving the ZRC ! For instance, we can construct the quotient ring of ZRC of SU(2) Q-system for L = 14, N = 7 spin chain within several minutes on a laptop. The systematic study of the Gröbner basis and the corresponding quotient ring of rational Q-systems and their applications will be presented in an upcoming publication [52].
After solving the BAE numerically or constructing the quotient ring analytically, one still need to compute the sum over all allowed solutions. This second step is basically trivial for numerical calculations but can be non-trivial for analytical approaches. The main reason is that the quantity we are dealing with, namely the sum rule is a complicated function of rapidities and the complexity grows exponentially with the number of magnon. Unless a simpler form of this quantity is given, this complexity is inherent to any analytical methods and should not be considered as the disadvantage of our algebro-geometric approach. One way to improve the efficiency is to decompose the quantity into smaller and simpler parts, compute the multiplication matrices for the smaller parts and then combine them together. The last step involves only manipulations of numerical matrices and can be done efficiently.

Higher loops
In this section, we discuss how to compute the sum rules at one-loop order. Consider the following sum where the summation runs over the solution of a set of BAE and possibly other additional selection rules that depend on an extra parameter λ. The prototype of this kind of equations we have in mind are the Beisert-Staudacher asymptotic BAE [53] and the zero momentum condition in N = 4 SYM where λ plays the rule of the 't Hooft coupling constant. For simplicity, we assume that the function F (u 1 , · · · , u N ; λ) are rational functions in rapidities {u j }. It depends on λ explicitly as well as implicitly through the rapidities.
We consider the weak coupling limit λ → 0 and develop a perturbative approach to compute the sum (4.43). The leading order is considered to be solved by our approach presented before and we use this knowledge to solve higher loop orders perturbatively. Let us first focus on one-loop order. We assume the summand allows a perturbative expansion in λ At one-loop, the contribution comes from F 0 and F 1 . We also perform a perturbative expansion of the Bethe roots Finally the sum F (λ) can also be expanded in λ The leading order F 0 is considered to be known and are interested in F 1 , which is simply given by where the sum is over leading order solutions. This implies that we only need the quotient ring of the leading order, which is known. The second term in (4.47) is explicit and are rational functions of {u (0) j }, which can be handled straightforwardly as before. The first term involves u where both p(u) and S(u, v) depend on λ. We can expand the above BAE in λ and obtain an approximated BAE valid at one-loop. Then we plug in the ansatz u j = u k . The important point is that the dependence of one-loop BAE in u (1) k is linear. Therefore we can regard u (0) j as constants and solve the linear problem for u (1) k . This can be done straightforwardly and we obtain u k . After we plug U j back into (4.47), the resulting expression is a rational function depending only on u (0) j and the sum can be performed as before. Generalization to higher loops orders is straightforward. The explicit part is easy to deal with. The main complexity of the implicit part comes from expressing u j , we need to solve n linear problems recursively. The procedure is straightforward to implement, but it will lead to increasingly complicated expressions as expected.
We have implemented our algorithm described here and applied it to the sum rule of OPE coefficient at one-loop order. For the example of L = 4, S = 4, l = 1, N = 2, we reproduce exactly [12]  Here, similar to the tree-level case, we can again use the symmetrization trick in (3.15) to simplify the Gröbner basis, companion matrix and the trace computation to get the same answer.
Similar to the leading order, in practice it is more efficient to plug in the solution of BAE numerically than performing an analytic calculation on the quotient ring at higher loop orders. In particular, solving the linear problems numerically are much easier than solving it analytically. However, the analytic method gives exact results without possible loss of accuracy and avoids other subtleties of numerical methods. We should emphasis that our approach here is the most straightforward method, but not necessarily the most efficient one. There is still a huge room to improve the efficiency using the algebro-geometric methods.

Conclusions, discussions and open questions
In this paper we introduced the powerful language of computational algebraic geometry to study Bethe ansatz equations. We developed new analytical methods based on algebraic geometry to tackle two kinds of problems in the framework of Bethe ansatz.
To investigate the completeness problem of Bethe ansatz, we developed efficient methods to count the number of solutions of BAE for fixed quantum numbers with additional constraints. Our method is based on Gröbner basis and the quotient ring and are much faster than solving BAE numerically.
We developed an analytical method to perform the sum of any rational function over all physical solutions of BAE for fixed quantum numbers without actually solving the BAE. We applied our method to calculate the sum rules for OPE coefficients in N = 4 SYM both at tree level and one loop. We obtained exact rational numbers and proved that the results are always rational.
The most prominent advantage of our methods is the conceptual beauty. In the algebrogeometric language, solving BAE is equivalent to constructing the quotient ring of BAE. While the BAE can only be solved numerically, the quotient ring can be constructed analytically and systematically. The quotient ring is a finite dimensional linear space which can be studied much further. Any physical quantity in terms of rapidities is represented as finite dimensional matrix called companion matrix in the quotient ring. Eigenvalues of the companion matrix correspond to the values of this quantity at the solutions of BAE and trace of the companion matrix leads to the sum over all physical solutions at which the physical quantity are evaluated. In addition, constructing quotient rings of BAE is much more efficient than finding explicit solutions of BAE, which makes our method appealing also in practical applications.
There are many open problems that one can pursue in the near future. We list some of them below.
• As we mentioned in section 4.3, Marboe and Volin recently proposed a new method to find physical solutions of BAE based on rational Q-systems [50,51]. The analogy of BAE in this approach is the zero remainder conditions (ZRC) which are much easier to solve numerically. It is very interesting to combine our methods with the rational Q-systems and study the quotient ring of the ZRC. Our preliminary studies have shown that constructing the quotient ring of ZRC is much more efficient than constructing the quotient ring of BAE. It is also naturally much faster than solving the ZRC numerically. If the physical quantities we are interested in are symmetric with respect to the rapidities, we can study them equivalently in the quotient ring of ZRC. This will further boost the efficiency of computing the trace of the corresponding companion matrices.
• Integrable models with symmetries of higher rank Lie algebras are solved by the socalled nested Bethe ansatz [54][55][56], the resulting nested BAEs involve both physical rapidities and auxiliary ones and are hence much more complicated. It is an important future direction to investigate these more challenging cases using algebro-geometric methods. Furthermore, summing over physical quantities at all physical solutions of nested BAE has important applications in the recent work of asymptotic four-point functions [13]. The summand in the generalized sum rules in [13] is symmetric with respect to rapidities and we can apply the quotient ring of ZRC mentioned in the previous point.
• A related question is how to construct companion matrices more efficiently. When the physical quantity under consideration is a complicated function of rapidities, the construction of the companion matrix can be quite tedious although straightforward. This is the main obstacle to the efficiency of our method. One possible way is to decompose the quantities into simpler parts. We can construct the corresponding companion matrices of the simpler parts and then combine them together. The latter step involves only operations on numerical matrices, which should be much easier to handle.
• Concerning the completeness problem of the Heisenberg XXX spin chain, it is desirable to have analytic formula for the various numbers of solutions of BAE with different constraints in terms of L and N. There are several relevant theorems for this type of counting such as the Bézout theorem and the more refined BKK theorem. Using these theorems in an ingenious way and combining some possible local analysis for the special cases, this ambitious goal does not seem to be impossible.
• It is also interesting to investigate the completeness problem of Heisenberg spin chains in higher spin representations where a similar conjecture like the one in (3.5) has been proposed [31]. The generalization to the cases with different boundary conditions [5,[57][58][59] is also an interesting problem.
• One particularly interesting direction is to generalize our current method to the quantum deformed XXZ spin chain. In this case, we have an additional parameter, namely the anisotropy to play with. The completeness problem was investigated in [10] (see also [60] for the generalization to the case of BAE in the asymmetric simple exclusion processes.) It is known from numerics 10 that the structures of solutions of BAE are different in different regimes of anisotropy (see for example [58,63]). It will be fascinating to see this kind of change in the structure of the quotient rings.
• Finally, we only applied the technique of Gröbner basis and resultants to two kinds of problems in the current paper. There are many other powerful tools in algebraic geometry as well as many interesting problmes in integrable models. A wider range of applications and a deeper mutual fertilization could be expected. One particularly interesting example is the computation of exact partition function of integrable lattice models such as six vertex model. A related question is computing grand partition functions of N = 4 SYM theory at one-loop [64].

A More on completeness of BAE
In this appendix, we discuss the four kinds of special solutions in more detail. The discussions below require some basic knowledge about algebraic Bethe ansatz, for which we refer to [3].
Coinciding rapidities If we solve BAE without any constraints, we indeed find solutions of the form {u, u, u 1 , · · · , u N }. They are legitimate solutions of BAE. In the case of coinciding roots, the BAE take a slightly different from which we derive below. Let us recall the RT T relation for XXX 1/2 spin chain Using these relations, one can derive the following result [34,65] where a i (λ, µ) and b i (λ, µ) (i = 1, 2, 3, 4) are some functions of f (µ, λ) and g(µ, λ). As we can see, due to the presence of coinciding rapidities, we have operators B ′ (u) = ∂ u B(u) as well as A ′ (u) = ∂ u A(u) and D ′ (u) = ∂ u D(u). In order the off-shell Bethe state be an eigenstate of the transfer matrix, one computes by moving the diagonal elements A(u) and D(u) 11 to the rightmost and acting on the pseudovacuum using the commutation relations (A.2). This will generate the so-called 'wanted terms' and 'unwanted terms'. By demanding the unwanted terms to vanish, one obtains the usual BAE. There are two modifications in the current case. First of all, the appearance of A ′ (v) and D ′ (v) lead to a ′ (u) and d ′ (u) which might modify the form of the cancelation conditions. In addition, the appearance of B ′ (u) means we need to impose cancelation conditions for the states involving B ′ (u). If we have more coinciding rapidities, from similar analysis we have more additional cancelation conditions.
In fact the cancelation conditions can be obtained most easily by demanding that the eigenvalue of transfer matrix is regular at the Bethe roots. Consider the solution BAE of a spin chain of length L in the spin-s representation with K+N magnons {u, u, · · · , u, u 1 , · · · , u N }. The eigenvalue of the transfer matrix is given by By construction, T (λ) is a polynomial in λ although it seems to have poles at λ = u, u 1 , · · · , u N . By requiring the residues of these 'poles' to vanish, we obtain the BAE. For µ = u j , (j = 1, · · · , N) we have Requiring λ = u is regular leads to the following conditions It was proved in [65] that for the 1D Bose gas where a(u) = e −iuL , d(u) = e +iuL , the BAE B j = R l = 0 do not have solutions for K ≥ 2. For the Heisenberg spin chain, it was found in [34] that there are no solutions with K ≥ 3 and the ones with more than one group of repeated roots such as {u, u, v, v, u 1 , · · · , u N }. However, one can find many solutions of the form {u, u, u 1 , · · · , u N }. Therefore, apart from the general believe that these solutions are not physical, there is no rigorous mathematical proof to this assertion as in the case of 1D Bose gas.
Solutions beyond the equator When looking for physical solutions, we usually restrict ourselves to the regime N ≤ L/2. The BAE itself is well defined also for N > L/2 and explicit solutions can be found. Why do we neglect these solutions ? The answer is that they are already included in the first case. To understand this, let us consider the N < L/2 magnon Bethe state of a spin chain of length L. The Bethe vector can be generated by acting N operators B(u) on the pseudovacuum where the rapidities should satisfy the BAE of N particles. This state has N down spins and L − N up spins. We can generate the eigenstate with the same amount of up spins and down spins by acting L − N operators C(v) on the flipped pseudovacuum Now the rapidities v 1 , · · · , v L−N should satisfy the BAE of L − N particles. As it turns out |Ψ = |Ψ , so (A.9) and (A.9) are merely two ways of constructing the same eigenstate. It is then clear that u = {u 1 , · · · , u N } and v = {v 1 , · · · , v L−N } should be related. This is indeed the case. To see this, one can define the Baxter polynomials It can be shown that the two polynomials satisfy the Wronskian relation, which implies that knowing one of the polynomials gives us the other one. The two polynomials are in fact two solutions of Baxter's T Q-relation which is a second order difference equation. The above analysis shows that we can safely restrict ourselves to one side of the equator N ≤ L/2. The other solutions lead to the same physical states.
Solutions at infinity The Bethe states which correspond to rapidities {u 1 , u 2 , · · · , u N } with none of the elements at infinity is the so-called highest weight state. This means The above relation is non-trivial but can be proved rather straightforwardly. The corresponding spin of this highest weight state is J = L 2 − N. As in quantum mechanics, we can use S − to lower the spins. For a spin-J representation, the dimension is 2J + 1. Therefore, for a highest weight state |u 1 , · · · , u N , the following states (S − ) n |u 1 , · · · , u N , n = 0, · · · , L − 2N (A.13) form a representation space of su(2) algebra. For the completeness of Bethe ansatz, it is thus expected that the number of physical solutions of N-particle BAE should be Then the total number of Bethe states is which is the dimension of the Hilbert space. The solution of BAE allows putting one or more excitations to infinity. Each rapidity at infinity correspond to acting an S − due to the fact Therefore solutions at infinity are allowed and are physical. To show the completeness of Bethe ansatz, we only need to count the solutions that correspond to primary states, the descendants of a primary state is easy to work out. Therefore when we count the solutions, we only count the ones corresponding to primary states.

Singular solutions
The solutions of BAE with two of the rapidities being ±i/2, namely {i/2, −i/2, u 3 , · · · , u N } (A. 17) are called singular solutions. To see that there is a problem at u = ±i/2, it is simplest to look at the eigenvalue in terms of the rapidities It is obvious that the function (u 2 + 1/4) −1 have two poles located at u = ±i/2. Therefore solutions containing u = ±i/2 are special. These solutions are more subtle than the ones we discussed before. The reason is that sometimes these solutions are physical and sometimes not. To see whether a solution is physical or not, one needs to perform a judicious regularization. Such analysis has been worked out in detail in the work of Nepomechie and Wang [43]. The conclusion of their analysis is that the solutions are physical if the remaining rapidities u 3 , · · · , u N satisfy the following equations The first equation is the usual BAE while the second one is an additional selection rule.

B OPE coefficients and sum rules in N = 4 SYM
In this appendix, we give more details about the OPE coefficients and sum rules in the main text. We mainly follow the discussion in [12]. The OPE coefficients can be obtained by computing three-point functions. In our case, we need to compute the three-point function with two BPS operators and one non-BPS operator in the SL(2) sector.
The three operators under consideration are the following. First we have two BPS operators which takes the following form O BPS 1 (x 1 ) = Tr (ZXXZ · · · )(x 1 ) + · · · (B.1) where Z and X are two complex scalar fields andZ,X are the corresponding complex conjugates. The third operator is a non-BPS and takes the following form The wave functions ψ(n 1 , n 2 , · · · , n S ) depend on the Bethe roots, namely the solution of Bethe ansatz equations. The operators O n 1 ,n 2 ,··· ,n S are given by where D is the covariant derivative projected to some light-cone direction D = D µ n µ with n 2 = 0.
For the two BPS operators, the lengths of the operators are defined as the total number of the scalar fields. We denote the lengths of BPS operators to be L 1 and L 2 and the number of scalar fields X (which is equal to the number of scalar fields ofX) to be N. We also define l = L 1 − N, which is the number of scalar field Z for operator O 1 . Let us denote the length (sometimes called twist, which is the number of scalar fields) of the non-BPS operator to be L 3 = L and the total number of covariant derivatives as S. Then we have the following relation The three-point functions of the three operators which we describe above is completely fixed up to a constant called the structure constant, which is the OPE coefficient that appears in the sum rule.
The explicit expression of C ••• u is given in (4.26), (4.27) and (4.28). The non-perturbative expression of the momentum and S-matrix are given by and σ(u j , u k ) is the so-called BES dressing phase [66]. The dressing phase is a rather complicated quantity but it will only start to contribute at three-loops.
We define and expand the sum rule as the follows sol. fixed L and S where y is an auxiliary variable. By computing the sum rule, one has predictions for the numbers P (n,m) S , which can also be obtained from four-point functions in the OPE limit. For more details, we refer to [12]. From the four-point function side, it is clear that P

C Method of resultant
In this appendix, we introduce another method to count the number of solutions of BAE with additional constraints. This method avoids the computation of Gröbner basis and uses another important object of computational algebraic geometry, which is the resultant.
Recall that the multi-variable resultant of the homogeneous polynomials F 0 , · · · , F n ∈ C[x 0 , · · · , x n ] is a uniquely defined polynomial in terms of coefficients of the coefficients of F i with the crucial property that whenever the equations F 0 = · · · = F n = 0 has a non-trivial solution, the so-called Macaulay resultant Res(F 0 , · · · , F n ) = 0 [25]. Our method is based on this fundamental property.
Suppose we have to solve n polynomial equations given by f 1 = · · · f n = 0 where f i ∈ C[u 1 , · · · , u N ]. The polynomials f i (u 1 , · · · , u n ) are not necessarily homogeneous. We then pick one of the variables, say u 1 (We can pick any u k ) and view it as a parameter. Then f i are polynomials depending on variables u 2 , · · · , u n . In order to define the resultant, we introduce another variable u 0 to homogenize the polynomials. Let us denote the homogenized polynomials by F i (U 0 , U 2 , · · · , U n ; u 1 ) 12 , (i = 1, · · · , n) and we have F i (1, u 2 , · · · , u n ; u 1 ) = f i (u 1 , u 2 , · · · , u n ; u 1 ). We can then compute the resultant of the polynomials F i (U 0 , U 2 , · · · , U n ; u 1 ) which is now a polynomial depending on u 1 . We then have q(u 1 ) = Res(F 1 , · · · , F n ).
(C.1) 12 We use capital letters to denote the variables and lower case ones to denote parameters, where The claim is that the number of solutions for the single variable polynomial q(u 1 ) = 0, or equivalently, the highest power of the polynomial q(u 1 ) gives the number of solutions for the original equations f 1 = · · · = f n = 0. 13 Let us illustrate our general procedure by a simple example. We consider the following equations f 1 = f 2 = f 3 = 0 where f i (u 1 , u 2 , u 3 ) is given by We view u 3 as a parameter and introduce another variable u 0 to homogenize the three polynomials, which leads to three homogenized polynomials F i (U 0 , U 1 , U 2 ; u 3 ), (i = 1, 2, 3) Suppose we find a root of q(u 3 ) = 0, denoted byū 3 . Then for u 3 =ū 3 , the equations F 1 = F 2 = F 3 = 0 have non-trivial solutions, which we denote by (U 0 , U 1 , U 2 ). The solution is projective. That is to say for fixedū 3 , if (U 0 , U 1 , U 2 ) is a non-trivial solution, then for any 13 Note that the original Macaulay resultant computation requires the number of equations equals the number of variables. In practice, we may have the situations for which the number of equations is larger then the number of variables. In these cases, the idea of Macaulay can also apply through the evaluation of several Macaulay resultants. For example, suppose that we have n + 1 equations f 1 = .
The main computation in this approach is the multi-variable Macaulay resultant. We find that so far the resultant computation for BAE is complicated and not as efficient as the Gröbner basis method. Since the resultant is given in terms of determinants of large sparse matrices, we expect that in the future, the special Gaussian elimination method optimized for Macaulay matrix can speed up the resultant computation drematically, and make this method applicable for complicated BAE. (For example, the GBLA algorithm described in [67] has a simple method of reducing large Macaulay matrices. However, the specific function for computing Macaulay resultant via GBLA algorithm is not available to the public yet.)

D Computation of Gröbner basis
A Gröbner basis can be computed by various algorithms like Buchberger [27], F4 [28] or F5 [29] algorithms. The classical Buchberger algorithm is the simplest (but may not be the most efficient) algorithm. To provide some intuitions of Gröbner basis computations, in this appendix we first briefly review Buchberger algorithm.
Given two polynomials f and g in a polynomial ring K[x 1 , . . . x n ] with a monomial order ≻, we can define the S-polynomials of f and g as, Here LCM means the least common multiplier, and LT means the leading term of a polynomial in the given monomial order. It is clearly that S(f, g) is a polynomial generated by f and g. So far this algorithm provides the Gröbner basis in the given monomial ordering. In many cases, this kind of Gröbner bases are enough for the practice. However, it is not in the "simplest form", the minimal reduced Gröbner basis. A minimal reduced Gröbner basis is a Gröbner basis such that the leading term from any polynomial in the basis cannot divide any monomial in other polynomials in this basis. The minimal reduced Gröbner basis of an ideal is unique for a given monomial ordering.
To determine the minimal reduced Gröbner basis for this example, we can do the following: note that LT(h 3 )|LT(h 1 ), LT(h 4 )|LT(h 2 ), so h 1 and h 2 are removed from the basis. Furthermore, Buchburger algorithm is simple and intuitive. However, it requires the reduction of many polynomials pairs and can be slow in the practice. In this paper, we mainly used 'slimgb' in the software Singular [30] and the C library 'gb' [67,68] and package 'FGb' [47] for computing Gröbner bases: • 'slimgb' is an improved Buchberger algorithm [69] which smartly picks up the polynomial pairs to reduce the size of intermediate results to speed up the computation.
• 'Fgb' is a modern Gröbner basis package written by Jean-Charles Faugère which applies F4 and F5 algorithm. It can reduce a lot of S-pairs at once and automatically drop the useless S-pairs in the computation.
• 'gb' is a new Gröbner basis C Library written by Christian Eder based on the GBLA algorithm [67] and fast linear algebra techniques.
Sometimes, the Gröbner basis computation over Q is slow. In this case, we can first calculate the Gröbner basis over finite fields Z/p 1 , . . . , Z/p k , where p 1 , . . . , p k are prime numbers. Then we can use Chinese remainder theorem and Farey fractions to lift the finitefield results to rational results. The step can be automated by the package 'modstd lib' in Singular.

E Maple and Mathematica codes
We attach "L12M5.mw", the Maple file for computing the Gröbner basis for nonsingular Bethe roots in SU(2) model with L = 12 and N = 5 with 'FGb' package, and also "SL L4S4.nb", the Mathematica file for computing structure constant in SL(2) model with L = 4 and S = 4, as computation examples.
To run "L12M5.mw", it is necessary to install 'Fgb' package for Maple first. To run "SL L4S4.nb", it is necessary to install Singular and furthermore to download the Mathematica packages for Singular interface and quotient ring basis computations. These Mathematica packages are included in this submission.