Cylinder partition function of the 6-vertex model from algebraic geometry

We compute the exact partition function of the isotropic 6-vertex model on a cylinder geometry with free boundary conditions, for lattices of intermediate size, using Bethe ansatz and algebraic geometry. We perform the computations in both the open and closed channels. We also consider the partial thermodynamic limits, whereby in the open (closed) channel, the open (closed) direction is kept small while the other direction becomes large. We compute the zeros of the partition function in the two partial thermodynamic limits, and compare with the condensation curves.


Introduction
Computing partition functions of integrable vertex models at intermediate lattice size is a hard problem. For small lattice size, the partition function can be computed simply by brute force. For large lattice size, where the thermodynamic limit is a good approximation, various methods are available, including the Wiener-Hopf method [1][2][3], non-linear integral equations [4,5] and a distribution approach [6]. At intermediate lattice size, brute force is no longer an option, and the thermodynamic approximation is inaccurate. In a previous work [7], three of the authors developed an efficient method to compute the exact partition function of the 6-vertex model analytically for intermediate lattice size. They considered the 6-vertex model at the isotropic point on the torus, i.e. with periodic boundary conditions in both directions. The method is based on the rational Q-system [8] and computational algebraic geometry (AG). The algebro-geometric approach to Bethe ansatz was initiated in [9], with the general goal of exploring the structure of the solution space of Bethe ansatz equations (BAE) and developing new methods to obtain analytic results in integrable models. The simplest example for such a purpose is the BAE of the SU(2)-invariant Heisenberg XXX spin chain with periodic boundary conditions. It is an interesting question to generalize these methods to more sophisticated cases such as higher-rank spin chains, quantum deformations and non-trivial boundary conditions. In the current work, we take one step forward in this direction and consider the partition function of the 6-vertex model on the cylinder. Namely, we take one direction of the lattice to be periodic and impose free open boundary conditions in the other direction. This set-up has several new features compared to the torus geometry already considered in [7].
First of all, to consider open boundary conditions for the vertex model, we put the model on a diagonal square lattice where each square is rotated by 45 • , as is shown in figure 1. The partition function on such a lattice can be formulated in terms of a diagonalto-diagonal transfer matrix [10], which does not commute for different values of the spectral parameter. Nevertheless, the R-matrix approach (the so-called Quantum Inverse Scattering Method) can be applied by reformulating this transfer matrix in terms of an inhomogeneous double-row transfer matrix [11] with suitable alternating inhomogeneities [12,13]. It turns JHEP06(2020)169 out that these inhomogeneities depend on the spectral parameter. As a result, the BAE depend on a free parameter ; hence, the Bethe roots are functions of this parameter, instead of pure numbers. In general, this new feature makes it significantly more difficult to solve the BAE. However, in the algebro-geometric approach, there is no extra difficulty, because the computations are purely algebraic and analytic -there is not much qualitative difference between manipulating numbers and algebraic expressions. Therefore, the AG computations can be adapted to cases with free parameters straightforwardly, which further demonstrates the power of our method.
Secondly, in the torus case, the computation of the partition function can be done in two directions which are equivalent. For the cylinder case, however, the computations of the partition function in the two directions or channels are quite different. In the open channel, we need to diagonalize the transfer matrix corresponding to open spin chains. The partition function is given by the sum over traces of powers of the open-channel transfer matrix, similarly to the torus case. In the closed channel, we diagonalize transfer matrices corresponding to closed spin chains, and the open boundaries become non-trivial boundary states. The partition function is thus given by a matrix element, between boundary states, of powers of the closed-channel transfer matrix.
For a given lattice size, the final results should be the same in both channels. Nevertheless, we may consider different limits. In the open (closed) channel, we can take the lattice size in the open (closed) direction to be finite and let the other direction tend to infinity. This partial thermodynamic limit has been studied in the torus case [7]. In the cylinder case, there are two different partial thermodynamic limits ("long narrow straw" and "short wide pancake"), which we study in detail in this paper. In these limits, it follows from the Beraha-Kahane-Weiss theorem [14] that the zeros of the partition function condense on certain curves in the complex plane of the spectral parameter.
The rest of the paper is structured as follows. In section 2, we give the set-up of the vertex model and its reformulation in terms of a diagonal-to-diagonal transfer matrix. In sections 3 and 4, we discuss the computation of the partition function in the two different channels, using Bethe ansatz and algebraic geometry. Section 5 is devoted to some general discussions on the algebro-geometric computations for the BAE/QQ-relation with a free parameter. In section 6 we present the partition functions which can be written in closed forms for any M or N in the two channels. These include M = 1 in the open channel and N = 1, 2, 3 in the closed channel. In section 7, we compute the zeros of the partition function close to the two partial thermodynamic limits (i.e., for very large aspect ratios) and compare them with the condensation curves, which we compute from the transfer matrix spectra. In appendix A, we give a brief introduction to the basic notions of computational algebraic geometry which are used in this paper. In appendix B, we give more details on the algebro-geometric computations. Appendix C to E contain the proofs of some statements in the main text. We collect explicit exact results in appendix F for small M and N , where we can perform the computation in both channels and make consistency checks.
Some of the results we obtained are too large to be presented in the paper. They can also be downloaded from the webpage (576 MB in compressed form): http://staff.ustc.edu.cn/~yzhphy/integrability.html JHEP06(2020)169

Set-up
We consider the 6-vertex model at the isotropic point on a (2M + 1) × 2N medial lattice, for positive integers M and N . We impose periodic boundary conditions in the vertical direction, and free boundary conditions in the horizontal direction; the geometry under consideration is a cylinder, as is shown in figure 1. The partition function on the lattice can be computed in two different channels.
Open channel. In the open channel, we define the diagonal-to-diagonal transfer matrix t D (u) =Ř 23 (u)Ř 45 (u) · · ·Ř 2M,2M +1 (u)Ř 12 (u)Ř 34 (u) · · ·Ř 2M −1,2M (u) , (2.1) shown in figure 2; by convention the direction of propagation (the "imaginary time" direction) is upwards in our figures. The subscripts label the spaces being acted upon, andŘ is related to the standard R-matrix of the isotropic 6-vertex model by where P is the permutation operator. Written explicitly, the R-matrix is given by  Closed channel. In the closed channel, the graph is rotated by 90 • degrees, as shown in figure 3. The rotated R-matrix which will be denoted by R c takes the following form Notice that this R-matrix does not satisfy the Yang-Baxter equation. In the closed channel, the partition function is no longer given by a trace, since periodic boundary conditions are not imposed in the vertical direction. Instead, the open boundary conditions give rise to non-trivial boundary states in the closed channel. The partition function in the closed channel is given by

JHEP06(2020)169
where U is the one-site shift operator U = P 12 P 23 · · · P 2N −1,2N , (2.9) andt D (u) is defined as whereŘ c ij (u) = P ij R c ij (u). The boundary state |Ψ 0 is given by where we have used the notation The result for the partition function of course does not depend on how we perform the computation, so we have To verify the correctness of our various computations (see below), we have explicitly checked this identity for small value of M and N . Our goal is to compute analytic expressions of Z(u, M, N ) explicitly for different intermediate values of M and N . When both M and N are large, the system can be well approximated by the computation in the thermodynamic limit. Here we instead focus on the interesting intermediate case where we keep one of M , N to be finite (namely, the one that determines the dimension of the transfer matrix) and the other to be large. For finite M (≤ 10) and large N (around a few hundred to thousands), we perform the computation in the open channel using (2.5); whereas for finite N and large M , we work in the closed channel using (2.8). We discuss the computation of the partition function in both channels from the perspective of Bethe ansatz and algebraic geometry.

Partition function in the open channel
In this section, we discuss the computation of the partition function in the open channel using Bethe ansatz and algebraic geometry. Using this method, we are able to compute the partition function for finite M ≤ 10 and large N (ranging from a few hundred to thousands).

Reformulation and Bethe ansatz
In order to apply the R-matrix machinery, the first step is to re-express the diagonal-todiagonal transfer matrix t D (u) (2.1) in terms of an integrable open-chain transfer matrix with 2M + 1 sites and with inhomogeneities {θ j } [11]. Let us define

JHEP06(2020)169
where the monodromy matrices are given by For our isotropic problem, the K-matrices are simply K + (u) = K − (u) = I. The eigenvalues Λ(u; {θ j }) of the transfer matrix t(u; {θ j }) (3.1), which can be obtained using algebraic Bethe ansatz [11], are given by The key point (due to Destri and de Vega [12]) is to choose alternating spectralparameter-dependent inhomogeneities as follows One can then show [13] that the diagonal-to-diagonal transfer matrix t D (u) is given by 1 noting that half of the R-matrices become proportional to permutation operators, and the spin chain geometry is transformed into the vertex one, see figure 4. Specifying in (3.3)-(3.4) the inhomogeneities as in (3.5), it follows that the eigenvalues Λ D (u) of t D (u) are given by where the {u k } are solutions of the Bethe equations JHEP06(2020)169 Here k = 1, . . . , K and K = 0, 1, . . . , M . Note that the BAE (3.9) depend on the spectral parameter u, which is an unusual feature. We observe that t D (u) has su(2) symmetry The Bethe states are su(2) highest-weight states, with spin For a given value of K, the corresponding eigenvalue therefore has degeneracy We conclude that the partition function (2.5) is given by where Λ D,K (u) is given by (3.8). Here sol(M, K) stands for physical solutions {u 1 , . . . , u K } of the BAE (3.9) with 2M + 1 sites and K Bethe roots. The number N (M, K) of such solutions has been conjectured to be given by [15] N (M, (3.14) In order to find the explicit expressions for the partition function (3.13), we need to find the eigenvalues Λ D,K (u). They depend on the values of rapidities which are solutions of the BAE (3.9). We encounter two difficulties. Firstly, the solution set of the BAE (3.9) contains some redundancy, since not all solutions are physical; therefore one needs to impose extra selection rules [15]. Secondly, generally Bethe equations are a complicated system of algebraic equations, which cannot be solved analytically. What is worse, our JHEP06(2020)169 BAE (3.9) depend on a free parameter u, which means that the Bethe roots are functions of u, thereby making the BAE even harder than usual to solve.
In order to overcome these two difficulties, we need new tools, namely the rational Q-system and computational algebraic geometry. These methods have been applied successfully in computing the torus partition function of the 6-vertex model [7]. The BAE can be reformulated as a set of QQ-relations, with appropriate boundary conditions [8]. The benefit of working with the Q-system is twofold. Firstly, it is much more efficient to solve the rational Q-system than to directly solve the BAE. Secondly, all the solutions of the Q-system are physical, so there is no need to impose further selection rules [16,17]. The rational Q-system, which was first developed for isotropic (XXX) spin chains with periodic boundary conditions [8], was recently generalized to anisotropic (XXZ) spin chains and to spin chains with certain open boundary conditions [17,18]. We briefly review the Q-system for open boundary conditions in section 3.2.
Turning to the second difficulty, finding all solutions of the BAE (or of the corresponding Q-system) is in general only possible numerically. However, it was realized in [9] that if the goal is to sum over all the solutions of the BAE/Q-system for some rational function f ({u j }) of the Bethe roots, then it can be done without knowing all the solutions explicitly. The idea is based on computational algebraic geometry. The solutions of the BAE/Q-system form a finite-dimensional linear space called the quotient ring. The dimension of the quotient ring is the number of physical solutions of the BAE/Q-system. A basis of the quotient ring can be constructed by standard methods using a Gröbner basis. Once a basis for the quotient ring is known, one can construct the companion matrix for the function f ({u j }), which is a finite-dimensional representation of this function in the quotient ring. Taking the trace of the companion matrix gives the sought-after sum. For a more detailed introduction to these notions and explicit examples in the context of toroidal boundary conditions, we refer to the original papers [7,9] and the textbooks [19,20].
The same strategy can be applied to the open boundary conditions. The new feature that appears in this case is the dependence on a free parameter u. While this creates extra difficulty for numerical computations, it does not cost more effort in the algebro-geometric approach. The reason is that the constructions of the Gröbner basis, the basis for the quotient ring and companion matrices are purely algebraic; and it does not make much qualitative difference whether we have to manipulate numbers or algebraic expressions. 2

BAE and Q-system
In this section, we review the rational Q-system for the SU(2)-invariant XXX spin chain with open boundary conditions [18]. Let us first consider the BAE with generic inhomo- where M a,s is the number of boxes in the Young tableau to the right and top of the vertex (a, s). The boundary conditions are chosen such that Q 2,s = 1, Q 1,s>K = 1 and Here Q 1,0 (v) is the usual Baxter Q-function, whose zeros are the Bethe roots. Comparing to the periodic QQ-relations [8], the main differences are an extra factor v that appears on the left-hand side of (3.16), and the degree of the polynomial of Q a,s which is twice the one for the periodic case. More details can be found in section 4.2 of [18]. For the Bethe equations (3.9), corresponding to the alternating inhomogeneities (3.5), we simply have 3 To solve the Q-system, we impose the condition that all the Q a,s functions are polynomials. This requirement generates a set of algebraic equations called zero remainder conditions (ZRC) for the coefficients c (k) a,s . In principle, one can then solve the ZRC's and find Q a,s , in particular the main Q-function Q 1,0 . The zeros of Q 1,0 are the Bethe roots {u k }, which are functions of the parameter u.
After finding the Q-functions, the next step is to find the eigenvalues Λ D (3.8), which in terms of Q-functions are given simply by . (3.20) Plugging these into (3.13), we finally obtain the partition function.

Algebraic geometry
In this subsection, we give the main steps for the algebro-geometric computation of the partition function: 1. Generate the set of zero remainder conditions (ZRC) from the rational Q-system; 3 Recall that the argument of the double-row transfer matrix t in (3.6) is u 2 , rather than u. Most steps listed above can be done straightforwardly, adapting the corresponding working of [7]. The only step that requires some additional work is step 4. The variables of ZRC are c (k) a,s which are coefficients of the Q-functions. From these variables, it is easy to construct the companion matrix of the Q-function. For fixed M and K, we denote the companion matrix by Q M,K . To find the companion matrix of Λ D , which is essentially the companion matrix of Λ (3.20) up to some multiplicative factors, the most direct way is to use homomorphism property of the companion matrix and write , (3.22) where T M,K (u) is the companion matrix for Λ D (u) with fixed M and K. Unfortunately, this method involves taking the inverse of the matrix Q M,K (u + i 2 ) analytically, which can be slow when the dimension of the matrix is large.
We find that a much more efficient way is to use the following T Q-relation In our case, we need to take L = 2M + 1 and z = u. To solve the T Q relation (3.23), we make the following ansatz for the two polynomials Notice that Q(u) is an even polynomial and only even powers of u appear, which is not the case for T (u). Plugging the ansatz (3.24) into (3.23), we obtain a system of algebraic equations for the coefficients {t 0 , t 1 , · · · , t 2L , s 0 , · · · , s K−1 }. In fact, solving these set of algebraic equations is yet another way to find the Bethe roots. For our purpose, we only solve the equation partially, namely we view {s 0 , · · · , s K−1 } as parameters and solve {t k } in terms of {s j }. This turns out to be much simpler since the equations are linear. We find that t k ({s j }) are polynomials in the variables {s j }. From ZRC and algebro-geometric JHEP06(2020)169 computations, we can find the companion matrix of s j which we denote by s j . Replacing s j by s j and the products by matrix multiplication in t k ({s j }), we find the companion matrix t k = t k ({s j }). Then the companion matrix of the eigenvalues of the transfer matrix is given by More details on the implementation of the algebro-geometric computations are given in appendix B.
Using the AG approach, we have computed the partition functions for M up to 6, with N up to 2048. We also calculated some partition functions with higher M and lower N . The results for 2 ≤ M, N ≤ 6 are given in appendix F.

Partition function in the closed channel
In this section, we compute the partition function in the closed channel. There are both simplifications and complications due to the presence of non-trivial boundary states. Indeed, the presence of boundary states imposes selection rules for the allowed solutions of the BAE. Firstly, it restricts to the states with zero total spin. This implies that the length of the spin chain must be even, which we denote by 2N ; and the only allowed number of Bethe roots is K = N . In contrast, for the periodic (torus) case [7], one must consider all the sectors K = 0, 1, . . . , N . Moreover, the Bethe roots must form Cooper-type pairs (4.33), which leads to significant simplification in the computation of the Gröbner basis and quotient ring.
This simplification comes with a price. Recall that the partition function in the closed channel takes the form of a matrix element given by (2.8). To evaluate this matrix element, we need the overlaps between the boundary states and the Bethe states. These overlaps are a new feature, which is not present in the open channel. They are complicated functions of the rapidities, which makes the computation of the companion matrix more difficult.

Reformulation and Bethe ansatz
To compute the expression (2.8) for the partition function in the closed channel, the first step is to rewritet D (2.10) in terms of integrable closed-chain transfer matrices. To this end, we observe that R c (u) (2.7) is related to R(u) by whereũ is the 'crossing transformed' spectral parameter defined bỹ The corresponding "checked" R-matrices are therefore related by For later convenience, we definẽ

JHEP06(2020)169
in terms of whicht D (2.10) is given bỹ It follows from (4.3) thatṼ where and the V (i) are the same as the correspondingṼ (i) , but with R's instead of R c 's: Let us now introduce integrable inhomogeneous closed-chain transfer matrices of length 2N where the monodromy matrices are defined in (3.2). Using crossing symmetry where 't 1 ' stands for transposition in the first quantum space, one can show that The transfer matrices (4.9) can be diagonalized by algebraic Bethe ansatz. We define the operators A, B, C, D as matrix elements of the monodromy matrix (3.2) as usual but note the shift in the spectral parameter. We consider the reference state or pseudovacuum The Bethe states and their duals are constructed by acting with B-and C-operators on the reference state 4 |u = B(u 1 ) · · · B(u K )|0 , u| = 0|C(u 1 ) · · · C(u K ) . JHEP06(2020)169 with eigenvalues Λ c (u; {θ j }) given by In contradistinction with (3.9) these BAE do not depend on the spectral parameter u, as is usually the case. The eigenvalues Λ c (u; {θ j }) of τ (u; {θ j }) are given by as follows from (4.11).
To make contact witht D , we again choose alternating spectral-parameter-dependent inhomogeneities The V 's (4.8) can then be related to the closed-chain transfer matrices (4.9) by which is similar to (3.6), see figure 5. In view of the relations (4.6) betweenṼ 's and V 's, we conclude thatt D (4.5) is given bỹ

JHEP06(2020)169
The expression (2.8) for the partition function in the closed channel can therefore be recast as where |Φ 0 is the so-called dimer state and we have also used the fact that U † Ω (1) U = Ω (2) . We now insert in (4.22) the completeness relation in terms of Bethe states (which are SU(2) highest-weight states) and their lower-weight descendants where sol c (N, K) stands for physical solutions u of the closed-chain BAE (4.17) with 2N sites and K Bethe roots. Moreover, the normalization factor is given by 25) and the ellipsis denotes the descendant terms. However, these descendant terms do not contribute to the matrix element (4.22), since the dimer state is annihilated by the spin raising and lowering operators Moreover, in view of the fact the overlap u|Φ 0 vanishes unless K = N . The matrix element (4.22) therefore reduces to where the sum runs over all physical solutions of the closed-chain BAE (4.17) in the K = N sector. We note that the expressions involving the eigenvalues are given by JHEP06(2020)169 as follows from (4.16) and (4.18). The normalization factor (4.25) is given by the Gaudin formula [21,22] where and (4.32) Overlaps similar to u|Φ 0 have been studied extensively, see e.g. [23][24][25][26][27][28][29][30], see also [31,32]. The cases of even and odd N must be analyzed separately.

Even N
Let us first consider even values of N . Interestingly, only Bethe states with "paired" Bethe roots of the form have non-zero overlaps [25,26]. Such Bethe states have even parity, see (E.5) below. For such Bethe states, the overlaps are given by (see appendix C) where with and K a (u) is defined in (4.32). We remark that, for these states,

JHEP06(2020)169
Moreover, we show in appendix D the relation The matrix element (4.28) therefore reduces to where we have used (4.29) to pass to the second equality, and the sum is over all physical solutions of the BAE (4.17) with paired Bethe roots (4.33), see (4.48) below.

Odd N
For odd values of N , the only Bethe states with non-zero overlaps have one 0 Bethe root, and all the other Bethe roots form pairs; i.e., the Bethe roots are of the form Such Bethe states have odd parity, see (E.6) below. The overlaps are now given by where H is the block matrix and G ± are now the N −1 2 × N −1 2 matrices given by a (u, v) and K a (u) defined as before, see (4.36), (4.32). Moreover, in (4.42), C is an N −1 2 -component column vector, C t is the corresponding row vector, and D is a scalar, which are given by We remark that, for these states, Moreover, The matrix element (4.28) now reduces to where we have used (4.29) to pass to the second equality, and the sum is over all physical solutions of the BAE (4.17) with paired Bethe roots (4.40), see (4.53) below.

BAE and Q-system
We now summarize the BAE and Q-systems in the closed channel. They are special cases of those for the spin chain with periodic boundary condition with length 2N and magnon number N .
Even N . For the paired Bethe roots (4.33), the closed-chain Bethe equations (4.17) reduce to open-chain-like Bethe equations where Q a,s (v) are even polynomial functions of v. In particular, the main Q-function is given by Moreover, Therefore, to obtain the ZRC for this case, we can simply take the ZRC for the generic periodic case and add the following constraints where k = 1, . . . , (N −1)/2. The corresponding QQ-relations are again given by (4.49), with Q 0,0 (v) given by (4.51). The Q-functions are odd polynomials in this case. In particular, the main Q-function takes the form Therefore, to obtain the ZRC in this case, we take the general ZRC for the generic periodic case and impose the conditions JHEP06(2020)169

Algebraic geometry
The procedure for algebro-geometric computations follows the same steps as in the open channel. As we mentioned before, the computation of the Gröbner basis and quotient ring is simpler. The complication comes from computing the companion matrices. The companion matrix of the transfer matrices Λ c (v; {θ j (u)}) can be constructed similarly from the T Q relation The most complicated part is the ratio of determinants in (4.39) and (4.47). These are complicated functions in terms of rapidities u.
As in the open channel, the natural variables that enter the AG computation are c (k) a,s . Therefore, in order to construct the companion matrices of the ratio of determinants, we need to first convert it to be functions c Odd N . For odd N , the result can be written as Similarly, we can do the symmetry reduction and write the result in terms of the elementary symmetric polynomials . . .

JHEP06(2020)169
They are related to the coefficients c There are two sources of complication worth mentioning. Firstly, computing the determinant explicitly and performing the symmetric reduction is straightforward in principle, but becomes cumbersome very quickly. It would be desirable to have a simpler form for these quantities. Secondly, the companion matrix of the quantity 1/D is the inverse of the companion matrix of D. Computing the inverse of a matrix analytically is also straightforward, but it has a negative impact on the efficiency of the computations when the dimension of the matrix becomes large. For the eigenvalues of the transfer matrix, we saw in (3.23) that the problem of computing inverses can be circumvented by using the T Q-relations. For the expression of the overlaps, it is not clear whether we can find better means to compute the companion matrix of the ratio N/D so as to avoid taking matrix inverses.
Using the algebro-geometric approach in the closed channel, we computed partition functions for N up to 7 and M up to 2048. The results for 2 ≤ M, N ≤ 6 are listed in appendix F.

Algebraic equation with free parameters
In this section, we discuss the Gröbner basis of the ZRC in the closed channel in more detail. This will demonstrate further the power of the algebro-geometric approach for algebraic equations, especially for cases with free parameters.
The system of algebraic equations we consider depends on a parameter u. This means that the coefficients of the equations are no longer pure numbers, but functions of u. As a result, the solutions also depend on the parameter u. As we vary the parameter u, the solutions also change. One important question is if there are any special values u where the solution space changes drastically. To understand this point, let us consider the following simple equation for x whose coefficients depend on the free parameter u At generic values of u, this is a quadratic equation with two solutions. However, when u = ±1, the leading term vanishes and the equation become linear. The number of solutions becomes one. Therefore at these 'singular' points, the structure of the solution space changes drastically. A similar phenomena occurs in the BAE of the Heisenberg spin chain. Consider for a moment the more general XXZ spin chain and take the anisotropy parameter (alias quantum group deformation parameter) q as the free parameter of the BAE. It is wellknown that the solution space is very different between generic q and q being a root of unity. The traditional way to see this is by studying representation theory of the U q (sl(2)) symmetry of the spin chain [33]. A more straightforward way to see this fact is by the algebro-geometric approach. We can compute the Gröbner basis of the corresponding JHEP06(2020)169 BAE/Q-system and analyze the coefficients as functions of q. We shall discuss this problem in more detail in a future publication.
Related to the current work, we consider the ZRC in the closed channel for the XXX spin chain. Here the free parameter is the inhomogeneity u. We want to know whether there are special singular points of u where the structure of solution space changes drastically. Recall that from elementary algebraic geometry, the number of solutions equals the linear dimension of the quotient ring. Furthermore, the quotient ring dimension is completely determined by the leading terms of the Gröbner basis. Therefore, to this end, we compute the Gröbner basis explicitly. For N = 3, the ideal can be written as g 1 , g 2 , g 3 where the elements of the Gröbner basis g i are given by Here we have chosen the ordering We see from (5.2) that the leading terms are independent of u. This implies that the dimension of the quotient ring C[s 0 , s 1 , s 2 ]/ g 1 , g 2 , g 3 , or equivalently the number of solutions, is independent of the value of u. Of course, the explicit solutions of the BAE will depend on the value of u, but there will always be 3 solutions to the ZRC for N = 3 at any value of u. Similarly, we can write down a slightly more non-trivial example for N = 4. The ideal is given by g 1 , · · · , g 6 where the Gröbner basis elements g i are given by The leading terms are in boldface letters. We see again these terms are independent of u.
For all the values of N which we compute, this is true. It would be nice to prove this for general N .

JHEP06(2020)169
Therefore from the algebro-geometric computation, we conclude that for any value of u, there exist solutions with definite parity (parity even/odd for even/odd N ). For fixed N , the number of solutions is the same for any value of u, which has been given in (4.56).
We end this section by the following comment. The conclusion that there always exist solutions with definite parity for any u is far from obvious from the ZRC or original BAE. It is also not easy to see this from numerical computations. On the contrary, it is a straightforward observation from the algebro-geometric computation. This shows again that algebro-geometric approach is a powerful tool to analyze the solution space of BAE.

Analytical results in closed form
In this section, we discuss the analytical results which can be written in closed forms for arbitrary N and M in the open and closed channel respectively.

Open channel
We first discuss the open channel. The partition function takes the same form as the torus case, which is written as the trace of the N -th power of the transfer matrix. If the eigenvalues of the transfer matrix can be found analytically, we can write down the partition function for any N . Here by analytical we mean more precisely expressible in terms of radicals. In the algebro-geometric approach, we first compute the companion matrix of the eigenvalue of the transfer matrix. The dimension of the companion matrix equals the number of physical solutions of the open channel BAE/Q-system. The eigenvalues of the companion matrix give the eigenvalues of the transfer matrices evaluated at each solution. From Galois theory, if the dimension of the companion matrix is less than 5, the eigenvalues can be expressed in terms of radicals. Therefore for the values of M where all the companion matrices have dimension less than 5, we can obtain the analytical expression for any N . This requirement is only met by M = 1. Already for M = 2, we need to consider the sectors K = 0, 1, 2, and for K = 1 and 2 the dimensions of the companion matrices are 4 and 6 respectively. For larger M , the dimensions of the companion matrices are even larger. We give the closed form expression for M = 1 and any N in what follows.
The M = 1 case. We need to consider K = 0, 1. For K = 0, the eigenvalue of the transfer matrix is given by The two eigenvalues of the transfer matrix in this sector are given by The closed-form expression of the partition function, taking into account the su(2) multiplicities (3.12), is then

JHEP06(2020)169
Let us make one comment on the comparison with the torus case. The closed-form results have been found up to 5 M = 6 in the torus case [7]. There we also used the fact that certain companion matrices can be further decomposed into smaller blocks, which implies the existence of non-trivial primary decompositions over the field Q. Physically, this primary decomposition is related to decomposing the solutions of BAE/Q-system according to the total momentum. In the cylinder case, however, the total momentum is automatically zero for all allowed solutions, due to the presence of the boundary. Therefore, further decomposition according to the total momentum is not possible in the current case.

Closed channel
The situation is more interesting in the closed channel. The expression for the partition function is qualitatively different from the torus case, since we have a new ingredient: the non-trivial overlap between Bethe states and the boundary state. To find the analytical expressions, we first compute the companion matrices, both for the transfer matrix and the overlaps. For the values of N where the dimensions of the companion matrices are less than 5, we can express the final result in terms of radicals for any M . This is satisfied by N = 1, 2, 3. The dimensions of the companion matrices are 1, 2, 3 respectively. We present the analytical results for these cases in what follows.
The N = 1 case. This case is somewhat trivial, but we give it here for completeness. There is only one allowed solution to the BAE, which is {0}. The eigenvalue of the transfer matrix is given by The contribution from the overlaps only comes from det H jk which is given by The partition function is thus given by (4.47).
The N = 2 case. This is the simplest non-trivial case where N is even. The Bethe roots take the form {u 1 , −u 1 }. There are two such solutions, which can be found straightforwardly by directly solving Bethe equations. Let us denote the companion matrices of Λ c ũ 2 , {ũ 2 } by T 2 (u) and companion matrix of the following factor Note that in the torus case, M denoted the length of the spin chain [7], while here, in the cylinder case, the length of the spin chain is given by 2M + 1.
Collecting all the results, the explicit closed-form expression for N = 2 and any M is given by One may check that this agrees in particular with (2.6) for M = 2. We see here that the eigenvalues take rather complicated forms in terms of radicals whose arguments are polynomials of u. Nevertheless, the final result is a polynomial, as it should be.
The N = 3 case. This is the simplest non-trivial case where N is odd. The results are bulky, therefore it is more convenient to write them in terms of smaller building blocks.
To this end, we recall the solution for cubic polynomial equations. Let us consider the following generic cubic equation where a = 0. We define Then the three solutions of the cubic equation (6.12) are given by

JHEP06(2020)169
For N = 3, the solutions of the Bethe equations take the form {u 1 , −u 1 , 0}. There are three physical solutions. Let us denote the companion matrix of Λ c ũ 2 , {ũ 2 } by T 3 (u) and the companion matrix of the following factor by F 3 (u). The partition function is given by where the trace is over the 3-dimensional quotient ring. One can check explicitly that T 3 (u) and F 3 (u) commute with each other and can thus be diagonalized simultaneously.
Let us denote their eigenvalues by λ T,i (u) and λ F,i (u), with i = 1, 2, 3. The characteristic equations of T 3 (u) and F 3 (u) take cubic forms where the coefficients are rational functions of u. The characteristic equations can be solved by radicals using (6.15). The relevant quantities are given as follows. For the eigenvalues of T 3 (u), we have and The eigenvalues {λ T,1 , λ T,2 , λ T,3 } are given by where C T is defined as in (6.14). For the eigenvalues of F 3 , we have and where

JHEP06(2020)169
The eigenvalues are given by 6 Finally, combining all the results, the closed form expression is given by

Zeros of partition functions
The study of partition function zeros is a well-known tool to access the phase diagram of models in statistical physics. The seminal works by Lee and Yang [34] and by Fisher [35] studied the zeros of the Ising model partition function, respectively with a complex magnetic field (at the critical temperature) and at a complex temperature (in zero magnetic field). But more generally, any statistical model depending on one (or more) parameters can be studied in the complex plane of the corresponding variable(s). In particular, the chromatic polynomial with Q ∈ C colors has been used as a test bed to develop a range of numerical, analytical and algebraic tools for computing partition function zeros and analyzing their behavior as the (partial) thermodynamic limit is approached [36][37][38][39][40][41][42]. Further information about the physical relevance of studying partition function zeros can be found in [43] and the extensive list of references in [36].
In the case at hand, we are interested in zeros of the partition function Z(u, M, N ) of the six-vertex model, in the complex plane of the spectral parameter, u ∈ C. As explained in section 2, the algebro-geometric approach permits us to efficiently compute Z(u, M, N ) close to the partial thermodynamic limits N M (open channel) or M N (closed channel), and more precisely for aspect ratios ρ := N/M of the order ∼ 10 3 and ∼ 10 −3 , respectively.

Condensation curves
An important result for analyzing these cases is the Beraha-Kahane-Weiss (BKW) theorem [14]. When applied to partition functions of the form (3.13) for the open channel, respectively (4.39) or (4.47) for the closed channel, it states that the partition function zeros in the partial thermodynamic limits (ρ → ∞ or ρ → 0, respectively) will condense on a set of curves in the complex u-plane that we shall refer to as condensation curves. In particular, the condensation set cannot comprise isolated points, or areas. By standard theorems of complex analysis, each closed region delimited by these curves constitutes a thermodynamic phase (in the partial thermodynamic limit). 6 Notice that the powers of ξ in (6.25) are slightly different from (6.21). The reason for this convention is to make sure that λ T i and λ F i correspond to the same eigenvector. Working directly with characteristic equations, it is not immediately clear which eigenvalues correspond to the same eigenvector. We establish the correspondence by making numerical checks. We choose u to be some purely imaginary numbers such that the arguments in the radicals are real and positive.

JHEP06(2020)169
To be more precise, let Λ i (u) denote the eigenvalues of the relevant transfer matrix (for the open or closed channel, respectively) that effectively contributes to Z(u, M, N ). For a given u, we order these eigenvalues by norm, so that |Λ 1 (u)| ≥ |Λ 2 (u)| ≥ · · · , and we say that an eigenvalue Λ i (u) is dominant at u if there does not exist any other eigenvalue having a strictly greater norm. Under a mild non-degeneracy assumption (which is satisfied for the expressions of interest here), the BKW theorem [14] states that the condensation curves are given by the loci where there are (at least) two dominant eigenvalues, |Λ 1 (u)| = |Λ 2 (u)|. It is intuitively clear that this defines curves, since the relative phase φ(u) ∈ R defined by Λ 2 (u) = e iφ(u) Λ 1 (u) is allowed to vary along the curve. Moreover, a closer analysis [36] shows that the condensation curves may have bifurcation points (usually called T-points) or higher-order crossings when more than two eigenvalues are dominant. They may also have end-points under certain conditions; see [36] for more details.
A numerical technique for tracing out the condensation curves has been outlined in our previous paper on the toroidal geometry [7]. It builds on an efficient method for the numerically exact diagonalization of the relevant transfer matrix, and on a direct-search method that allows us to trace out the condensation curves. We refer the reader to [7] for more details, and focus instead on a technical point that is important (especially in the closed channel) for correctly computing the condensation curves for the cylindrical boundary conditions studied in this paper.
One might of course choose to obtain the eigenvalues by solving the BAE, either analytically or numerically. However, the Bethe ansatz does not provide a general principle to order the eigenvalues by norm. It is of course well known that in many, if not most, Betheansatz solvable models, for "physical" values of the parameters the dominant eigenvalue and its low-lying excitations are characterized by particularly nice and symmetric arrangements of the Bethe roots, and hence one can easily single out those eigenvalues. However, we here wish to examine our model for all complex values of the parameter u, and it is quite possible -and in fact true, as we shall see -that there will be a complicated pattern of crossings (in norm) of eigenvalues throughout the complex u-plain. To apply the BAE one would therefore have to make sure to obtain all the physical eigenvalues and compare their norm for each value of u. By contrast, the numerical scheme (Arnoldi's method) that we use for the direct numerical diagonalization of the transfer matrix is particularly well suited for computing only the first few eigenvalues (in norm), so we shall rely on it here. We shall later compare the computed condensation curves with the zeros of partition functions obtained using Bethe ansatz and algebraic geometry.
The reader will have noticed that above we have twice referred to the diagonalization of a "relevant" transfer matrix. By this we mean a transfer matrix whose spectrum contains only the eigenvalues that provide non-zero contributions to Z(u, M, N ), after taking account of the boundary conditions via the trace (2.5) in the open channel, or the sandwich between boundary states (2.8) in the closed channel. These contributing eigenvalues correspond to the physical solutions in (3.13) for the open channel, or in (4.39) and (4.47) for the closed channel. A "relevant" transfer matrix is thus not only a linear operator that can build up the partition function Z(u, M, N ), but it must also have the correct dimension, namely K N (M, K) given by (3.14) in the open channel, or N (N ) given by (4.56) in the JHEP06(2020)169 closed channel. Ensuring this is an issue of representation theory. We begin by discussing it in the open channel, which is easier.

Open channel
The defining ingredient of the transfer matrix is theŘ-matrix. Using (2.2)-(2.3), it readš with a(u) = u + i, b(u) = u and c(u) = i. The most immediate transfer matrix approach is to let the diagonal-to-diagonal transfer matrix t D (u) given by (2.1) act in the 6-vertex model representation, that is, on the space {| ↑ , | ↓ } ⊗2M +1 of dimension 2 2M +1 .
If we constrain to a fixed magnon number K, the dimension reduces to 2M +1 K . This is larger than N (M, K) given by (3.14), because we have not restricted to su(2) highestweight states. Therefore, each eigenvalue would appear with a multiplicity given by (3.12). Since each eigenvalue actually does contribute to Z(u, M, N ), dealing with this naive representation provides a feasible route to computing the condensation curves (and this was actually the approach used in [7]). However, the appearance of multiplicities is cumbersome and impedes the efficiency of the computations.

Temperley-Lieb algebra
To overcome this problem, notice that in the more general XXZ model with quantum-group deformation parameter q, the integrableŘ-matrix may be taken aš for certain coefficients α, β depending on u and q. Here I denotes the identity operator and E i is a generator of the Temperley-Lieb (TL) algebra. The defining relations of this algebra, acting on L = 2M + 1 sites, are 3) where i, j = 1, 2, . . . , L − 1 and the parameter δ := q + q −1 . A representation of E i , written in the same 6-vertex model representation as (7.1), reads

JHEP06(2020)169
By taking tensor products, one may check that this satisfies the relations (7.3). We can match with (7.2) by taking The trick is now that there exists another representation of the TL algebra having exactly the required dimension N (M, K). The basis states of this representation are link patterns on L sites with d := L − 2K defects. A link pattern consists of a pairwise matching of L − d = 2K points (usually depicted as K arcs) and d defect points, subject to the constraint of planarity: two arcs cannot cross, and an arc is not allowed to straddle a defect point. We show here two possible link patterns for L = 5 and d = 1 (hence M = 2 and K = 2): and The TL generator E i acts on sites i and i + 1 by first contracting them, then adding a new arc between i and i + 1. This can be visualized by placing the graphical representation E i = on top of the link pattern. If a loop is formed in the contraction, it is removed and replaced by the weight δ. If a contraction involves an arc and a defect point, the defect point moves to the other extremity of the arc. If a contraction involves two distinct arcs, the opposite ends of those arcs become paired by an arc. For instance, the action of E 1 on the two link patterns in (7.6) produces and δ× (7.7) Recall from (3.11) that the spin s associated with the K-magnon sector in the chain of L = 2M + 1 sites reads s = L 2 − K. The generators E i can decrease s by contracting a pair of defects and replacing them by an arc. It is however possible to define a representation of the TL algebra in which s is fixed, by defining the action of E i to be zero whenever there is a pair of defects at sites i and i + 1. In the literature on the TL algebra, these representations in terms of link patterns with a conserved number of defects are known as standard modules and denoted W s . Meanwhile, in TL representation theory, the partition function in the open channel is no longer written in terms of a trace as in (2.5). Instead it is written as a so-called Markov trace (1 + 2s) q tr Ws t D (u) N . (7.9)

JHEP06(2020)169
We have here defined the q-deformed numbers where U p (x) denotes the p-th order Chebyshev polynomial of the second kind. This result (7.9) can be proved by using the quantum group symmetry U q (su (2)) enjoyed by the spin chain in the open channel [44], or alternatively by purely combinatorial means [45]. The factors (1 + 2s) q appearing in (7.9) account for the multiplicities in the problem. In the limit q → −1 corresponding to the XXX case of interest, the q-deformed numbers become (n) q = n for n odd, and (n) q = −n for n even. The latter minus sign can be eliminated at the price of an overall sign change of the partition function (7.9), since only even n = 1 + 2s occur in the problem. This corresponds to q → −q, so that the quantum group symmetry U q (su(2)) becomes just ordinary su(2) in the limit. The multiplicities then become (1 + 2s) q = 1 + 2s = 2M + 2 − 2K, in agreement with (3.11).
In conclusion, we see that not only do the link pattern representations of the TL algebra lead to the correct dimensions N (M, K), but they also account for the correct su(2) multiplicities 2M + 2 − 2K of eigenvalues in the XXX spin chain.

Results
We have computed the condensation curves by applying the numerical methods of [7] to the transfer matrix t D (u) given by (2.1). The latter is taken to act on the representation given by the union of link patterns on L = 2M + 1 sites with K ∈ {0, 1, . . . , M } arcs and d = L − 2K defects.
The results for the condensation curves with M = 2, 3, 4, 5 are shown in figure 6. The curves are confined to the half-space Im u ≤ 0, and they are invariant under changing the sign of Re u. Therefore it is enough to consider them in the fourth quadrant: Re u ≥ 0, Im u ≤ 0. The condensation curves display several noteworthy features:   both with finite-size corrections proportional to 1/M . We conclude that the leftmost point of the connected component converges to the same value as the end-point, namely u = −i. This kind of "pinching" is characteristic of a phase transition [34,35]; note however that the limit u → −i of the XXX model is singular and does not present a critical point in the usual sense.
We now compare the condensation curves with the partition function zeros. The partition functions Z(u, M, N ) were first computed from the algebro-geometric approach, for M = 2, 3, 4, 5 and N = 1024, which corresponds to a very large aspect ratio ρ ∼ 10 3 . The zeros of Z(u, M, N ) were then computed by the program MPSolve [46] [47], which is a multiprecision implementation of the Ehrlich-Aberth method [48,49], an iterative approach to finding all zeros of a polynomial simultaneously.
The resulting zeros are shown in figure 7, as red points superposed on the condensation curves of figure 6. The agreement is in general very good, although some portions of the condensation curves are very sparsely populated with zeros; in those cases the zeros are still at a discernible distance from the curves, in spite of the large aspect ratio. We have verified that the agreement improves upon increasing ρ. 7 JHEP06(2020)169

Closed channel
In the closed channel theŘ-matrix can be inferred from (2.7) and (2.2). It readš still with a(u) = u + i, b(u) = u and c(u) = i. However, as we shall soon see, it is convenient to apply a diagonal gauge transformation D = diag(1, −1) in the left in-space and the right out-space ofŘ c ; that is,Ř c 12 → D 1Ř c 12 D 1 . This has the effect of changing the sign of c(u) while leaving the partition function unchanged: the gauge matrices square to the identity at the intersections betweenŘ-matrices when taking powers of the transfer matrixt D (u) given by (2.10). To complete the transformation, the first and last row of gauge transformations have to be absorbed into a redefinition of the boundary states Ψ 0 | and |Ψ 0 appearing in (2.8).

Temperley-Lieb algebra
As in the open channel, we can rewrite theŘ-matrix in terms of TL generators (7.4): To match (7.13), with c(u) = −i after the gauge transformation, we must now set In the closed channel, the TL algebra is defined on L = 2N sites. The goal is now to find a representation having the same dimension N (N ), see (4.56), as the number of physical solutions appearing in the closed-channel expressions of the partition function, (4.39) and (4.47). This issue is more complicated than in the open channel.
As a first step, we let the TL generators act on the basis of link patterns, as before. Since the boundary states restrict to zero total spin, the only allowed number of Bethe roots is K = N (see section 4). This implies that the link patterns are free of defects (d = 0). The transfer matrixt D (u) is then given by (2.10) with (7.14), where the TL generators E i act on the link patterns as described in section 7.2.1. To reproduce the partition function (2.8) we also need to interpret the boundary state (2.11) within the TL representation. The natural object is the quantum-group singlet of two neighboring sites which is represented in terms of link patterns as a short arc joining the neighboring sites. We should however remember at this stage the gauge transformation that allowed us to switch the sign of c(u). To compensate this, we need to insert a minus sign for a down-spin in the second tensorand, to obtain With q = −1, this is proportional to |ψ 0 of (2.11).

JHEP06(2020)169
On the other hand it is easy to check from (7.4) that the TL generator E i is nothing but the (unnormalized) projector onto the quantum-group singlet. Therefore, just as E i = , the initial boundary state |Ψ 0 can be represented graphically by the defect-free link pattern in which sites 2j − 1 and 2j are connected by an arc, for each j = 1, 2, . . . , N . Similarly, the final boundary state Ψ 0 | is interpreted as the TL contraction of the corresponding pairs of sites. With these identifications, we have explicitly verified for small N and M that the TL formalism produces the correct partition functions, such as (2.6).
With the spin-zero constraint imposed, the TL dimension is thus equal to the number of defect-free link patterns on L = 2N sites. This is easily shown to be given by the Catalan numbers Although this is smaller than the dimension of the 6-vertex-model representation constrained to the S z = 0 sector, viz. 2N N , it is not as small as (4.56), so further work is needed.
The transfer matrix and the boundary states are also symmetric under cyclic shifts (in units of two lattice spacings) of the L = 2N sites. This symmetry can be used to further reduce the dimension of the transfer matrix. Indeed, after acting witht D (u) we project each link pattern obtained onto a suitably chosen image under the cyclic group Z N . In this way each orbit under Z N is mapped onto a unique representative link pattern. The dimension of the corresponding rotation invariant transfer matrix then reduces to [50]  But one can go a bit further, since the transfer matrix and boundary states are also invariant under reflections. This gives rise to a symmetry under the dihedral group D N . The dimension of the rotation-and-reflection invariant transfer matrix then becomes [50]  The process of imposing more and more symmetries and reducing the dimension of the relevant transfer matrices might be realized at the ZRC level by imposing more and JHEP06(2020)169 more constraints on the Q-functions. Consider the ZRC of a closed spin chain with length L = 2N and N magnons. The corresponding Bethe states are in the S z = 0 sector. In order to restrict to the parity symmetric solutions, we need to impose the condition Q(u) = Q(−u) for any u. This leads to further constraints to the ZRC and reduces the number of allowed solutions down to N (N ). Since Q(u) is a polynomial of order N , the constraint Q(u) = Q(−u) can also be imposed by Q(x k ) = Q(−x k ) at N different values. Now the main observation is that at certain values of x k , the constraints have a clear physical meaning. For example, taking x 1 = i/2, the constraint Q(i/2) = Q(−i/2) is equivalent to which restricts to the solutions with zero total momentum. It is therefore an interesting question to see whether the dihedral symmetry can be realized in this way. If so, at which further value(s) of x k would we need to impose Q(x k ) = Q(−x k ) ?
We do not presently know if and how one can identify a TL representation whose dimension equals N (N ) = N N/2 given by (4.56). It certainly appears remarkable at this stage that as already noticed in [50]. It is also worth pointing out that N (N ) can be interpreted as the number of defect-free link patterns on 2N sites which are symmetric around the mid-point. We leave the further investigation of this question for future work.

Results
We have computed the condensation curves in the closed channel, using the TL link-pattern representations identified above, namely using: 1) spin-zero (i.e., defect free) link patterns, 2) spin-zero link patterns with cyclic symmetry Z N , and 3) spin-zero link patterns with dihedral symmetry. Figure 8 shows the results using only the spin-zero constraint. In the regions enclosed by curves of blue color there is a unique dominant eigenvalue, whereas in the regions enclosed by green (resp. yellow) color the dominant eigenvalue has multiplicity two (resp. three). Obviously the sought-after representation of dimension N (N ) is expected to be multiplicity-free, so the corresponding condensation curve should be free of green and yellow branches. Nevertheless, the curves in figure 8 are expected to correctly produce the condensation curves of partition function zeros for any system described by the transfer matrixt D (u) and with boundary states that impose only the spin-zero symmetry, while breaking any other symmetry (e.g., by imposing spatially inhomogeneous weights).
Next we show in figure 9 the results using both the spin-zero and the cyclic symmetry Z N . When compared to figure 8 it can be seen that many branches of the curves are unchanged. However all of the yellow and some of the green curves have now disappeared, reflecting the fact that the eigenvalues which were formerly dominant inside the regions enclosed by green and yellow colors have now been eliminated from the spectrum, since they do not correspond to Z N symmetric eigenstates. The curves in figure 9 should give the correct condensation curves for systems having the spin-zero and cyclic symmetries.  But those with N = 4, 5 should even provide the correct results for the full "Cooper-pair" symmetry (4.33), since the TL dimension is then equal to the number of physical solutions.

JHEP06(2020)169
Finally we depict in figure 10 results using the spin-zero and dihedral symmetry D N . For N = 6 the condensation curve is identical to the one found with cyclic symmetry, meaning that none of the four eliminated eigenvalues (when going from dim Z N (6) = 28 to dim D N (6) = 24) was dominant anywhere in the complex u-plane. It should provide the correct result for the full "Cooper-pair" symmetry (4.33) if the elimination of four more eigenvalues (going from dim D N (6) = 24 to N (6) = 20) turned out to be equally innocuous. The N = 7 curve with dihedral symmetry has only branches corresponding to multiplicity-free eigenvalues, so it may also apply to the full paired symmetry, although a greater amount of eigenvalues are redundant in this case.
All the curves contain an end-point u e close to the origin for which we have found the following results: L=14 T-points End-points Figure 9. Condensation curves for partition function zeros on a (2M + 1) × 2N cylinder, in the limit M → ∞ (closed channel), for a system exhibiting spin-zero and cyclic symmetry Z N . The panels show, in reading direction, the cases N = 4, 5, 6, 7.  To finish this section, we now compare the condensation curves with the actual partition function zeros. This is done for N = 4, 5, 6, 7 in figure 11. For the partition function zeros, we have M = 2048, except for N = 7 where we have only M = 1024; this ensures an aspect ratio ρ < 10 −2 in all cases. The agreement with the condensation curves appears excellent, with the possible exception of the bubble-shaped region with −5 < Im u < −2 in the N = 7 case.  Figure 11. Comparison between the partition function zeros on a (2M + 1) × 2N cylinder, with M = 2048, and the corresponding condensation curves in the M → ∞ limit (closed channel). The panels show, in reading direction, the cases N = 4, 5, 6, 7. To overcome these challenges, we have further developed the application of algebro-geometric methods to the Bethe ansatz equations in various directions. We have incorporated recent developments in integrability such as rational Qsystems for open spin chains [18] and the exact formulae for overlaps between integrable boundary states and Bethe states [23][24][25][26][27]. We have also developed powerful algorithms to perform the algebraic geometry computations, such as the construction of Gröbner bases and companion matrices, in the presence of a free parameter.

Conclusions and discussions
Equipped with these new developments, we obtained the following exact results for the cylinder partition function Z(u, M, N ). We studied the partial thermodynamic limit of the partition function in both channels using the exact results. In particular, we computed the zeros of the partition functions in these limits and found that they condense on certain curves. The condensation curves in the partial thermodynamic limit can be found by a numerical approach based on the BKW theorem. This numerical approach has been applied in the torus case and was further developed in the current context by taking into account the new features, especially in the closed channel. Comparing the distribution of the zeros obtained from the exact partition function and the condensation curve obtained from the numerical approach, we found nice agreement and were able to shed light on several interesting features. The condensation curves in both the open and closed channels were found to involve very intricate features with multiple bifurcation points and enclosed regions. We believe that the further study of these curves might be of independent interest. There are many other questions which deserve further investigation. One of the most interesting directions is to compute the partition function for the q-deformed case. For generic values of q (namely, when q is not a root of unity), Q-systems for both closed and open chains have been formulated in a recent work [18]. This should provide a good starting point for developing an algebro-geometric approach, since QQrelations are more efficient than Bethe equations and give only physical solutions. It would presumably be easier to first study the torus case where the relevant Bethe equations are those of the periodic XXZ spin chain. After that, one could move to the more complicated cylinder case.
We have focused here on the cylinder geometry with free boundary conditions. The case of fixed boundary conditions (with two arbitrary boundary parameters) may now also be in reach, using the new Q-system [51].
From the perspective of algebro-geometric computations, it will be desirable to sharpen the computational power of our method. For instance, we will try to apply the modern implements of Faugère's F4 algorithm [53], which is in general more efficient than Buchberger's algorithm.

A Basic notions of computational algebraic geometry
In this appendix, we give a brief introduction to some basic notions of computational algebraic geometry which are used in the main text.

A.1 Polynomial ring and ideal
Polynomial ring. Let us start with the notion of polynomial ring which is denoted by A K [z 1 , . . . , z n ] or A K for short. It is the set of all polynomials in n variables z 1 , z 2 , . . . , z n whose coefficients are in the field K. In our case, the field is often taken to be the set of complex numbers C or rational numbers Q.
Ideal. An ideal I of A K is a subset of A K such that 1. f 1 + f 2 ∈ I, if f 1 ∈ I and f 2 ∈ I, 2. gf ∈ I, for f ∈ I and g ∈ A K . Importantly, any ideal I of the polynomial ring A K is finitely generated. This means, for any ideal I, there exists a finite number of polynomials f i ∈ I such that any polynomial F ∈ I can be written as We can write I = f 1 , f 2 , . . . , f k . Here the polynomials {f k } are called a basis of the ideal.

A.2 Gröbner basis
As mentioned before, an ideal is generated by a set of basis {f 1 , . . . , f k }. The choice of the basis is not unique. Namely, the same ideal can be generated by several different choices of bases Notice that in general k does not have to be the same as s. For many cases, a convenient basis is needed. For solving polynomial equations and the polynomial reduction problem, it is most convenient to work with a Gröbner basis.

JHEP06(2020)169
We can introduce the Gröbner basis by considering polynomial reduction. A polynomial reduction of a given polynomial F over a set of polynomials {f 1 , . . . , f k } is given by where r is a polynomial that cannot be reduced further by any of the f i . The polynomial r is called the remainder of the polynomial reduction. One important fact is that the polynomial reduction is not unique for a generic basis {f 1 , . . . , f k }. As a simple example, let us take f 1 = y 2 − 1 and f 2 = xy − 1, and consider F (x, y) = x 2 y + xy 2 + y 2 . The polynomial reduction of F (x, y) over f 1 and f 2 can be performed in two different ways F (x, y) = f 1 + (x + y)f 2 + (x + y + 1).
As we can see, the remainders are r 1 = 2x + 1 and r 2 = x + y + 1 respectively. Therefore for a generic basis, the remainder of the polynomial is not well-defined. A basis of an ideal {g 1 , . . . , g s } is said to be a Gröbner basis if the polynomial reduction is well-defined in the sense that the remainder is unique. Now we move to the more formal definition of the Gröbner basis.
Monomial ordering. To define a Gröbner basis, we first need to define monomial orders in the polynomial ring. A monomial order ≺ is specified by the following two rules • If u ≺ v then for any monomial w, we have uw ≺ vw.
• If u is non-constant monomial, then 1 ≺ u.
Leading term. Once a monomial order ≺ is specified, for any polynomial f , we can define the leading term uniquely. The leading term, which is denoted by LT(f ), is defined as the highest monomial of f with respect to the monomial order ≺.
Gröbner basis. A Gröbner basis G(I) of an ideal I with respect to the monomial order ≺ is a basis of the ideal {g 1 , . . . , g s } such that for any f ∈ I, there exists a g i ∈ G(I) such that LT(f ) is divisible by LT(g i ). A Gröbner basis for a given ideal with a monomial order can by computed by standard algorithms such as Buchberger algorithm [52] or the F4/F5 algorithms [53]. For an ideal I, given a monomial order ≺, the so-called minimal reduced Gröbner basis is unique.

A.3 Quotient ring and companion matrix
Quotient ring. Given an ideal I, we can define the quotient ring A K /I by the equivalence relation: f ∼ g if and only if f − g ∈ I.
Polynomial reduction and Gröbner basis provide a canonical representation of the elements in the quotient ring A K /I. Two polynomials F 1 and F 2 belong to the same element

JHEP06(2020)169
If M g is an invertible matrix, we can actually define the companion matrix of the rational function f /g by The companion matrix is a powerful tool for computing the sum over solutions of the polynomial system (A.5). As we mentioned before, the dimension of Q I equals the number of solutions of (A.5). Let us denote the N solutions to be ( ξ 1 , . . . , ξ N ). Then we have the following important result

B More details on AG computation
In this section, we summarize the algorithm of our algebra-geometry based partition function computation for this paper.
• We first compute the Gröbner basis and quotient ring linear basis of the TQ relation equations (3.23) and QQ relation equations (4.51). Note that the TQ and QQ relations contain the free parameter u. It is possible to compute the corresponding Gröbner basis analytically in u via sophisticated computational algebraic-geometry algorithms, like "slimgb" in the software Singular [54].
However, we find that it is more efficient to set u to some integral value, and compute the Groebner basis. The computation is done with the standard Gröbner basis command "std" in Singular. In this approach, we maximize the power of parallelization since the Groebner basis running time for different values of u is quite uniform.
• Then we compute the power of companion matrices (T M,K ) N . Although from the Gröbner basis it is straightforward to evaluate (T M,K ) N , T M,K is usually a dense matrix and the matrix product is a heavy computation. Instead, we postpone the matrix computations to the end, and evaluate the polynomial power F N first. Here F is the corresponding polynomial of T M,K . After each polynomial multiplication step, we divide the polynomial by the ideal's Gröbner basis to save RAM usage by trimming high-degree terms. To speed up the computation, we apply the binary strategy, i.e., F N = F N/2 F N/2 . After F N is calculated, a standard polynomial division computation provides the companion matrix power (T M,K ) N . So the partition function for a particular u value is obtained.
• In previous steps, u is set as an integral number. To get the analytic partition function in u, we have to repeat the computation, and then interpolate in u. We know that for both the closed and open channel partition functions, the maximum degree in u is 2M N . Hence, we compute the partition function with 2M N +1 integer values. These results are then interpolated to the analytic partition function in u.
The interpolation is carried out with the Newton polynomial method.

JHEP06(2020)169
The whole computation is powered by our codes in Singular. The parallelization is implemented in the Gröbner basis and companion matrix power steps, for different integer values of u's. The interpolation step is not parallelized, although it is also straightforward to do so in the future.
We remark that through the computations, the coefficient field is chosen to be the rational number field Q. We observe that the resulting analytic partition function contains large-integer coefficients, so finite-field techniques may not speed up the computation.

C The overlap (4.34)
The overlap (4.34) can be deduced from results in the paper by Pozsgay and Rákos [27] (based on [23,24]), to which we refer here by PR. Their R-matrix is given by PR (2.6), which has the same form as ours (2.3), except with Moreover, they work with the "quantum monodromy matrix" given by PR (2.33)

D The relation (4.38)
We show here that the relation (4.38) follows from two simpler lemmata.

E Parity of states with paired Bethe roots
Following [55,56], the parity operator Π in the closed channel (length 2N ) is defined by where X n is any operator at site n ∈ {1, 2, . . . , 2N }, and is given by The dimer state |Φ 0 (4.23) is an eigenstate of parity with eigenvalue (−1) N which is consistent with the fact that the overlaps Φ 0 |u are nonzero only for Bethe states with paired Bethe roots (4.33), (4.40).