On the Q operator and the spectrum of the XXZ model at root of unity

The spin-1/2 Heisenberg XXZ chain is a paradigmatic quantum integrable model. Although it can be solved exactly via Bethe ansatz techniques, there are still open issues regarding the spectrum at root of unity values of the anisotropy. We construct Baxter's Q operator at arbitrary anisotropy from a two-parameter transfer matrix associated to a complex-spin auxiliary space. A decomposition of this transfer matrix provides a simple proof of the transfer matrix fusion and Wronskian relations. At root of unity a truncation allows us to construct the Q operator explicitly in terms of finite-dimensional matrices. From its decomposition we derive truncated fusion and Wronskian relations as well as an interpolation-type formula that has been conjectured previously. We elucidate the Fabricius-McCoy (FM) strings and exponential degeneracies in the spectrum of the six-vertex transfer matrix at root of unity. Using a semicyclic auxiliary representation we give a conjecture for creation and annihilation operators of FM strings for all roots of unity. We connect our findings with the 'string-charge duality' in the thermodynamic limit, leading to a conjecture for the imaginary part of the FM string centres with potential applications to out-of-equilibrium physics.


Conclusion 54
A Quantum sl 2

Introduction
To study the dynamics of a quantum many-body system, it is of vital importance to know the full spectrum, i.e. all eigenstates of the Hamiltonian. For instance, with the knowledge of the spectrum it is possible to calculate the density matrix, which is the central object to study the entanglement properties of many-body systems [1]. Furthermore, one can use the spectrum to study the quantum quench problem, a paradigmatic example of out-of-equilibrium dynamics following the logic of the Quench Action [2,3]. The spectrum of quantum Hamiltonian is also closely related to the partition function of the corresponding classical statistical physics model, which can be used to detect phase transitions [4].
However, it is very difficult to obtain the full spectrum, especially for strongly interacting systems. Some models can be mapped into a free fermion (parafermion) model through a non-local mapping [5][6][7][8], resulting in relatively simple spectra. However, generally such a construction does not exist for an interacting many-body system.
Fortunately, there are strongly interacting models amenable to exact methods. Using quantum integrability [9,10] it is possible to obtain the spectrum. For these models, the transfer matrices are the generating functions of the conserved (quasi-)local charges that contribute to the dynamics in the thermodynamic limit. In this setting the generalized Gibbs ensemble (GGE) [11] has proven a crucial tool to study quench problems in integrable models [12], with a close relation to the transfer matrix approach. For quantum integrable models Bethe-ansatz techniques provide an exact characterisation of the spectrum. The eigenvalues of the transfer matrix are obtained via a set of rapidity parameters whose physical values, called Bethe roots, obey the Bethe equations [13]. The latter follow in a straightforward way from a difference equation known as Baxter's TQ relation [9,[14][15][16]. The problem of determining whether the solutions of the Bethe equations provide a complete set of eigenstates has drawn attention in both mathematics and physics [17][18][19][20][21][22][23][24]. 1 The results in this article can help to understand the completeness at root of unity by describing and illustrating different structures present in the spectrum.
It is in principle possible to obtain the full spectrum of the spin-1 2 XXZ chain at root of unity by solving Baxter's TQ relation. The main problem is that there are (a large number of) degenerate eigenstates of the transfer matrix, which cannot be easilly interpreted as solution to the (functional) TQ relation [25][26][27][28][29]. One of the difficulties specific to the root of unity case, is that given a solution of the TQ relation which corresponds to an eigenstate, it is possible to multiply the Q function by certain polynomial factors left undetermined by the TQ relation, which provide other solutions called descendants. When these new solutions turn out to correspond to true eigenstates, they can be interpreted as the presence of bound states called 'exact strings' (Fabricius-McCoy strings) that do not scatter with other excitations [16,20,[30][31][32][33][34][35][36]. In Ref. [32], Fabricius and McCoy proposed to organise the descendants into irreducible representations of the loop algebra of sl 2 , cf. [37][38][39]. They can be obtained by acting with creation operators on highest weight states. These representations have dimension 2 n and are characterised by a degree-n polynomial called the Drinfeld polynomial which, unlike the regular eigenvalues of the Q operator, shows when the descendants are present.
In this article, we propose an operatorial approach to this problem. Namely, building on Refs. [36,40] we define a two-parameter transfer matrix that is used to construct the Q operator, i.e. the solution to the matrix TQ relation. Once the Q operator is obtained, all the physical solutions to Bethe equations can be extracted as the zeroes of the eigenvalues, including the solutions associated to exact strings. Baxter found a Q operator for the sixvertex model by directly solving the matrix TQ relation [9]. Our Q operator is closer to that proposed by Korff [33]. It can be easily implemented on a laptop for system size N ≤ 16 when the denominator 2 of the root of unity is less than 10. We illustrate with several explicit examples the appearance of Fabricius-McCoy strings as well as Bethe roots at infinity. In addition, we prove in an elementary way some conjectured results in Refs. [32,34,41]. Closely related to the Q operator, we obtain a second solution of the TQ relation corresponding to the polynomial P of Pronko and Stroganov [25]. The eigenvalues of the P and Q operators both contain factors associated to Fabricius-McCoy strings. The product of these factors coincides with the Drinfeld polynomial of [32,39], see also Section 8.5. In this setting, the degeneracy of the loop-algebra multiplets is recovered from the different ways to decompose the Drinfeld polynomial into two factors belonging to Q and P, respectively, cf. Section 7.1. We present conjectures for the creation and annihilation operators for eigenstates associated to the degeneracies in Section 9.
Furthermore, we find that the existence of Fabricius-McCoy strings is closely related to the quasilocal Z charges [42][43][44] that lie at the origin of the fractal structure of the spin Drude weight in the spin- 1 2 XXZ model at root of unity [45]. The 'string-charge duality' [46], 1 The Bethe ansatz has been proven to be complete for the XXZ spin chain with generic inhomogeneities and twist [18]. At a practical level these results are sufficient as one can calculate physical quantities for the inhomogeneous model and take the homogeneous limit at the end. In this work we focus on the homogeneous XXZ spin chain. We do not address completeness at a rigorous level here. a functional relation that links the root density of bound states to the generating functions of conserved charges based on the thermodynamic Bethe ansatz, can be illustrated by considering the spectrum with Fabricius-McCoy strings at root of unity. In addition, we give a conjecture about the string centres of Fabricius-McCoy strings based on the string hypothesis.
Outline. We start in Section 2 with some examples to illustrate the problems that one encounters at root of unity. The main body of the paper can be divided into two parts. Sections 3-6 deal with the formalism. After a short introduction to the framework of the quantum inverse scattering method in Section 3, we define the two-parameter transfer matrices and derive their factorisation properties in Section 4. These results allow us to construct the Q operator explicitly as well as transfer matrix fusion and Wronskian relations for arbitrary anisotropy parameter ∆ in Section 5. In Section 6 we prove the truncated transfer matrix fusion and Wronskian relations at root of unity, leading to an interpolation formula that has been previously conjectured in Refs. [32,34,41].
Sections 7-10 deal with the applications to understanding the spectrum of the XXZ spin chain (and six-vertex model) at root of unity. Equipped with the exact form of the Q operator, we demonstrate the descendant tower structure with numerous examples in Sections 7-8. We present a conjecture for the Fabricius-McCoy string creation and annihilation operators in Section 9. Discussions on the thermodynamic limit of solutions with Fabricius-McCoy strings follow in Section 10.
We conclude in Section 11. Some technical details are provided in the Appendices.

Motivation
Before we start with the technical part, let us describe features that motivate our study of the spectrum of the XXZ spin chain at root of unity.
Note that the representation (2.5) breaks the symmetry of (2.1) under parity, i.e. spatial reflection. We will also encounter the parity-conjugate (opposite) representation of U q (sl 2 ), which we denote byS z = S z andS ± = S ± | q →q −1 . We call such representations on the spin-chain Hilbert space (C 2 ) ⊗N 'global'. See Appendix A for more about U q (sl 2 ) and its representations.
Here the momentum is defined such that e ip is the eigenvalue of the (twisted, right) translation operator, see Appendix B.2. The Bethe equations (2.7) can be rewritten as where the scattering phase for two Bethe roots is S(p m , p n ) = − e i(pm+pn) − 2 ∆ e ipn + 1 e i(pm+pn) − 2 ∆ e ipm + 1 = sinh(u m − u n − η) sinh(u m − u n + η) .
(2.11) By slight abuse of notation we will denote the final expression by S(u m , u n ). Let us stress that, more precisely, by 'Bethe roots' we will mean the roots of the Q function, i.e. the eigenvalue of Baxter's Q-operator, see Section 5.5. Each set {u m } of Bethe roots in the strip −π/2 < Im u ≤ π/2 3 specifies an eigenstate. The distribution of Bethe roots in the complex plane is crucial to understand the thermodynamic properties of integrable systems, leading to thermodynamic Bethe ansatz (TBA). As we will see in Section 10 the study of the spectrum at finite system size N may already shed light on the origin of string-charge duality [46] which is important to study dynamics in the thermodynamic limit.
The main question that we set out to understand is what is the structure of the spectrum of the XXZ spin chain at root of unity? Despite of its exact characterisation this is not yet understood systematically. In the following two subsections we discuss phenomena that motivate our study. Although we focus on the spin-chain perspective, note that these phenomena are equally relevant for the root-of-unity case of the six-vertex model, whose transfer matrix gives rise to the XXZ Hamiltonian (2.1), see Section 3.

Bethe roots at infinity
When ∆ = ±1 the spectrum exhibits degeneracies reflecting the model's isotropy (sl 2 -invariance). Although there are additional degeneracies in the spectrum of H due to parity invariance, the latter are lifted when taking into account the momentum -or any other parity-odd charge generated by the transfer matrix, or indeed the transfer matrix itself. Let us focus on the degeneracies in the joint spectrum and return to isotropy. The lowering operator j σ − j of sl 2 can be thought of as adding a magnon with vanishing quasimomentum; indeed, it can be shown that a Bethe-ansatz vector has highest weight iff p m = 0 mod 2π [10]. The isotropic limit of (2.9) (rescale u m η u m and let η → 0) shows that p m = 0 corresponds to u m = ±∞. That is, for the XXX spin chain Bethe roots at infinity signal descendant states. Now switch on the anisotropy parameter ∆ = ±1. Then j σ ± j are no longer symmetries, so one might expect there to be no more Bethe roots at infinity. However, numerical solutions of the Bethe equations (e.g. via the recipe from Appendix C.1) show that infinite Bethe roots are in fact present for the XXZ spin chain.
To see when Bethe roots at infinity occur we turn to the Bethe equations (2.7). Write n ±∞ for the number of Bethe roots at ±∞, so that of course n +∞ + n −∞ ≤ M . Note that p m → ∓iη for u m → ±∞, while S(u m , u n ) → e ∓2η as long as u n is finite or goes to ∓∞. Let us assume that the roots at +∞ do not scatter (S = 1) amongst each other, and likewise for those at −∞. Then the Bethe equation for the infinite root u m = ±∞ becomes [20,30,55] exp ± N − 2 (M − n ±∞ ) η = exp(−iφ). (2.12) Let us examine the possible values of n ±∞ ∈ N in each regime of the XXZ spin chain. For η ∈ R, i.e. in the gapped regime (|∆| > 1), and twist φ ∈ R the only solutions are This implies that physical solutions with Bethe roots at infinity only exist when M > N/2 ('beyond the equator') and N is even, and that the infinite Bethe roots come in pairs. Next consider the gapless regime (|∆| < 1). When η ∈ i (R\π Q), so not at root of unity, there are more possibilities to allow for Bethe roots at infinity, as long as the twist φ is tuned to an (even or odd, depending on the parity of N ) integer multiple of i η: (2.14) In particular, the number of roots at +∞ and −∞ do not coincide if φ = 0. Finally, at root of unity η = iπ 1 / 2 the condition becomes In this case there is an additional condition that we will discuss momentarily, see (2.17). The meaning of n ±∞ = 0 for the XXZ model can be understood from the algebraic Bethe ansatz (see Section 3): Bethe roots at infinity correspond to applications of the lowering operators of the global U q (sl 2 ) algebra (see Appendix A.1). Namely, when N is even and φ vanishes, each eigenstate beyond the equator (say at M > N/2) is obtained from a Bethe eigenvector at M = N − M < N/2 by M − N/2 = N/2 − M applications of the parityinvariant product S −S− . Up to an overall factor the result is the spin-reverse of the Bethe vector we started with. That is, we find that If N is odd the XXZ spin chain is still invariant under a global spin flip, reversing ↑ ↔ ↓ everywhere, but there are no infinite roots and we have not been able to find a simple relation between the Bethe roots {v m } M m =1 and {u m } M m=1 on the two sides of the equator. At root of unity we have to be a little more careful since S − andS − are nilpotent. Here an additional requirement is needed to ensure that (2.16) is nonzero: In particular, in the periodic case (φ = 0) there is at most one nonzero solution to (2.15) in the range (2.17), leaving only three possible scenarios: The machinery that we shall develop in Sections 4-6 will help to understand the structure present in the spectrum due to such infinite Bethe roots. The conclusion (2.18) will be useful in Section 8.1 when we discuss applications to the spectrum of the XXZ model.

Fabricius-McCoy strings
At root of unity -including, but certainly not limited to, the free-fermion point ∆ = 0, η = iπ/2 -the spectrum of the XXZ spin chain has many degeneracies [37]. Fabricius and McCoy realised [30,31] This describes a bound state, which as a whole does not scatter with other Bethe roots: A FM string is not only 'transparent' in scattering, but has vanishing energy; indeed, it does not carry any local charge generated by transfer matrix T 1/2 (u), except that it may carry momentum π as we explain Section 8. It actually carries a specific type of quasilocal charges called the Z charges [42,43,45]; the physical implications of this will be discussed in Section 10.2. When any FM string is present amongst the Bethe roots for a given eigenstate, it is not possible to determine the location of the string centre α fm by solving the (functional) TQ relation using method in Appendix C.1. In Ref. [31] the authors derive an equation for the string centres; we find (see Section 8.5) that its solutions correspond to eigenvectors of H that are not eigenvectors for the Q operator (cf. the following example), and in particular the FM strings obtained following [31] do not correspond to zeroes of the Q function. In this article we present a way that enables us to determine α fm in a way that remedies this. A conjecture for the imaginary part of the FM string centres is given in Section 10.2.
Example. To preview the kind of results we are able to get with our methods consider a homogeneous spin chain with N = 6, φ = 0 and ∆ = 1/2 = cos(π/3) so 1 = 1 and 2 = 3. At M = 3 there are two degenerate states with the same eigenvalues of T s (u), 2s ∈ Z >0 . These eigenvalues moreover coincide with those of |↑↑↑↑↑↑ and |↓↓↓↓↓↓ up to a sign. Using the techniques from Section 6 we find that the two degenerate states have FM strings whose centres we can determine exactly: where the Bethe roots are found analytically using the truncated two-parameter transfer matrix that we will introduce in Section 6. The degeneracies signal some symmetry acquired by the XXZ spin chain at root of unity. Namely, in this case a there is a representation of the loop algebra of sl 2 [37][38][39]. In short: although the 2 th powers of S ± andS ± vanish one can regularise these operators to obtain generators of the loop algebra. Its lowering operators produce eigenvectors that correspond to the FM string (2.21). In fact, the eigenspace is spanned by any two out of the four vectors lim η→iπ/3 That is, the Bethe vectors on the right in (2.22) are nontrivial linear combinations of the two 'loop descendants' on the left. In particular, the latter are not eigenvectors of the two-parameter transfer matrixT(x, y) that we will introduce in Section 6.2, and hence neither of the Q operator.

Basics of the QISM
Let us briefly review the quantum inverse scattering method (QISM). We start with the Lax operator associated to a single spin-1 2 site. Consider an auxiliary space V a , which is an irreducible representation of U q (sl 2 ). Let us work with generators S ± a along with K a = q S z a . In terms of this formulation the commutation relations (2.6) read Although the sites of the spin chain will always have spin 1 2 , we will consider various representations V a , summarised in Appendix A.2. The Lax operator on V a ⊗ C 2 is where the matrix acts on the spin-1 2 representation at site j and the entries are operators on the auxiliary space. Importantly, the Lax operator obeys the RLL relations is a second copy of the auxiliary space. We will parametrise the entries of the R-matrix in such a way that the algebraic Bethe ansatz immediately gives the usual Bethe equations, which requires the shift by η/2 in (3.3). In case the auxiliary space is the spin-1 2 representation the R-matrix which can be expressed as Figure 1: Graphical representation of the transfer matrix with N sites. The (cyclic) thick red line is the auxiliary space and the black lines represent the spin-1 2 spaces. The twist operator is labelled by φ. and the Lax operator (3.2) is given by the same matrix with u u − η/2. 4 Now consider N spin-1 2 sites, each with its own (local) Lax operator L aj (u). The mon- where we included φ through the twist operator E a (φ). We will only consider diagonal twists, so that the R-matrix commutes with E a (φ) E b (φ). When V a is the spin-1 2 irrep we have E a (φ) = diag(e iφ , 1); see Appendix B.1 for the other representations that we will use. It is easy to show (using the 'train argument') that the monodromy matrix obeys the RLL relations (3.3) too. There are several things we get out of this setup.
First of all we can construct the Hamiltonian (2.1) and its symmetries. The transfer matrix T is the trace of monodromy matrix M a over the auxiliary space T(u, φ) = tr a M a (u, φ), (3.6) provided this trace can be taken -which is not obvious when V a is infinite dimensional; we will get back to this. The RLL relations 6 imply that this is a one-parameter family of commuting operators, As a consequence any expansion in u generates a hierarchy of commuting operators that do not depend on u. Particularly important charges are obtained by taking logarithmic derivatives with respect to the spectral parameter u: In particular, when the auxiliary space has spin 1 2 -the transfer matrix of the six-vertex model -this includes the twisted translation operator (j = 1) and the XXZ Hamiltonian (2.1) (j = 2), see Appendix B.2. More generally, when V a ∼ = C 2s+1 is the spin-s irrep of U q (sl 2 ) with s ≥ 1 the resulting conserved charges are quasilocal [42,43,45].
Next, the spectrum of the transfer matrix with V a ∼ = C 2 , and in particular the XXZ spin chain, can be characterised via the algebraic Bethe ansatz. Keeping V a ∼ = C 2 the monodromy 4 When Va is two-dimensional the Lax matrix (3.2) is symmetric in the sense that P L P = L, with P the permutation on C 2 ⊗ C 2 . This allows us to reverse the roles of a and j in (3.2). One obtains the matrix in (3.4) with u replaced by u − η/2 using S ± a = σ ± a and Ka = exp(η σ z a /2) = diag(e η/2 , e −η/2 ). 5 Inhomogeneity parameters can be effortless included as usual. Since we are interested in the spectrum of the homogeneous spin chain we omit the inhomogeneities from the start to keep the notation light. 6 A more general construction of the R matrix is given in Section 4.3, cf. Eq. (4.23). matrix can be written as a matrix in auxiliary space, with entries that act on the spin-chain Hilbert space (C 2 ) ⊗N . Let us stress that this point of view is opposite to that in (3.2), where we were thinking of the Lax operator as a 2 × 2 matrix on the physical space associated to site j, with entries that were operators on the auxiliary space. When V a = C 2 and N = 1 the two perspectives happen to coincide, see Footnote 4 on p. 11. The operators on the diagonal in (3.9) give rise to the (twisted) transfer matrix, The operators in the off-diagonal entries of (3.9) are used for the algebraic Bethe ansatz. Namely, let |Ω = |↑↑ · · · ↑ be the pseudovacuum state, killed by C(u). A routine computation using the RLL relations shows that is an eigenvector of (3.10) with eigenvalue provided the parameters {v m } satisfy the Bethe equations (2.7). The energy and momentum in (2.8)-(2.9) follow from (3.8).
Remark. When the spectrum of transfer matrices of a system at root of unity contains eigenstates associated with FM strings, the construction (3.11) is not enough to obtain all the eigenstates [20,30,31]. However, it is still possible to label the eigenstates as |{u j } M j=1 where {u j } M j=1 are interpreted as the zeroes of the eigenvalue of the Q operator constructed in Section 6.2, see also Sections 5.5 and 6.6. We will use this definition to label all the eigenstates of Q operator. In the absence of FM strings they can be explicitly constructed via the algebraic Bethe ansatz (3.11); for the case with FM strings see our conjectures in Section 9.

Transfer matrices
By varying the dimension of the auxiliary space one obtains different transfer matrices. In general, for a choice of auxiliary space V a , local Lax operator L aj (u) on V a ⊗ C 2 and twist operator E a (φ) on V a one gets a monodromy matrix and corresponding transfer matrix as in (3.5)-(3.6). In (3.9)-(3.10) we encountered the example where V a = C 2 is a spin-1/2 space, yielding the ('fundamental') s = 1/2 transfer matrix of the six-vertex model.
In order to remember which auxiliary space is used in the definition of a particular transfer matrix we will from now on decorate each T (with e.g. a subscript, superscript, tilde) to keep track of the choice of V a that was traced over; here we deviate from the usual meaning of subscripts that we used so far. That is, if V a is some U q (sl 2 ) irrep R then we will write For example, from now on T 1/2 denotes the fundamental transfer matrix (3.9)-(3.10). We will need the following choices of R, whose details can be found in Appendix A.2: • The unitary spin-s representation with s ∈ 1 2 Z ≥0 . Here V s = C 2s+1 . We denote the monodromy and transfer matrices by M s and T s . This generalises the case s = 1/2 considered so far, and leads to the quasilocal charges mentioned above.
• The highest-weight spin-s representation with 'complex spin' s ∈ C. We denote its transfer matrix by T hw s .
-For generic q (not at root of unity) this representation is infinite dimensional. One has to take care that the trace in (3.6) makes sense in this case.
-If s ∈ 1 2 Z ≥0 it can be truncated to a (2s + 1)-dimensional (sub)representation. The result coincides with the unitary spin-s representation up to a gauge transformation (conjugation). In particular, having taken the trace over the auxiliary space, the truncated transfer matrix is equal to T s . See Section 5.1.
-At root of unity both the lowering and raising operators of the complex-spin highest-weight representation become nilpotent, allowing for another truncation to a 2 -dimensional subspace for any s ∈ C. We will denote this truncated transfer matrix byT s . It will be introduced in Sections 6.1-6.2.
Let us already stress an important difference between the two truncated transfer matrices. Since 2 varies wildly as η = iπ 1 2 runs through iπ Q the spectrum of the truncated transfer matricesT s , with 2 -dimensional auxiliary space, will not be continuous in η. This is illustrated in Appendix E.2. Instead, the transfer matrices T s for s ∈ 1 2 Z have spectra that vary smoothly with η.
A generalisation (obtained by fusion in the auxiliary spaces) of the RLL relations (3.3) guarantees that each of the above transfer matrices commute amongst themselves, like in (3.7), as well as with each other (for different s, s ). As a result, they share the eigenvectors produced by the algebraic Bethe ansatz (3.11).
The transfer-matrix fusion relations are a system of equations that show how the higherspin transfer matrices T s can be constructed from T 1/2 . This is usually proven by fusion, taking the tensor product of several spin-1 2 representations (putting several monodromy matrices on top of each other) and projecting onto the spin-s submodule, see e.g. [56,57]. We will get back to the fusion relations in Section 5.3 and 6.4.

Q operators
There have been numerous endeavours to understand the relations between the different transfer matrices and how they can be used to characterise all eigenstates of the XXZ spin chain. The most famous of these was found by Baxter in the 70s [14][15][16], where he constructed the Q operator that satisfies the matrix TQ relation (with respect to T 1/2 ) for the eightvertex model, which is closely related to the XYZ spin chain. A similar construction can be performed for the six-vertex model and the XXZ spin chain [9]. The eigenvalues of the Q operator are called Q functions, and their zeroes are precisely Bethe roots, i.e. solutions to the Bethe equations (2.7), yielding eigenvectors via the Bethe ansatz as in (3.11). Baxter constructed the Q operator by solving the matrix TQ relation directly. For the purpose of numerically obtaining Bethe roots it is much easier to solve the functional TQ relation, i.e. the relation between the eigenvalues of T 1/2 and the Q operator. We summarise how this works in Appendix C.1. At root of unity, however, the functional TQ relation is not enough to determine the full spectrum.
In the following sections we will construct the Q operator explicitly and prove the matrix TQ relation and the transfer matrix fusion relation. We use a new approach that is based on the factorisation and decomposition of transfer matrix T hw s associated to an auxilary space that is an infinite-dimensional complex-spin representation [40]. This construction works for any anisotropy parameter ∆ ∈ R. In the case of root of unity we prove truncated fusion relations for the transfer matricesT s using the same method. These truncated fusion relations will in turn allow us to prove a conjecture of [32,34,41]. This enables us to construct the Q operator explicit at root of unity, and elucidate the exponential degeneracies, which are closely related to the FM string from Section 2.3 but cannot be resolved via the functional TQ relation.

Factorisation of T hw s
We start with an important technical result: the spin s ∈ C highest-weight transfer matrix gives rise to a two-parameter transfer matrix that admits a useful factorisation. In this and the following sections we consider arbitrary q; we will specialise to root of unity later.

Factorisation of Lax operator
The factorisation of the spin-chain transfer matrix has a long history starting with the chiral Potts model [58,59]. In particular it was used for the XXX model [60], for the XXZ chain [61] and for the XXZ chain with nonconservative boundary conditions [40]. We follow here the approach of [40].
Before we start to prove the factorisation of the transfer matrix T hw s we need to study the factorisation of the Lax operator when the auxiliary space V a is an infinite-dimensional complex spin-s representation. Although for T hw s we are interested in the 'half-infinite' highestweight (Verma) module, see (A.9) in Appendix A.2, our starting point is a bigger space. Although the monodromy matrix can be defined using this internal space, its trace cannot be taken as discussed in Appendix B.1 and needs to be truncated. Let V a be the infinitedimensional Hilbert space with orthonormal basis |n a , n ∈ Z. It decomposes as a direct sum where V + a is spanned by |n a with n ≥ 0, which is the space that we are after, while V − a is the span of |n a with n < 0. For this auxiliary space we will implicitly assume that the trace in (3.6) is only over V + a . Let us write |n n | a := |n a a n | for the matrix basis. Consider the following operators: (4.1) If we think of V a as the sequence space 2 by identifying |n a with the sequence δ n with entries (δ n ) i = δ ni then X a is the right shift. The two operators (4.1) form a Weyl pair, W a X a = q X a W a . In other words, W a counts the weight and X a raises it. We will also use the adjoint X † a = n |n n + 1| a . On V a it is the inverse of X a . In particular X † a commutes with X a on V a .
The half-infinite space V + a is preserved by both of (4.1), but not by X † a . The role of the latter is taken over by the adjoint of the restriction X a | V + a , which we denote by On V + a we have Y a X a = 1 while X a Y a = n>0 |n n| a = 1. Together, (4.1) and (4.2) can be used to give a highest-weight representation of U q (sl 2 ) on V + a with spin s ∈ C: on V + a : See also Appendix A.2. Now we turn to the Lax operator L aj (u, s) associated to a site j. Let us introduce two 'spectral parameters' x and y as simple combinations of u and s, x := u + 2s + 1 2 η, y := u − 2s + 1 2 η, (4.4) that will be convenient when deriving the transfer matrix fusion relations. Starting with the big auxiliary space V a the Lax operator can be decomposed into a product of operators separating the dependence on these spectral parameters (see also appendix B of Ref. [62]): where the two by two matrices u(x) and v(y) are defined as The factorisation (4.5) coincides with the one introduced in [61]. To understand (4.5) we compute the product on the right-hand side. The result is where we simplified the top-left entry using W a X a = q X a W a and X † a X a = 1. Importantly, the four matrix elements preserve the subspace V + a . This is obvious for all but the top-right entry, for which the point is that X † a (W a − W −1 a ) |n a = (q n − q −n ) |n − 1 a has vanishing prefactor when n = 0. Thus the effect of the restriction is to replace X † a by Y a . In view of (4.3) and q = e η we recognise the entries of (4.7) as on V + a : This shows that the right-hand side of (4.5) is the same as that in (3.2). The point of this discussion is to show that (4.5) can be used instead of (3.2) when computing the transfer matrix provided that we restrict the trace to V + a . This will be useful in Section 4.2.

Intertwiners
For convenience let us denote the product on the right-hand side of (4.5) by L aj (u x , v y ).
Exchanging u x and v y we instead obtain In the second line the transpose is both in the auxiliary space and in the physical space; note that (on the auxiliary space) X T a = X † a while the diagonal operator W T a = W a is symmetric. This time the matrix elements clearly preserve V − a , which allows us to consider their action on the quotient V /V − a ∼ = V + a . Practically this just means that we treat all |n a with n < 0 as zero. We can thus view (4.9) as acting on V + a by substituting Y a for X † a in (4.9). Then we construct the monodromy matrix as in (3.5) and finally take the trace over V + a as in (3.6) to obtain the transfer matrix. Let us show that the result is the same when we first take the product of (4.9) as in (3.5) and then restrict to V + a to take the trace. Let P + a denote the orthogonal projection of V a onto V + a . Consider P + a times a product of the matrix elements of (4.9). Since the latter preserve V − a we can replace each factor in the product by P + a times that factor without changing the result. Therefore, the restricted trace of such a product in V + a amounts to trace the product of the projected matrix elements and the effect is to replace X † a by Y a in (4.9). It is straightforward to intertwine L aj (u x , v y ) with L aj (v y , u x ) when both are viewed as acting on V + a : where the solution for the intertwiner is  Note that F a (x, y) only depends on x − y = (2s + 1)η, and not on the original spectral parameter u. It is well defined for generic values of this quantity, namely for x−y / ∈ η Z⊕2πi Z.
To verify that F a (x, y) does the job compare (4.7) and (4.9). It is clear that the intertwining relation holds for the entries on the diagonal. For the remaining two entries (which are related by transposition) the relation follows from Now we introduce a copy V b ∼ = V a with its own spectral parameter u and spin s ∈ C, or equivalently x , y as in (4.4). Consider the product L aj (u x , v y ) L bj (v y , u x ). We can construct an intertwiner G ab (y, y ) that interchanges v y and v y in this product: (4.12) To solve this we first consider the big space V a ⊗ V b and write the Lax operators as products like in (4.5). We look for G ab (y, y ) = g y,y (X a X † b ) in the form of a power series g y,y in Z ab := X a X † b . This commutes with X † a on the left and with X b on the right. Further using (4.13) Using (4.6) this reduces to the functional equation g y,y (q z) (1 + z e −y+y ) = g y,y (q −1 z) (1 + z e y−y ), (4.14) which is solved by Since X a X † b preserves the subspace V + a ⊗ V − b the intertwiner does so too. Restricting to V + a and passing to the quotient Similarly, we define an intertwiner H ab (x, x ) which interchanges u x and u x : This time the roles of V a and V b are reversed and we seek a power series in X † a X b . Proceeding as before, now taking the quotient V a /V − a ∼ = V + a and restricting V b to V + b , we find

Two-parameter transfer matrix
Now we reparametrise the transfer matrix with infinite-dimensional complex spin-s auxiliary space as a two-parameter transfer matrix: The important parameters are x, y and u, s, related by (4.4); the twist φ is just a spectator. Due to the existence of the intertwiners this two-parameter transfer matrix satisfies 20) and in particular forms a family of commuting operators, The proofs of (4.20) are routine, using a variation of the argument that establishes (3.7). Since it might nevertheless be instructive to see it in the present setting let us show the first equality in (4.20). We start with site j, for which the preceding implies By the 'train argument' this readily extends to the two-parameter monodromy matrix: Now assume that H ab (x, x ) is invertible; this is true so long as x − x / ∈ −η Z ≥0 ⊕ 2πi Z. We multiply from the left by F a (x, y) −1 H ab (x, x ) −1 F a (x , y) and take the trace over V + a ⊗ V + b . By the cyclicity of the trace the conjugation by F a (x , y) H ab (x, x ) −1 F a (x, y) −1 drops out on the right-hand side, and we arrive at the desired equality.
More precisely, F a is well defined for s / ∈ 1 2 Z ≥0 ⊕ 2πi Z, and we furthermore need x − x / ∈ −η Z ≥0 ⊕ 2πi Z to ensure that H ab is invertible. Thus the preceding argument establishes the first equality in (4.20) only for almost all values of x, y. However, the Lax operator, and therefore the two-parameter transfer matrix, are continuous in these two spectral parameters. Thus the conclusion holds in full generality by continuity. (The situation is analogous in the standard proof of commutativity of ordinary transfer matrices from the RLL relations; there the R matrix is only invertible for almost all values of the spectral parameter.) The exchange of y and y is shown analogously. Let us note that, together, the intertwiners give rise to an R matrix where P ab is the permutation operator between the auxiliary spaces V + a and V + b . Using the properties of the intertwiners one can show that this is the R matrix for which the Lax operator (4.5) satisfies the RLL relation, Via the train argument this gives a direct proof of the commutativity (4.21). A more detailed investigation of the properties of this R matrix is beyond the scope of the present paper.

Factorisation of two-parameter transfer matrix
The property (4.20) implies that the two-parameter transfer matrix T(x, y, φ) can be factorised into two parts that only depend on the spectral parameter x or y, respectively. Namely, where provided that T(x 0 , y 0 , φ) is invertible. We will return to this invertibility (for the truncated case) in Section 7.4. It further allows us to change the value of y 0 at will: note that Let us therefore suppress the dependence on y 0 . A similar argument applies to P x 0 ,y 0 (y, φ), whose dependence on x 0 , y 0 will from now on be suppressed too.
According to (4.20) the operators (4.26) commute with themselves (at different values of the spectral parameter) and with each other: (4.28)

Matrix TQ relation and transfer matrix fusion relation
It is time to harvest the fruits of our labour. We will show that Q from (4.26) is nothing but Baxter's Q operator, satisfying the matrix TQ relation with twist φ, see (5.10). Moreover, P obeys a very similar matrix 'TP relation', see (5.11). In particular, in the periodic case (φ = 0) Q, P are 2 linearly independent solutions of the matrix TQ relation. We will furthermore derive the transfer matrix fusion relations as well as an interpolation formula that expresses the half-integer spin transfer matrix in terms of Q operators. As in the previous section q is arbitrary.

Decomposition of highest-weight transfer matrix
In order to demonstrate that the Q operator in (4.26) is indeed the same as Baxter's Q operator, satisfying matrix TQ relation, we shall use a decomposition [40,60,[63][64][65][66] of the two-parameter transfer matrix T(x, y, φ), when specialising the complex spin to 2s ∈ Z ≥0 . Recall from Section 4.1 that the auxiliary space is spanned by |n for n ≥ 0; this is what we denoted by V + a in Section 4.1. When 2s ∈ Z ≥0 we can decompose it as V + a ⊕ V a , where V + a is the span of all |n with n > 2s, while the finite-dimensional piece V a is spanned by |2s , · · · , |1 , |0 . The latter is certainly preserved by the diagonal operator K a as well as by S − a , whose block-triangular form is shown in Fig. 2. Since 2s ∈ Z ≥0 the operator S + a preserves V a too. Indeed, as [2s − n] q = 0 for n = 2s the coefficient of |2s + 1 2s| a in S + a vanishes: this entry is marked in red in Fig. 2. That is, all of K a , S ± a are of block lower triangular form. The (2s + 1) × (2s + 1) blocks that act on V a differ from the unitary spin-s representation by a simple gauge transformation (conjugation).
The * represent non-zero entries. The square orange (pink) block acts on V + a (resp. V a ), while the rectangular blue block maps V + a to V a . Note that we order the basis decreasingly, · · · , |1 , |0 , cf. Appendix A.2. Now we go towards the transfer matrix. Since the Lax operator is built from K a , S ± a , see (3.2), for 2s ∈ Z ≥0 it assumes a block lower triangular form with respect to the decomposition V + a ⊕ V a too. Let us indicate its block structure, paralleling that in Fig. 2, by Here we can think of L a j as a square infinite matrix acting on V + a , L a j as a square matrix on V a , and (for want of a better notation) L a j as a rectangular matrix sending V + a to V a ; all with entries that are operators acting at site j. For the monodromy matrix we consider more sites. Note that the blocks on the diagonal only 'talk' amongst themselves: Thus the monodromy matrix inherits the block triangular form By taking the trace we obtain the transfer matrix Since the unitary spin-s representation differs from that on V a by a gauge transformation, tr a M a is nothing but the transfer matrix T s for spin s ∈ 1 2 Z ≥0 . As V + a ∼ = V + a moreover tr a M a is another complex-spin highest-weight transfer matrix! Accounting for the twist and the correct value of the new complex spin we arrive at the decomposition

Generalised Wronskian and matrix TQ relation
By (4.4) and (4.19) we can rewrite (4.25) as In these terms the decomposition (5.5) reads This is the generalised Wronskian relation. Let us consider some examples. For s = 0 the operator T 0 (u, φ) is a scalar, which is independent of the twist according to our choice of the latter (see Appendix B.1). Thus it can be denoted as T 0 (u) = sinh N (u). We therefore obtain the Wronskian relation Note that T 0 (u) is independent of the twist φ while Q and P in the right-hand side do depend on φ. When s = 1/2 we obtain a relation for the fundamental (six-vertex) transfer matrix: Multiplying both sides by Q(u, φ) and using (5.8) to get rid of the P operator we find We have recovered Baxter's matrix TQ relation [9]! If we instead multiply by P(u, φ) and use (5.8) to eliminate the Q operator we analogously obtain a matrix 'TP relation' Note the different positions at which the twist e iφ appears in (5.10) and (5.11). More precisely, the TP relation for the rescaled operator e iφu/η P(u, φ) is nothing but the TQ relation with opposite twist −φ. In the periodic case, φ = 0, the Q and P operators are two solutions to the matrix TQ relation (5.10) that are linearly independent as their Wronskian, the right-hand side of (5.8), is nonzero (T 0 = 0).

Transfer matrix fusion relations
Interestingly, the decomposition of the two-parameter transfer matrix T(x, y, φ) further allows us to derive the transfer matrix fusion relations [56,57,67,68]. As the name suggests these are typically derived by fusion in the auxiliary space (tensoring auxiliary spaces and projecting onto any irrep via suitably related spectral parameters). In the present setup the simple derivation goes as follows. 7 We multiply both sides of (5.7) from the left by T 1/2 (u± 2s+1 2 η, φ) and apply (5.10)-(5.11). Collecting terms with the same T 0 and using (5.7) for each we obtain These two equations, one for the upper signs and one for the lower signs, are the transfer matrix fusion relations. Taken together for all half-integer values of the spin s the functional form of these relations, i.e. the analogous relations for the eigenvalues T s , comprises a system of difference relations called a T-system. Together with a reformulation known as the Y-system it is of vital significance for physical applications like the thermodynamic Bethe ansatz. See Ref. [69] for a thorough review and further references.

Interpolation formula
When s ∈ 1 2 Z ≥0 the transfer matrix T s can be expressed in terms of Q; let us derive this in our framework. Rewrite (5.7) as The meaning of the fractions for multiplication with the inverse is unambiguous since the operators involved all commute. Specialising to s = 0 we likewise rewrite (5.8) as We translate the arguments of (5.14) by kη, multiply by e ikφ and sum over k from 0 to 2s. The telescoping sum on the right-hand side yields the right-hand side of (5.13). In this way we obtain the interpolation-type formula [56,57] T s (u) = Q u + 2s + 1 2

Structure of the eigenvalues of Q and P
Let us now study the properties of the Q and P operators, and their eigenvalues, on their respective spectral parameters. It is convenient to use multiplicative spectral parameters r := e x and t := e y . As always we start from the Lax operator, which we here write as We start with the dependence on t. Let us say that a Laurent polynomial f (t) is a 'trigonometric polynomial of degree n' if t n f (t) is a polynomial in t 2 of degree n. From (4.7) we see that L 11 a , L 21 a are trigonometric polynomials of degree one in t while L 12 a , L 22 a are independent of t. It follows that on a vector with S z = N/2 − M the two-parameter monodromy matrix M a (x, y, φ) acts by a matrix whose entries are trigonometric polynomials of degree N −M in t. Thus the same holds for the two-parameter transfer matrix T(x, y, φ). Regarding r we see that L 11 a , L 12 a are both degree zero while L 21 a , L 22 a are trigonometric polynomials in r of degree one. This means that when acting to the left on (the dual of) the M -particle subspace, M a (x, y) has entries that are trigonometric polynomials in r of degree M . As before this carries over to T(x, y, φ). Since the latter moreover preserves the value of M this remains true when acting to the right as well. The upshot is the operator T(x, y, φ) acts on the M -particle subspace, with S z = N/2 − M fixed, by a matrix with entries that are trigonometric polynomials of degree M in r and degree N − M in t. In terms of the Q and P operators the conclusion is that, on this subspace, Q(x, φ) acts by a matrix consisting of trigonometric polynomials in r of degree M , and P(y, φ) likewise by a matrix of trigonometric polynomials in t of degree N − M . These properties are inherited by the eigenvalues since the (joint) eigenvectors are independent of r, t in view of the commutativity (4.28).
Let us make this a little more concrete. By the commutativity of all transfer matrices the joint eigenvalues of the Q and P operator are given by the algebraic Bethe ansatz. For any on-shell Bethe state, i.e. (3.11) -see also the remark following that equation -subject to the Bethe equations (2.7), we have Baxter realised that the TQ equation determines the eigenvalues Q in familiar terms. Namely, by the above it is a trigonometric polynomial of degree M in t := e u . Let us denote the zeroes by t m : According to the matrix TQ relation (5.10) the eigenvalues obey the functional TQ relation Picking u so that t = t m is a zero of Q, the left-hand side of (5.19) vanishes. Recalling that T 0 (u) = sinh N (u) the result reduces to the Bethe equations (2.7) in multiplicative form, allowing us to identify t m = e um as the multiplicative version of the Bethe roots. That is, the zeroes of the eigenvalues of the Q operator are precisely the Bethe roots [9]. This observation is used in the numerical recipe for finding the Bethe roots in Appendix C.1. In the presence of any Bethe root at infinity, which does occur for XXZ as reviewed in Section 2.2, the form of the eigenvalues has to be modified a little. Indeed, Therefore, Bethe roots at infinity show up in the eigenvalues of the Q operator: if we rearrange the u j so that the infinite roots are last then One can similarly show that the eigenvalues of the P operator P are of the form (5.21) but with M replaced by N − M . The functional version of the TP relations (5.11) gives rise to Bethe-type equations for the N − M zeroes of the eigenvalues P (v, φ) for an M -particle eigenvector: sinh These are precisely the Bethe equations for a Bethe state |{v n } N −M n=1 of the XXZ model with opposite twist −φ. Note that the algebraic Bethe ansatz (3.11) only uses the pseudovacuum |↑ · · · ↑ and the B operator B(u), which is independent of the twist φ, cf. (3.9). Therefore, the off-shell Bethe state |{v m } N −M m=1 does not depend on the twist either; the latter only enters on shell, i.e. upon imposing the Bethe equations. When we impose (5.22) the Bethe vector is not an eigenstate of T s (u, φ), but rather of T s (u, −φ). In particular, in the periodic case (φ = 0) it can be interpreted as a Bethe vector beyond the equator (if M < N/2, so that N − M > N/2). A detailed calculation is presented in Appendix C.3.

Truncated transfer matrix at root of unity
Now we specialise the anisotropy parameter to a root of unity (2.3). Here the infinitedimensional auxiliary space V + a has a finite-dimensional subspaceṼ a , whose size depends on the root of unity, and which is preserved by U q (sl 2 ). Truncating toṼ a allows for another decomposition of the two-parameter transfer matrix that leads to a proof of a conjecture of [32,34,41] and to truncated fusion relations.
Importantly, the truncation enables us to construct the Q operators explicitly, only using finite-dimensional matrices at all intermediate steps. In practice we can do this for all eigenvectors a spin chain with N ≤ 16 and thus obtain the full spectrum of the XXZ spin chain with arbitrary twist φ, revealing the conditions for the appearance of exponential degeneracies that have been observed before [20]. As we will see later this has significant consequences for the thermodynamic limit.

Truncation and intertwiners at root of unity
At root of unity η = iπ 1 2 the matrix elements of the Lax operator (4.7) in the auxiliary space acquire the periodicity: where we recall that ε = q 2 = e iπ 1 . As the periodicity suggests, only a finite part of the Lax operator is really relevant: we can truncate the auxiliary space to a finite-dimensional subspace. This goes via a variation of the construction from Section 5.1, as follows.
a is the span of |n a with n ≥ 2 andṼ a be the span of |n a for 0 ≤ n ≤ 2 − 1. At root of unity we have [ 2 ] q = 0 so one of the entries of S − a vanishes as illustrated in Fig. 3. This time all generators of U q (sl 2 ) preserve the infinite-dimensional subspace V + a . We are interested in the finite-dimensional subspaceṼ a . Like in Section 4.2 we can get there by taking the quotientṼ a ∼ = V + a /V + a . In this way we get a finite-dimensional representation of U q (sl 2 ) onṼ a , see also Appendix A.2. More concretely, all of K a , S ± a are block upper triangular with respect to the decomposition V + a = V + a ⊕Ṽ a , see again Fig. 3. This property is inherited by the Lax operator where L a j can be viewed as an infinite square matrix on V + a , L a j as an 2 × ∞ rectangular matrix mappingṼ a to V + a , andL a j as an 2 × 2 matrix onṼ a . The entries of each of these Figure 3: The decomposition of S − a (left) and S + a (right) at root of unity with 2 = 3, where * represents non-zero elements of the matrices. As in Fig. 2 the bottom-right entry corresponds to |0 0|. are 2 × 2 matrices acting at site j. The truncation to the 2 -dimensional spaceṼ a amounts to treating all |n a with n ≥ 2 as zero.
Now focus on the 2 -dimensional auxiliary spaceṼ a , where we drop the double prime. To describe the truncated Lax operatorL aj more explicitly definẽ Replacing all operators on the auxiliary space in the factorised formula (4.5) by these truncated versions yields an expression with the same structure as (4 .7): which coincides with the right-hand side of (3.2). This time, however, all four matrix entries are finite matrices, with size depending on the value 2 of the root of unity. Notice that the truncation allows us to keep s ∈ C arbitrary. The truncated intertwiners arẽ The intertwining relations are as before, where the truncationL aj (v y , u x ) =L aj (u x , v y ) T tõ V a arises by restriction (rather than taking a quotient).

Decomposition of two-parameter transfer matrix at root of unity
The truncated Lax operator (6.4) can be used to construct monodromy matrices that yield truncated two-parameter transfer matricesT(x, y) =T s (u, φ) where x, y were defined in (4.4).
These truncated two-parameter transfer matrices coincide with the "Q s (exp ηu)" obtained by Korff in Ref. [34], who conjectured that they form a two-parameter family of commuting matrices. Since s ∈ C is unrestricted, the spectral parameters x, y, which were defined in (4.4), are now really just two free parameters. Using (6.5) as before we see that it obeys exchange relations just as in (4.20): and in particular forms a family of commuting operators, As in Section 4.4 the truncated two-parameter transfer matrix therefore factorises as where the truncated Q and P operators are defined as We once more suppress the dependence on x 0 , y 0 ; we will get back to this in the remark in Section 7.4. In terms ofT s (u, φ) the factorisation takes the form (6.10)

Truncated Wronskian and TQ relations
Repeating the arguments from Section 5.1 for an arbitrary complex spin s ∈ C we have the decomposition of the highest-weight transfer matrix at root of unity where up to the twist factor, T hw s− 2 (u, φ) coincides with the matrix T hw s (u, φ) restricted to the basis |m + 2 with m ≥ 0. With (6.1), we simplify the decomposition (6.11) intõ When 2s ∈ Z ≥0 , the decompositions (5.5) and (6.12) yield a decomposition of T s in terms ofT s : This argument relies on the convergence of the trace defining T hw s (u, φ) and requires |e iφ | < 1. In Appendix D we give another proof of this relation valid for an arbitrary twist. From (6.10) we obtain the truncated Wronksian relation (6.14) the left-hand side of (6.14) vanishes. Comparing (6.15) to the conditions (2.15) for the existence of Bethe roots at infinity we see that when (6.15) is satisfied there exist certain numbers M of down spins for which (2.15) is satisfied too.
Proceeding exactly as before we obtain TQ and TP relations that look the same as (5.10)-(5.11), now involving truncated matrices and

Truncated fusion relations
Proceeding as in Section 5.3, but using (6.10) instead of (5.7), we readily obtain fusion-like relations forT s : In analogy to the cases of general q, we will call these the truncated fusion relations. However, we stress that in the present case the internal auxiliary spaces have the same dimensions 2 for all truncated transfer matrices in (6.18), unlike for the fusion relations (5.12).

Interpolation formula: proof of a conjecture
Fabricius and McCoy [32], Korff [34], and more recently De Luca et al. [41] -see Eq. (S22) in the supplementary material of Ref. [41] -conjectured that the complex spin transfer matrix eigenvalues can be expressed in terms ofQ, similarly for the situation of half-integral spin representations from Section 5.4. In our notation, after introducing the dependence on the twist, the formula reads We can prove this using arguments like those leading to (5.15). Rewrite (6.10) as where we introduce a common factor on both sides for convenience. We also cast (6.14) specialised to s = 0 in the form Translating u by (k − s)η, multiplying both sides by e ikφ and summing over k from 0 to 2 − 1, we get From the parity of Q and P -see the discussion following (5.16) -it follows that the second ratio in the second line simplifies to ε N times the first ratio in that line. As long as the twist is not commensurate, i.e. does not obey (6.15), we can cancel 1 − ε N e i 2 φ on both sides of the equation. In fact, since the two sides of (6.22) are continuous in the twist, the result holds for any twist. (An alternative proof, that does not rely on this cancellation or completeness, can be obtained using the decomposition (6.11) like before rather than the simplification (6.12).) So we have Multiplying both sides byQ u + 2s+1 2 η , and using (6.10) we arrive at (6.19).

Structure of eigenvalues ofQ andP
Similar to the discussion in Section 5.5 the eigenvalues of the truncated Q operator, can be expressed asQ with v m = log t m obeying the Bethe equations (2.7) due to the truncated TQ equations (6.16).
In particular the eigenvalues Q(u, φ) andQ(u, φ) of a given eigenvector have the same zeroes.
Since we are at root of unity, the eigenvalues of the truncated Q operator on the M -particle sector are quasiperiodic too:Q Likewise, we can show that the eigenvalues of the truncated P operator are of the form (5.21) but with M replaced by N − M , and

Connection to the work of Frahm et al
Frahm, Morin-Duchesne and Pearce [70] observed that the zeroes of the eigenvalues of the Q operator appear as part of zeroes of eigenvalues of the transfer matrix whose argument is shifted as T ( 2 −1)/2 (u + 2 η/2, φ) at root of unity. In our framework this is clear too.

Applications to XXZ at root of unity: general results
Now we turn to applications of the general formalism developed in Sections 3-6 to the spectrum of the XXZ spin chain (and the six-vertex model's transfer matrix) at root of unity. In this section we use the formalism to derive several general results. This and other features of the spectrum at root of unity will be illustrated by numerous explicit examples in Section 8.

Preliminaries
Let us first define a few useful concepts. Recall that the Q operator and all transfer matrices commute with each other. We will call a joint eigenvector of this family of operators a state. We will assume that the Bethe ansatz is complete, so that any state is of the form |{u m } M m=1 , see (3.11) and subsequent remark, with {u m } M m=1 a solution to the Bethe equations (2.7). Consider for a moment the periodic isotropic Heisenberg XXX spin chain. Here infinite rapidities (vanishing quasimomenta) can be added to the Bethe roots to get sl 2 -descendants (see e.g. §1.1.3 in [10]). Unlike when a finite root is added, infinite rapidities do not change the values of the other Bethe roots, and the eigenvalues for the transfer matrix stays the same as well. Likewise, for the XXZ model at root of unity there are states |{u m } M m=1 for which certain special roots can be added to get another state without affecting the values of {u m } M m=1 or eigenvalues. We will call states that are minimal in this sense 'primitive'; it is similar to a highest-weight condition except that we do not use representation theory to characterise it. Namely, we call a state |{u m } M m=1 primitive if no nontrivial subset of {u m } M m=1 corresponds to a physical state whose eigenvalue for T s differs at most by a sign. (For the reason why we allow for a sign see Section 7.2.) Primitive states have no FM strings. Recall that by n ±∞ we denote the numbers of Bethe roots at ±∞, respectively. Typically, namely for twist φ / ∈ {0, π}, a primitive state has no roots at infinity either. 8 In the (anti)periodic case φ ∈ {0, π} a primitive state may have n ±∞ = 0 as long as n ∓∞ = 0, see (2.18) and Section 8.1.3. 8 States with n±∞ = 0 but n∓∞ = 0 do occur at twist φ / ∈ {0, π}. In Section 8.2 we will show that they have the same Ts-eigenvalues as certain primitive states at twist −φ that have the same finite Bethe roots, and of which they should be considered descendants.
Let us remark that it is possible for two primitive eigenstates to be degenerate, namely when N is odd and the two states are related by the spin flip operator j σ x j : we will discuss this case in Section 7.4.
Any state that is not primitive is a descendant of some primitive state: it satisfies where |{u m } M m=1 is primitive and the additional Bethe roots {w n } M n=M +1 consist of FM strings and roots located at ±∞. Away from the isotropic points descendant states only exist when condition (6.15) is satisfied.
Finally, for φ / ∈ {0, π} away from the (anti)periodic points, a primitive state and its descendants might be eigenstates for Hamiltonians with opposite twist, φ and −φ, see Appendix C.3. The sign of the twist in T s (u, ±φ) is fixed accordingly. We will further illustrate this in Section 8.2.
Remark. We have used the term 'descendant' here, because the corresponding 'descendant tower' that we will introduce in Section 8 resembles the descendant structures of Lie algebras. We use 'primitive' rather than 'highest weight' to stress that we do not use representation theory to characterise these states. (One could similarly use 'derived state' rather than 'descendant', but we will refrain from doing so. We hope that this will not cause too much confusion.) In fact, primitive states do have highest weight for the loop algebra of sl 2 [39], cf. [32]. From our perspective this algebraic perspective is not completely satisfactory since the loop-algebra action does not commute with the two-parameter transfer matrix or Q operator: as the example in Section 2.3 illustrates, the loop-algebra descendants are no longer eigenvectors of the latter. We do believe that there exist a deeper algebraic structure that is compatible with the Q operator, which would allow for unambiguous use of the terms 'highestweight' and 'descendant'. We intend to investigate the underlying algebraic structure in the future.

Impact of FM strings on transfer-matrix eigenvalues
With this terminology let us show that a descendant state |{v m } M m =1 of a primitive state |{u m } M m=1 has T s (u, ±φ)-eigenvalues that differ by at most a sign. This means that the momentum of the primitive state and its descendant may differ by π while all other (quasi-)local charges generated by T s (u, ±φ) are identical.
Write the eigenvalues of T 1/2 and the Q operator for the primitive state |{u m } M m=1 as For simplicity we assume that |{u m } M m=1 does not contain any Bethe roots at ±∞. By Section 5.5 the eigenvalue of the Q operator is a trigonometric polynomial of degree M , and satisfies the functional TQ relation is also an eigenvector of the Q operator, where for φ / ∈ {0, π} the sign of the twist ±φ has to be chosen to match that of the Hamiltonian, see Section 8.2. From the definition (2.19) of an FM string we note that the eigenvalue of the Q operator on the descendant state is where n ±∞ and n fm are the numbers of roots at ±∞ and FM strings, respectively, of the descendant state. Observe that Q (u) satisfies the TQ relation where we recall that ε = q 2 ∈ {±1}. Hence, the eigenvalue of T 1/2 for the descendant state Since the T 1/2 -eigenvalue of the descendant state differs by at most a sign from that of the primitive state, the descendant has the same local charges generated by T 1/2 as the primitive state, except that its momentum might differ by π. More generally, from the transfer matrix fusion relations (5.12) we obtain (7.9) We will give several explicit examples of the descendant states associated with primitive states, and the descendant towers that they form, in Sections 8.1-8.2.

FM strings at commensurate twist
In Section 6.3 we saw that, at root of unity, when the twist φ is commensurate as in (6.15) the truncated Wronskian relations (6.14) trivialise in the sense that the eigenvalues ofQ and P are proportional to each other. This can only be achieved when the eigenvalues ofP are related to those ofQ through adding or subtracting Bethe roots at infinity or FM strings, cf. (6.14): eigenstates with FM strings must occur when the twist is commensurate. Conversely, if FM strings are present amongst the zeroes of the eigenvalues ofQ orP the truncated Wronskian relations (6.14) vanish, since for s = 0 the left-hand side of (6.14) does not contain zeroes of FM strings while the right-hand side does. This means that descendants containing FM strings can only exist when the twist φ is commensurate. In conclusion, FM strings can only, and necessarily do, occur at commensurate twist. From (6.12), wheneverT s (u, φ) has nonzero eigenvalues it follows that the eigenvalues T hw s (u, φ) must have poles if 1 − ε N e i 2 φ = 0 to compensate the vanishing prefactor. This corroborates Refs. [25,28] in which it was shown that terms proportional to log t can arise when solving the functional TQ relations, yielding eigenvalues that are quasi-polynomials in t [28]. As it turns out, the appearance of these quasi-polynomials is also closely related to the FM strings and their string centres, cf. (6.15). A detailed discussion of this will appear in a forthcoming article [71]. When (6.15) is satisfied, a quantisation condition for the centres of FM strings based on (6.23) was obtained in Refs. [32,34,36]. Let us review their arguments. For an eigenstate where the eigenvalues are trigonometric polynomials in t = e u . Let us decompose the Q function in to a 'regular' part Q r (u), consisting of n r factors whose zeroes are Bethe roots that are neither FM strings nor infinite, a 'singular' part Q s (u), consisting of n fm FM strings, along with factors accounting for Bethe roots at infinity as in (5.21): where the number of down spins M = n r + 2 n fm + n +∞ + n −∞ . The proportionality signs in (7.11) indicate equality up to t-independent factors. Because the Wronskian (6.14) vanishes, the eigenvalue P (u) has the same regular zeroes as Q(u): wheren ±∞ andn fm are the number of Bethe roots at ±∞ and FM strings, respectively, present amongst the zeroes of the P operator. The total number of zeroes of Q(u) and P (u) is equal to the system size N : N = 2 n r + 2 n fm + 2nfm + n +∞ + n −∞ +n +∞ +n −∞ . (7.13) Consider the functional form of the interpolation formula (6.23), (7.14) Using (7.11)-(7.12) and the fact that Q s (u + η) = ε nfm Q s (u) this can be rewritten as e ikφ T 0 u + (k + 1/2) η Q r (u + k η) Q r u + (k + 1) η e i(2k+1)(n −∞ −n +∞ )η .

(7.15)
This is a 'quantisation condition' for the string centres of FM strings. Indeed, (7.15) implies that for any two states belonging to the same descendant tower, i.e. sharing the same regular part of their Q functions, the combination Q s (u) P s (u) contains the same zeroes. In other words, the string centres of FM strings for any state belonging to a descendant tower are determined solely by the regular part of the Q functions. Moreover, string centres of FM strings are free within a descendant tower in the sense that adding or removing FM strings whose string centres are given by the zeroes of Q s (u) P s (u) results in other eigenstates within the same descendant tower. We can add FM strings from the zeroes of Q s (u) P s (u) to the primitive state in order to generate the descendant states. We present explicit examples of the resulting tower structures in Sections 8.1-8.2.

Primitive degenerate eigenstates
When the system size N is odd and condition (6.15) for a commensurate twist is satisfied, it is possible to have two primitive eigenstates with degenerate eigenvalues of T s [72]. Namely, the two eigenstates are related by reversing all spins, When this happens, the eigenvalues ofT s are zero for both states, In this case, the matrixT(x, y, φ) =T s (u, φ) can not be inverted for any x, y ∈ C. Despite this, the decomposition (6.20) still applies, where rather than via (6.9) the P operator is defined through the TQ relation for the state beyond the equator. However, it is no longer possible determine the eigenvalues of the Q operatorQ for the degenerate primitive states. In practice this is not a problem since, as both states are primitive, we can use the numerical recipe in Appendix C.1 to determine the zeroes of the Q functions (Bethe roots) for both states by solving the functional TQ relation numerically.
Curiosity at supersymmetric point. In general the number of degenerate primitive eigenstates at root of unity and commensurate twist increases as the (odd) system size N grows. The exception to this rule is the special point η = 2π 3 , φ = 0 (or η = π 3 , φ = π, cf. Appendix C.2), for which there are only two degenerate primitive eigenstates for any odd N . (There are no degenerate primitive eigenstate when N is even.) These values of η correspond to the supersymmetric point ∆ = − 1 2 . For commensurate twist φ = 0 and odd system size N , the antiferromagnetic ground states are always doubly degenerate, and have been well studied. We call them Razumov-Stroganov (RS) states, from the conjecture made in Ref. [73] and later proven in Ref. [74]. RS states are closely related to lattice supersymmetry [75,76]. (Note that in our convention for the Hamiltonian the RS states are the highest excited states in the spectrum.) Interestingly, RS states are the only two primitive states that have degenerate eigenvalues for T 1/2 at ∆ = −1/2, for any odd N . Let us denote them by |RS 1 , in the sector with M = (N − 1)/2 down spins, and |RS 2 = N j=1 σ x j |RS 1 , with M = (N + 1)/2. Their eigenvalues for T 1/2 are The Bethe roots can be obtained using numerical method in Appendix C.1. For instance, for N = 5 we have and (7.20) One naturally wonders what the structure of the truncated two-parameter transfer matrix T(x, y, φ) is at ∆ = 1/2, and how it relates to the RS states. We will postpone this question to future work.

Q functions for fully polarised states
Since all transfer matrices and the Hamiltonian commute with the total magnetisation S z , the two fully polarised states |↑↑ · · · ↑ and |↓↓ · · · ↓ are eigenstates of all transfer matrices, and therefore of the Q and P operators. Let us study this in more detail.
For the pseudovacuum |↑↑ · · · ↑ on top of which magnon excitations are built it is easy to see that for any system size N and twist φ, with Q function that does not depend on u or φ. From the definition of the two-parameter transfer matrix (3.6) and (6.4) we havẽ where the Q function is which is consistent with (7.15). (Note that one needs to compare Q ⇓ (u, φ) with P ⇑ (u, −φ).) From the final expression in (7.23) we see that and since ε = q 2 = ±1, we conclude that if 1 − ε N e i 2 φ = 0 then Q ⇓ (u + η, φ) = e iφ Q ⇓ (u, φ). (7.25) This is precisely the commensurate-twist condition (6.15) for the appearance of FM strings. When φ = 0 and the condition (6.15) is satisfied Q ⇓ (u, φ) is a trigonometric polynomial in z = t ξ 2 = e ξ 2 u with ξ = 1 (ξ = 2) if 1 is even (odd). In that case the Bethe roots, i.e. the zeroes of (7.23), consist of only FM strings and pairs +∞, −∞.

Applications to XXZ at root of unity: examples
The physical Hilbert space (C 2 ) ⊗N consists of all joint eigenstates of the transfer matrices.
Assuming the completeness of the Bethe ansatz these states can be distinguished by their Bethe roots. Degeneracies are known to occur amongst the eigenvalues of the transfer matrices T s when the condition (6.15) for commensurate twist is satisfied [20,30,31]. At first sight, this might resemble the degeneracies due to the SU (2), or sl 2 = (su 2 ) C , symmetry at the isotropic point. Yet such degeneracies are not expected away from the isotropic point. In this section we will show with several concrete examples how to construct the descendant towers (Hasse diagrams) that link all different eigenstates with the same eigenvalues (possibly up to a sign for T 1/2 , cf. Section 7.2), for the transfer matrices. Moreover, as was observed by Fabricius and McCoy [32], we will illustrate in Section 8.1.1 that these degeneracies grow exponentially. This has significant consequences in the thermodynamic limit, which will be discussed in Section 10.2. The discussions in this section are less rigorous than the previous sections. We have used the following method. In order to study the spectrum of XXZ spin chain at root of unity we need to know the Bethe roots associated to the eigenstates of the model, which is equivalent to knowing the Q functions, i.e. the eigenvalues of the Q operator, for these eigenstates. Making use of the decomposition of the truncated two-parameter transfer matrixQ(x, φ) =T(x, 0, φ), we construct the Q operator for the XXZ spin chain at root of unity explicitly for system size N ≤ 16. This is possible thanks to the truncation of the auxiliary space at root of unity. All the explicit examples in Section 8.1-8.2 are obtained through this procedure. In many instances, analytic expressions can be obtained using symbolic algebra software.

Descendant towers in periodic case
Let us first consider twist φ = 0 vanishes, i.e. periodic boundary conditions, and illustrate how descendant states can be found from a primitive state together with FM strings or pairs +∞, −∞. The result will be a descendant tower, consisting of all eigenstates with degenerate (possibly up to a sign for s = 1/2) eigenvalues of T s for all 2s ∈ Z >0 . At φ = 0 every eigenstate |{u m } M m=1 } (M < N 2 ) has at least one eigenstate with the same eigenvalue for T s , namely the spin-flipped state beyond the equator, We will refer to |{u m } M m=1 } and |{v m } N −M m =1 } as states on the opposite side of the equator with respect to each other. As we will illustrate, these two states are closely related when (6.15) is satisfied, forming the top and bottom of a descendant tower that includes various intermediate states.
It is easy to verify (7.25). The zeroes of this Q function 9 are of the form , α fm 4 = α 4 + i Now let us for a moment consider the free-fermion case ∆ = 0 (η = iπ/2), φ = 0, which was an important example for [37]. In that case the S -matrix reduces to S = −1, so the Bethe equations decouple with Bethe roots taking values of the form α n ± iπ/2 with α n ∈ [0, 2π). Reality of the energy requires any solution to consists of complex-conjugate pairs α n + iπ/2, α n − iπ/2. The very simple resulting spectrum can be understood as consisting of FM strings of length 2 = 2, and the free-fermion nature is visible in that any string can be added or removed without affecting the location of the other strings. See e.g. [36] for the details. More generally, for any root of unity, (7.15)  This observation immediately yields the following descendant tower. Given a primitive state, construct the corresponding eigenstate beyond the equator and find its Bethe roots as above. The descendant tower is obtained from the primitive state by adding FM strings one by one, with magnetisation M jumping by 2 each time, until we reach the eigenstate beyond the equator at the bottom. For example, the descendant tower obtained in this way from the pseudovacuum that we considered above is illustrated in Fig. 4. This structure is easily verified by explicitly constructing the Q operator, and one sees that the eigenvalues of the transfer matrices T s are degenerate (up to a sign) for all descendant states.
In particular it follows that the number of descendant states within the magnetisation sector with number of down spins fixed to M = M 0 + n 2 is equal to n fm n , i.e. grows binomially in the number n fm of FM strings for the state beyond the equator. By adding up all descendant states at the occurring values of M we find that the total number within the descendant tower is n total = nfm n=0 n fm n = 2 nfm , (8.6) in agreement with the loop-algebra prediction, cf. Eq. (1.19) of Ref. [32]. The the descendant tower is exponentially large in the number n fm of FM strings present in the state beyond the equator corresponding to the primitive state. Figure 4: Illustration of a descendant tower from the pseudovacuum for N = 12, ∆ = 1/2, φ = 0. Arrows of different colours correspond to the addition of different FM strings.

Descendant towers with pairs of roots at infinity
Next we turn to an example of a descendant tower that is generated by adding FM strings as well as pairs +∞, −∞ to the Bethe roots of a primitive state. Recall that Bethe roots at infinity correspond to applications of the lowering operators of U q (sl 2 ) (see Appendix A.1). At root of unity 2 applications of either of S ± orS ± already gives zero, so the total numbers of Bethe roots at +∞ and −∞ must be smaller than 2 : we need n ±∞ < 2 . As we will illustrate, pairs of Bethe roots at +∞, −∞ play a similar role as FM strings. We continue with the example N = 12, ∆ = 1 2 (η = iπ 3 ), φ = 0. Let us now start from a primitive state with one down spin (M 0 = 1), |{u 1 } , say for the solution to Bethe equations (2.7) given by Consider again the corresponding state beyond the equator, |{v m } 11 m =1 ∝ 12 j=1 σ x j |{u 1 } .
We obtain its eigenvalue for the Q operator by solving (7.15), yielding the Q function The zeroes of this Q function are u 1 , two FM strings, along with (cf. Section 5.5) two pairs of roots +∞, −∞: Here the FM strings are centred at with real parts that have numerical values satisfying the quantisation condition (7.15). Note that n ±∞ = 2 < 2 = 3. Unlike for FM strings, Bethe roots at infinity cannot added one by one: {u 1 , ±∞} does not satisfy the Bethe equations, as it violates the condition (2.12). On the other hand, using (A.7) we can easily check that is an eigenstate with the same eigenvalues for T s as |{u 1 } . That is, a pair of Bethe roots +∞, −∞ can be viewed as a 'bound state': these two roots always appear together in order to satisfy condition (2.15). 10 Like FM strings, pairs +∞, −∞ do not scatter with other magnons. Therefore, when the corresponding descendant state beyond the equator contains one pair of infinite Bethe roots as in our example, the total number of states in the descendant tower is as illustrated in Fig. 5. Almost all (primitive) states at M = 1 give rise to a descendant tower of this form, only differing in the locations of the FM string centres. There are two exceptions, to which we turn next.

Mirroring descendant towers
There is one more interesting feature present in the spectrum for N = 12, ∆ = 1 2 (η = iπ 3 ), φ = 0. At M = 1 down spin we find from (2.15) that there are two eigenstates, namely |{+∞} and |{−∞} , which have an infinite Bethe root. These are primitive: their T s -eigenvalues are different from those of the pseudovacuum, reflecting that we are not at the isotropic point. We will concentrate on |{+∞} in this section; the situation for |{−∞} is analogous.
In general, in the (anti)periodic case φ ∈ {0, π}, the existence of an eigenstate of the form |{u m } M m=1 ∪ {n × +∞} for some n > 0 implies the existence of another eigenstate |{u m } M m=1 ∪ {( 2 − n) × −∞} with the same eigenvalues of T s (up to a sign), whose Bethe roots can be found from the Bethe equations (2.7). (Likewise, the presence of a state |{u m } M m=1 ∪{n×−∞} implies that of |{u m } M m=1 ∪ {( 2 − n) × +∞} .) One can verify that the state beyond the equator corresponding to |{+∞} does not belong to the same descendant tower as |{+∞} , unlike for the examples presented in the previous sections. The eigenvalue for the Q operator on the corresponding state beyond the equator, |{v m } 11 m =1 } ∝ 12 j=1 σ x j |{+∞} , is obtained through (7.15), yielding (8.14) Its zeroes contain three FM strings and two Bethe roots at −∞, where the FM strings have centres α fm n = α n + i π 6 , α 1 = −0.404723313, α 2 = 0.075767627, α 3 = 0.613080368, (8.16) which satisfy the quantisation condition (7.15). The resulting descendant-tower structure is as follows. The primitive state |{+∞} gives rise to descendants via the addition of the above FM strings. This yields a descendant tower that accounts for half of the states with the same eigenvalues (possibly up to a sign) of the transfer matrices T s . The other half of the states form a 'mirroring tower', obtained from the first tower by spin reversal. At the top of the mirroring tower we find the primitive eigenstate |{2 × −∞} , whose eigenvalues for T s are the same as those of |{+∞} except for a sign for s = 1/2. Spin reversal relates states on opposite sides of the equator as Fig. 6 for an illustration.

Descendant towers for nonzero commensurate twist
Now we turn to commensurate twist (6.15) with φ / ∈ {0, π} away from the (anti)periodic case. It will be instructive to consider the two commensurate twists ±φ simultaneously. Like with the mirror-pair of descendant towers in Section 8.1.3 we need to consider two copies of descendant towers. This time, however, the two copies include states at twist φ as well as states at opposite twist −φ. Indeed, in Appendix C.3 we show that the transfer matrices T s (u, φ) with twist φ are related to those with opposite twist by flipping all spins, is an eigenstate of (the Hamiltonian, and more generally) the transfer matrices T s (u, −φ) whenever |{u m } M m=1 is so for T s (u, φ), cf. (C.17). (Observe that φ = 0, π are the only two values for which φ = −φ mod 2π.) We study the example N = 6, ∆ = 1 2 (so η = iπ/3) and commensurate twist φ = ± 2π 3 . We wish to understand the descendant tower for primitive states with M 0 = 1 down spin. These 'twisted magnons' are discussed in detail in Appendix B.2. Consider |{u 1 } with u 1 = arctan tan(π/18) √ 3 , which is an eigenvector for the Hamiltonian with twist φ = 2π 3 . Its quasimomentum, see Eq. (B.9), is Next consider the primitive state |{u 2 } with u 2 = −arctan tan(π/18) √

3
, where we decorate the ket by a prime to indicate that it is an eigenstate for the Hamiltonian with twist φ = − 2π 3 . 11 This second state has opposite quasimomentum These two primitive states have the same eigenvalues for the energy (cf. Appendix B.2) and, as ε nfm = (−1) 0 = 1 in (7.8), the same 'twisted' momentum p = p + φ/N = p − φ/N = π mod 2π. We note, however, that these two states are not degenerate for the transfer matrices. For example, the fundamental transfer matrix T 1/2 acts by 64 a 1 t 6 + a 2 t 4 + a 3 t 2 + a 4 + a 5 t −2 + a 6 t −4 + a 7 t −6 |{u 1 } , (8.22) in terms of t = e u as usual, while To construct the descendant towers we turn to the corresponding states beyond the equator. For |{u 1 } this is 6 j=1 σ x j |{u 1 } ∝ |{u 1 , −∞, α fm 1 , whose eigenvalue for the Q operator The state beyond the equator corresponding to |{u 2 } is |{u 2 , +∞, α fm 2 , with eigenvalue for the Q operatorQ(t, φ) given by Here the string centres of the two FM strings are The descendant towers can now readily be constructed using the 'free fermion'-like property as in Section 8.1.1. The resulting tower structure is depicted in Fig. 7, where in general states come in pairs that have the same (possibly up to a sign, as always) eigenvalues for all transfer matrices e ∓isφ T s (u, ±φ), where the sign in ±φ is determined by the Hamiltonian for which the state is an eigenvector.
A key difference with the (anti)periodic case is that, in order to construct the descendant tower by considering the state beyond the equator corresponding to the primitive state, we are led to consider the systems with twists φ and φ = −φ simultaneously. The states within either tower in Fig. 7, e.g. |{u 1 } and |{u 1 , −∞} , have degenerate (possibly up to a sign) eigenvalues for all e ∓isφ T s (u, ±φ), with sign of the twist determined by the states. On the other hand, between the two different towers, the eigenstates are only degenerate for the momentum (possibly up to a shift by π) and energy: as we saw for the two primitive states this does not extend to the higher charges generated by the transfer matrices. This is the second important difference with the (anti)periodic case from Section 8.1.3, where all eigenstates in Fig. 6 have degenerate (up to a possible sign) eigenvalues for T s . To conclude this discussion we remark that the descendant-tower structure for twists ±φ / ∈ {0, π} can get more complicated for roots of unity with 2 > 3.

Full spectrum at root of unity: an example
Although we have not been able to find an algorithmic description of the structure of the spectrum for given values of the spin-chain parameters N, ∆, φ, our numerical investigations of numerous examples suggest that the full spectrum can be understood case by case in terms of the descendant towers that we have described in this section. Let us demonstrate this for a system with N = 10, ∆ = 1 2 (η = iπ 3 ) and φ = 0. The resulting description of the full spectrum is given in Table 1.
We stress that, due to the vanishing Wronskian relation (6.14) and the 'free fermion'like behaviour of FM strings, each primitive state is at the top of a descendant tower, with descendants that are obtained by adding FM strings or pairs of roots +∞, −∞. By (2.18) the states come in three categories, classified by the number of Bethe roots located at infinity: i) n +∞ = n −∞ (= 0 for primitive states), ii) n +∞ > n −∞ = 0, iii) n −∞ > n +∞ = 0.
Start with |↑ · · · ↑ = |{∅} . From the zeroes of the Q function (7.23) we find (cf. Footnote 9 on p. 36) that the corresponding state beyond the equator can be written as |↓ · · · ↓ ∝ |{α fm 1 , α fm 2 , 2 × ±∞} , with FM strings of length 2 = 3. We proceed as in Section 8.1.2. Draw the candidate descendant tower using the 'free-fermion' property. Note that (2.15) and (2.17) exclude all potential descendant states in this tower with one pair +∞, −∞ of infinite roots: all four infinite roots must be added at the same time. The primitive state at M = 0 thus gives rise to a descendant tower with two descendants at M = 3 and at M = 7, and one descendant at each of M = 4, 6, 10.
All ten states |{u 1 } at M = 1 must be primitive. Their Bethe root are finite in view of (2.15) and (2.17). The corresponding states beyond the equator have room for two FM strings and one pair of infinite roots (n +∞ = n −∞ = 1, allowed at M = 9). The descendant tower is as in Fig. 5, except that there is just a single pair of infinite roots. All intermediate descendants occur: the pair of infinite roots at M = 3, 6, 9 is allowed by (2.15) and (2.17). Each M = 1 vector thus sits at the top of a descendant tower containing one descendant at each of M = 3, 7 and two descendants at M = 4, 6.
The 45 states at M = 2 are also primitive, and must have finite Bethe roots. The corresponding states beyond the equator allow for two FM strings. By the free-fermion property each descendant tower contains four states.
Next we turn to M = 3, where 120 − 12 = 108 primitive states remain. Here it starts to become a little more complicated to predict the structure by hand, as infinite Bethe roots (n ±∞ = 1) may occur. Due to parity invariance there must be equally many states in the classes (ii) and (iii). It is possible to find out how many such states occur. Indeed, suppose that n +∞ = 1, n −∞ = 0 (the opposite case is treated analogously). Then two Bethe roots remain to be determined. Their Bethe equations effectively acquire a twist due to the presence of the infinite root, cf. the discussion preceding (2.12). One can solve these Bethe equations, and delete all solutions that themselves contain infinite roots. Alternatively one can directly use the eigenvalues of the truncated two-parameter transfer matrixT(x, y, φ) to compute the eigenvalues of the Q operatorQ numerically. In either case, the corresponding states beyond the equator need four more Bethe roots. For the primitive states with n +∞ = n −∞ = 0 these extra roots come from n +∞ = n −∞ = 2, while for the primitive states with n ±∞ = 1 and n ∓∞ = 0 these descendants have n ±∞ = 2, n ∓∞ = 0 plus one FM string.
One proceeds analogously for the 210 − 12 = 198 primitive states left at M = 4. The corresponding states beyond the equator need two more Bethe roots, which come from adding a pair +∞, −∞ for the primitive states in class (i). The situation is a bit more tricky for the primitive states in the other classes. For class (ii) they are of the form |{u 1 , u 2 , +∞, +∞} , but the corresponding states beyond the equator are only allowed to have one infinite root. In this case one of the infinite roots is removed to make place for an FM string: the descendants are |{u 1 , u 2 , α fm , +∞} .

Deformations of FM strings
At the start of Section 7.3 we showed that FM strings occur precisely when the system is at root of unity with commensurate twist (6.15). Let us further substantiate this by studying what happens to the Bethe roots constituting an FM string under small changes to the twist φ or the anisotropy parameter η. (Observe that such deformations immediately break the relation (6.15), as 2 fluctuates wildly when the value of η varies through i Q ⊂ i R.) Turning on a small twist φ (Appendix E.1) leads to smooth changes of the Bethe roots associated to FM strings, and the same is true for the eigenvalues of the Q operator. This is closely related to the string hypothesis in the thermodynamic limit, as explained in Section 10.  ) and φ = 0, organised in terms of primitive and descendant states with n +∞ = n −∞ , n +∞ > n −∞ = 0, n −∞ > n +∞ = 0.
Bethe roots associated to the FM strings (Appendix E.2): unlike before, FM strings cannot be continuously reconstructed when η is deformed, and the eigenvalues of the Q operator change drastically. This is expected from the representation theory of U q (sl 2 ), which is very sensitive to the precise root of unity. However, as one would expect on physical grounds, physical quantities such as the eigenvalues of the transfer matrices T s do change continuously as η is varied. Further investigations are required in order to fully understand this behaviour.

Connection to the work of Fabricius and McCoy
In [31] Fabricius and McCoy derived an equation for the string centre of FM strings at zero twist, see (1.11) therein, by linearising the Bethe equations (2.7) around an arbitrary root of unity. In other words, these string centres are obtained by continuity from their values in an infinitesimal neighbourhood of the root of unity. Let us explain why this is not in contradiction with the results from Section 8.4 and Appendix E.2.
By 'Bethe roots' we mean the zeroes of the eigenvalue of Q operator constructed from the two-parameter transfer matrix. We thus determine the FM string centres at commensurate twist (6.15) from the zeroes of the Q operator. Such zeroes associated with FM strings are among the zeroes of the Drinfeld polynomial (parameters of the evaluation representation). Note that the Drinfeld polynomial Y (v) in (1.42) of Ref. [32], cf. (4.2) therein and [39], is proportional to Q s P s in our notation, as can be seen by comparing (1.42) in [32] with our (7.15).
Next, our roots associated with FM strings do not satisfy (1.11) in [31]. This can most easily been seen when N is a multiple of 2 and 1 is even (i.e. L a multiple of N and r even in the notation of [31]), so that (1.11) from [31] reduces to (1.10) therein. In that case the denominator of Y (v) is one, and the zeroes of Y (v) coincide (up to a simple change of notation) with the zeroes of (7.23), which we identify with the roots constituting the FM string. This would coincide with the solutions of [31] if one would ignore the second factor in parenthesis in (1.10).
The solutions to (1.11) in [31] can be used to construct the eigenvectors of H. Indeed, Fabricius and McCoy proposed a creation operator for FM strings, see (1.38) and (1.41)- (1.42) in [32]. In examples we find that, taking the spectral parameter therein to by a solution to (1.11) in [31], the result is an eigenvector of H and T 1/2 , but not of the two-parameter transfer matrix or Q operator. (Note that we cannot take the spectral parameter to be the FM string centres as found by our methods, cf. Section 8.5, since these are zeroes of Y (v) which appears as a denominator in (1.38) from [32].) Instead, the FM string creation operator from [32] yields a linear combination of the eigenvectors of the Q operator. In Appendix E.2 we illustrate this with an example.
In the next section we propose creation and annihilation operators for FM strings that do yield eigenvectors of the Q operator.

Conjectures for FM creation and annihilation operators
Vernier et al. [36] proposed to use semicyclic representations in order to construct degenerate states with different magnetisation. These representations have also been used to construct quasilocal charges [77]. Building on these ideas we present conjectures for explicit constructions of the FM string creation and annihilation operators that commute with the XXZ transfer matrix while changing the magnetisation in steps of 2 .

Case q 2 = +1
We start with roots of unity obeying ε = q 2 = +1. Pick the semicyclic representation for the auxiliary space: V sc a is 2 dimensional, has spin s ∈ C and furthermore depends on a parameter β ∈ C, see Appendix A.2. Replace X a from (4.1), (6.3) by X sc a = X a + β x y/(x − y) |0 l − 1| a , which still obeys the commutation relation W a X sc a = q X sc a W a . Substituting this in (4.7) modifies the matrix entry S + a = L 21 a / sinh η to S +, sc while keeping the other three matrix elements in (4.7) unchanged. In this way we obtain a Lax operator L sc aj , and the usual construction (3.5)-(3.6) gives the transfer matrix T sc s (u, φ, β) = tr a L sc aN (u, β) · · · L sc a1 (u, β) E a (φ) (9.2) depending on s ∈ C and β ∈ C. Here the twist E a (φ) acts on V sc a by the usual expression (B.1). The matrix elements of (9.2) change the magnetisation (N − 2M )/2 by −n 2 and are proportional to β n for positive n ∈ Z >0 . For example, in a chain of length N = 2 , the expansion of the transfer matrix contains a term tr a [(S +, sc a ) N ] N j=1 σ − j that changes the magnetisation by − 2 , creating 2 magnon excitations.
The matrices T sc s (u, φ) do not commute with each other at different values of u, but do commute with T 1/2 (u, φ) when the twist is commensurate. To see this, note that the construction of the semicyclic Lax operator guarantees that the following version of the RLL relation (3.3) with (six-vertex) R-matrix (3.4) holds true for any β: This is an identity of operators on V j ⊗ V k ⊗ V sc a . Let us reinterpret it by thinking of the spin-1/2 space labelled by k as an auxiliary space, V k V b . The symmetry property R jk = R kj of (3.4) allows us to view the R-matrix as a Lax matrix L bj , whereas L sc ak takes the role of an R-matrix R sc ab . Reversing the two sides of (9.3) and changing v → u − v we arrive at an RLL relation on V sc Here a corresponds to the semicyclic auxiliary space, b to a spin-1/2 auxiliary space, and k to a spin-1/2 physical space. The train argument implies that the twisted semicyclic transfer matrix (9.2) commutes with the twisted fundamental transfer matrix (3.5)-(3.6), provided R sc ab commutes with the tensor product of the twists matrices. This requires the twist to be commensurate, e i 2 φ = 1, cf. (6.15).
By fusion (9.2) further commutes with T s (u, φ) for any 2 s ∈ Z >0 : Since the twisted semicyclic transfer matrix T sc s (u, φ, β) changes the magnetisation of an eigenstate of the Q operator in steps of 2 , it mixes states that are degenerate for T s (u, φ) within each descendant tower.
One can use T sc s (u, φ, β) to construct eigenstates of the Q operator itself. Because the part of T sc s (u, φ, β) of first order in β changes the magnetisation of eigenstates of the Q operator by − 2 , we make Conjecture 1. When ε = q 2 = +1 the linearisation in β of (9.2) at s = ( 2 − 1)/2, is the creation operator for the FM string {α fm } = {u, u + iπ 2 , · · · , u + iπ 2 −1 2 }. The spectral parameter can be taken to be any Bethe root from the FM string, e.g. u as in (9.6).
Note that at β = 0 and 2s = 2 − 1 the semicyclic Lax operator coincides with the Lax operator whose auxiliary space is the 2 -dimensional highest-weight representation. By (9.2) the operator (9.6) can thus be expressed more explicitly as where we write e n,n a = |n n | a for the matrix units on V a . The construction of the FM-string creation operator (9.6) is just like that of the generating function of the quasilocal Y charges proposed in Ref. [77], except that the transfer matrices are evaluated at different values of s. Therefore B fm (u), like those Y charges, commutes with T s (u) with 2s ∈ Z >0 but not withT s (u) when 2s ∈ C \ Z >0 .
Next we turn to the annihilation operator. Let us denote the parameter of the cyclic representation by γ. If we repeat the same construction with the transposed Lax matrix L aj (v y , u x ) from (4.9), this time replacing the entryŜ − a := L 12 a (v y , u x )/ sinh η bŷ we obtain a semicyclic transfer matrix that changes the magnetisation by positive multiples of 2 . Like before we arrive at Conjecture 2. When ε = q 2 = +1 the linearisation in γ at s = ( 2 − 1)/2, tr a L aN (u) · · · L a,j+1 (u) e 2 −1,0 a σ + j L a,j−1 (u) · · · L a1 (u) E a (φ) , (9.13) annihilates the FM string with Bethe roots {u, u + iπ 2 , · · · , u + iπ 1 −1 2 }. We have again checked this in many examples including the one above. We find that the consistency condition C fm (α fm n ) B fm (α n ) |{u m } m ∪{α fm n } = |{u m } m ∪{α fm n } = B fm (α n ) C fm (α fm n ) |{u m } m ∪ {α fm n } holds so long as |{u m } m ∪ {α fm n } and |{u m } m ∪ {α fm n } are eigenstates of the Q operator, even though we do not know the commutation relations between B fm (u) and C fm (v) in general. We postpone such 'off-shell' relations to future investigations.
Recall that for ε = −1 the commensurate twist φ depends on the system size N through the condition e i 2 φ = (−1) N from (6.15). Therefore on V sc Combining (9.14) and (9.17) the train argument, illustrated in Figs. 8 and 9, yields , provided it exists, and taking the trace over the 2 -dimensional auxiliary space we see that the semicyclic transfer matrix (9.16) commutes with the fundamental transfer matrix in the sense that With the aid of the fusion relation we obtain the commutation relations In particular, the semicyclic transfer matrix commutes with T s for integer s .
Remark. Another way to get this result is to evaluate the monodromy matrix M sc a with an 2 2 -dimensional semicyclic auxiliary space V a when ε = −1. One can verify that tr a M sc if an extra twist E = 2 1 −1 n=0 |n + 2 mod 2 2 n| a is included. In either way we are led to Conjecture 3. For any ε = q 2 = ±1 the linearisation in β of (9.16) at s = ( 2 − 1)/2, creates an FM string with Bethe roots {u , u + iπ 2 , · · · , u + iπ 2 −1 2 }. The spectral parameter in (9.21) is related to the FM string by u = u if ε = +1 and u = u − iπ 2 2 when ε = −1. This conjecture has been successfully tested on all examples in Section 8, as well as for various cases when 2 is even.
Let us focus on ε = −1 again. Taking the derivative of (9.20) with respect to the parameter β we obtain the (anti)commutation relations The anti commutation of B fm (u) with T 1/2 (v) for ε = −1 might be surprising, but can be understood as follows. When ε = −1 an FM string adds π to the momentum without affecting the energy or other charges. This requires B fm (u) to anticommute with the translation operator, and to commute with all higher conserved charges (logarithmic derivatives of T 1/2 ). This is guaranteed by the anticommutation. The annihilation operator for FM strings can be defined like in (9.13). Generalise the semicyclic transfer matrix (9.12) to arbitrary ε ∈ {+1, −1} by including staggered parameters ±γ like in (9.16). Then we propose Conjecture 4. When ε = q 2 = ±1 the operator ε j tr a L aN (u) · · · L a,j+1 (u) e 2 −1,0 a σ + j L a,j−1 · · · L a1 (u) E a (φ) (9.23) annihilates an FM string with Bethe roots {u , u + iπ 2 , · · · , u +iπ 2 −1 2 }. The spectral parameter is again related to the FM string by u = u if ε = 1 and u = u − iπ 2 2 for ε = −1.

FM strings and Z charges
In Section 7.2 we have seen that all states within a descendant tower have the same eigenvalues (up to a possible minus sign) for T s . This may lead one to wonder if those states can be distinguished at all using charges that are (quasi)local in the thermodynamic limit. The answer to this question is positive: the (exponentially many) degeneracies can be lifted by taking into account the quasilocal Z charges [42][43][44][45] that can be constructed at root of unity as logarithmic derivatives of the truncated transfer matrixT s from Section 6. The quasilocality of the Z charges in the thermodynamic limit was demonstrated in Refs. [42][43][44][45]. The Z charges were used to study out-of-equilibrium phenomena such as quantum quenches in Ref. [41]. Here we focus on their role in distinguishing states that are degenerate for T s . In Section 6 we already studied the key ingredient for the construction of the Z charges, viz. the truncated two-parameter transfer matrixT(x, y, φ) =T s (u, φ) at root of unity. Similar to Ref. [41] we define the generating function of the Z charges as Like in (3.8) the Z charges are the coefficients of the expansion around u = η 2 , The Z charges are able to tell apart each member of the descendant tower. Let us illustrate this for the example N = 6, ∆ = 1 2 (η = iπ 3 ) and φ = 0 from Section 2.3.Consider the descendant tower formed by where the string centres for the two FM strings are The generating function (10.1) acts on the descendant tower by The eigenvalues of Z(u) are different for each eigenstate. Now consider the thermodynamic limit N → ∞. There is a subtlety in the presence of FM strings. For instance, when φ = 0 and η = iπ 1 / 2 with 1 odd, (6.15) implies that for odd N there are no FM strings in the spectrum, while for even N there are (exponentially many) states associated with FM strings. The spectrum of systems at finite size is therefore sensitive to the parity of N , and it is not a priori clear whether the thermodynamic limit is well defined. However, numerics suggests that at odd N there are states in the spectrum that differ by terms that vanish as N → ∞ to give rise to the same asymptotic degeneracies as obtained when taking the limit via even N . Thus the result in the thermodynamic limit should be independent of the way in which the limit N → ∞ is taken.
The role of the Z charges in separating the degeneracies in the presence of FM strings has important implications for the thermodynamic limit: Z charges have to be included when constructing the generalised Gibbs ensemble (GGE) for the XXZ spin chain at root of unity. Z charges are also important to obtain non-vanishing spin Drude weight with XXZ model at root of unity [42,43]. This is further supported by the discussion in the remainder of this section.
10.2 TBA, string-charge duality, and a conjecture for string centres of FM strings One of the cornerstones in the study of the thermodynamic properties of quantum integrable models is the thermodynamic Bethe ansatz (TBA), which has been used extensively to study out-of-equilibrium problems such as transport, quantum quenches and generalised hydrodynamics.
In order to solve the TBA equations for the XXZ spin chain at root of unity one relies on the string hypothesis [78], which stipulates the types of bound states that the model permits. This is a powerful tool that has yielded numerous results for thermodynamic properties of the XXZ model. A general review can be found in [79]. In the following we will use the notation of Ref. [46].
In [46], Ilievski et al. made the remarkable observation that the root densities of the Bethe strings (bound states) are in one-to-one correspondence with the generating functions of (quasi)local charges (3.8) of the model. This string-charge duality is of great use and convenient for the study of expectation values of (quasi)local charges, especially for quantum quenches [12]. The string-charge duality furthermore provides the root densities for the longtime non-equilibrium steady states [46,80].
When we are at root of unity we have to take into account the generating function of the Z charges to obtain the root densities of the 'last two string types' in the list of Takahashi [78,79]. This result has been derived in [41,46]. 12 Moreover, the total root densities of the 'last two string types' are the same from the TBA calculation (see Eq. (9.66) of [79]), revealing a deep connection between the two.
The 'last two string types' are related to FM strings too. In fact, two Bethe strings, one of each of the last two string types, that have coinciding real parts of their string centres (cf. Footnote 9 on p. 36) together form an FM string. For example, consider ∆ = 1 2 (η = iπ 3 ), for which the string hypothesis says that the last two string types are (2, +) and (1, −). Call a bound state with n magnons whose Bethe roots have the same real part an 'n-string'. Then '(2, +)' denotes a 2-string with even parity (centred at the real axis, i.e. a complex conjugate pair), and '(1, −)' a 1-string with odd parity (centred at Im u m = π 2 ≡ − π 2 ). For η = iπ 3 the imaginary parts of the 2-string with even parity are ± π 6 while the imaginary part of the 1-string with odd parity is π 2 . On the other hand, according to the examples in Section 8, FM strings for η = iπ 3 can be expressed as But this can be viewed as two strings, one of type (2, +) and one of type (1, −), with coinciding real part of the string centre. We expect that for any system size with commensurate twist FM strings can be decomposed in terms of the last two string types of the string hypothesis in this way, even at finite system size. Indeed, for any principal root of unity η = iπ 2 the last two string types are ( 2 − 1, +) and (1, −) [78]. We have investigated FM strings for all 2 ≤ 6 at various N and found that they can all be viewed as consisting of strings of the last two string types with equal real parts of the string centres. It even appears to hold for non-principal root of unity, see Appendix F for examples.
By studying the string centres of FM strings as the value of φ approaches the commensurate twist (6.15) as in Appendix E.1 we are led to Conjecture 5. For any finite system size N and root of unity η = iπ 1 2 with commensurate twist, any Fabricius-McCoy string has string centre with imaginary part This conjecture has been confirmed for all examples in Section 8.
In the thermodynamic limit, the conjecture (10.8) is compatible with the known TBA results [79]. Due to the exponentially large degeneracies of the descendant tower, FM strings in the thermodynamic limit will contribute to the thermodynamic quantities. Through conjecture (10.8) the density of FM strings has been included already in the root densities of the 'last two string types'. It implies that we do not need to include any new functional relation to the TBA formalism in [79].
Combining the conjecture with the knowledge that the generating function of the Z charges is crucial to obtain the root densities of 'last two string types', we conclude that the Z charges are directly related to the FM strings.
Remark. There is a key difference between the thermodynamics of the 'descendant states' in our article and the sl 2 -descendants in spin-1/2 Heisenberg XXX spin chain. For the latter the sl 2 -descendant states are associated with Bethe roots at infinity (vanishing quasimomenta), and do not enter in the TBA calculations. However, for the XXZ spin chain at root of unity the 'descendant states' in our terminology (see Section 7.1) do enter the TBA calculation through the density of FM strings, which is in principle determined by the density of the last two string types.

FM strings and spin Drude weight
One of the most important physical consequences of the quasilocal Z charges is the nonvanishing high-temperature spin Drude weight 13 of the XXZ model at root of unity, due to the non-commutativity between the spin flip operator j σ x j and the Z charges [42,43,45]. It can be considered as a manifestation of the exponentially many degeneracies in the thermodynamic limit.
In Appendix E.2 we demonstrate that perturbing the anisotropy parameter η away from root of unity can change the structure of the Bethe roots dramatically, especially for states with FM strings. On the other hand, the spin Drude weight [45] is also known to change significantly under such perturbations. This hints at a connection between the existence of FM strings and the fractal structure of the spin Drude weight.
Another hint for the intimate relation between FM strings and non-vanishing spin Drude weight at root of unity comes from the domain-wall quench, i.e. the time evolution of an initial state |↑ · · · ↑↑↓↓ · · · ↓ , for the XXZ model at root of unity. Here ballistic spin transport (nonvanishing spin Drude weight) can be treated analytically using generalised hydrodynamics [83]. The right half of the system, i.e. the fully polarised state |↓↓ · · · ↓ , has Q function given in Section 7.5 for finite size. In the thermodynamic limit, according to the TBA, this fully polarised state consists of a filled 'Fermi sea' with Bethe roots of the last two string types [79,83], cf. Section 10.2. In this case, each pair of Bethe strings with one string from each of the last two string types has the same real parts of the string centres. From the definition of the density of Bethe strings it follows that the densities of the last two string types must coincide. The conjecture (10.8) implies that the Bethe roots of |↓↓ · · · ↓ consists solely of FM strings. For the domain-wall quench the ballistic spin transport from the right half of the system is solely carried out by the quasilocal Z charges [83]. Notice that, even though the FM strings do not contribute to the dynamics, the 'FM strings' of the right half of the system are not true excitations from the perspective of the whole system. Thus they do contribute to the dynamics, as combinations of Bethe strings of the last two string types, resulting in the domain-wall melting phenomenon.

Conclusion
We have studied the full spectrum of the transfer matrices associated to the quantum spin-1 2 Heisenberg XXZ chain, focussing on root of unity with arbitrary twist. To this end we constructed the Baxter's Q operator, and the P operator, from the factorisation of a twoparameter transfer matrix (4.25). The eigenvalues of the Q operator, i.e. Q functions, are polynomials whose zeroes encode the physical solutions of the Bethe equations (2.7).
As a by-product of our construction we rederived the matrix TQ relation (5.10) and transfer-matrix fusion relations (5.12) from a decomposition of the two-parameter transfer matrix, providing a simplification of the conventional approach. At root of unity we further derived truncated transfer-matrix fusion and Wronskian relations from the two-parameter transfer matrix with auxiliary space truncated to a finite-dimensional space. We also proved an interpolation-type formula conjectured in Refs. [32,34,41].
Equipped with these algebraic tools we obtained analytic results about the spectrum at root of unity (Section 7). These results enabled us to demonstrate the presence of descendant towers for the XXZ model at root of unity with commensurate twist, for with explicit examples for the various scenarios that occur (Section 8). We elucidated the exponential growth of the degeneracies at root of unity. Since we can construct the Q operator explicitly at root of unity, via a trace over a finite-dimensional auxiliary space, we obtained analytic results for rather large system size compared to previous works. From the quantisation condition (7.15) we found that FM strings associated with the descendant states behave like free fermions within each descendant tower (Section 8.1.1). We have found new semicyclic transfer matrices that satisfy unconventional RLL relations, see e.g. (9.18), from which we conjectured an explicit expression for the creation and annihilation operators of FM strings (Section 9).
Even though our main results and discussions concentrate on systems with finitely many sites, we moreover compared our results with recent works on the thermodynamic limit (Section 10). We explained the relation between the truncated two-parameter transfer matrices and the quasilocal Z charges, which are of crucial importance in many applications such as quantum quenches and spin transport at root of unity. Inspired by the string-charge duality we found a connection between the last two string typesà la Takahashi, the Z charge and the FM strings in the XXZ model at root of unity. This led us to a conjecture about the imaginary part of the string centres of FM strings based on the string hypothesis.
Outlook. Several interesting aspects remain to be explored. First of all we hope to apply the construction of the two-parameter transfer matrix to other quantum integrable models, such as the XXZ spin chain with other boundary conditions as well as their higher-spin generalisations. It is not known how the construction should be generalised to integrable models with higher-rank symmetry, where there exist several Q operators.
We would like to understand more properties of the FM strings in the thermodynamic limit, e.g. whether quasilocality of the Z charges has a deep relation to FM strings (cf. Section 10). In addition, the semiclassical limit of the domain-wall quench in the gapless regime yields a similar ballistic spin transport behaviour [84,85] to the quantum counterpart [83]. However, it is not clear what the semiclassical limit of the Z charges would be since the classical Landau-Lifshitz field theory does not have an underlying quantum-group structure. It would be desirable to clarify the relation between the mechanisms in both quantum and classical regime when a qualitative classical-quantum correspondence of spin transport [84] applies.
The structure of descendant towers from Section 8 implies the existence of a hidden Onsager algebra at least for a part of the physical Hilbert space [86]. At the free-fermion point ∆ = 0, and at η = iπ 3 for the integrable spin-1 XXZ chain, the Onsager algebra plays a crucial role. The properties of the Onsager algebra may allow one to obtain results like in Ref. [36]. The relation between the Onsager generators and the semicyclic transfer matrices from Section 9 needs to be elucidated. In particular we would like to discover algebraic structures like those in [36] in the extensively studied, yet not completely understood, XXZ spin chain.
The full spectrum of the XXZ model at root of unity has potential physical applications. For instance, equipped with our results, it would be interesting to calculate the partition function of the six-vertex model at root of unity analytically using algebraic geometry, following recent works on the isotropic case [87,88]. Furthermore, other physically relevant quantities such as the density matrix and overlaps with a generic state can be extracted from the full spectrum at finite system size, which could shed light on the studies of quantum entanglement and non-equilibrium dynamics in exact solvable models.
A Quantum sl 2 The quantum group U q (sl 2 ) is the unital associative algebra with generators S + , S − along with K (which is invertible) subject to the commutation relations We take the coproduct to be S ± → S ± ⊗ K −1 + K ⊗ S ± and K → K ⊗ K. There is a counit and antipode, see e.g. Eqs. (1.2)-(1.4) in Ref. [47]. A good (but hard to find) introduction is [89].
In this appendix we summarise the representations of U q (sl 2 ) that we will use. For each case it is easy to check that the commutation relations (A.1) hold and that K = exp(η S z ). The spin of an irrep is defined by the eigenvalue [s] q [s + 1] q of the quantum Casimir operator which generates the centre of U q (sl 2 ).

A.1 Global representation
The physical Hilbert space (C 2 ) ⊗N of the spin chain comes with two 'global' representations. When N = 1 the representation is given by S ± = σ ± and K = q σ z /2 . For N > 2 repeated application of the coproduct gives the (reducible) representation from (2.5). By reversing the factors (taking the opposite coproduct) we obtain another (reducible) representation: All of these operators can be obtained from the entries (3.9) of the monodromy matrix is the ('braid') limits u → ±∞. In particular, the B operator from the QISM is closely related to the above spin-lowering generators. To see this we write the Lax operator (3.2) with spin-1 2 auxiliary space (see Footnote 4 on p. 11) in the form In the ('braid') limits u → ±∞ the diagonal entries behave as Therefore, using the definition of monodromy matrix (3.9), we have One similarly recovers the spin-raising operators from C, and K, K −1 from either of A, D.

A.2 Auxiliary representations
In this article we use various choices for the auxiliary space, which is always a representation of U q (sl 2 ). We summarise the key ingredient in this appendix. First of all we use the finite-dimensional unitary spin-s representation of U q (sl 2 ). Denote the (orthonormal) basis of C 2s+1 by |n for n = 0, 1, · · · , 2s. Then the generators are given by When s = 1/2 this gives the case N = 1 of (A.3) when we identify |1 = |↑ and |0 = |↓ . Note that we are thus led to a decreasing ordering of the basis: we prefer conventions such that S ± |n ∝ |n ± 1 and the basis is labelled by non-negative integers; then the matrices match the standard choice provided we order the basis as |2 s , |2 s − 1 , · · · , |0 . Next we use the complex spin-s highest weight representation of U q (sl 2 ). It is defined on an infinite-dimensional Hilbert space with orthonormal basis n| indexed by n ∈ Z ≥0 . The generators are given by [2s − n] q |n + 1 n|, This is related to the transfer matrices T hw s (4.19). In Figs. 2 and 3 we write down explicit matrices; let us stress that -albeit perhaps somewhat unusual in an infinite-dimensional context -as above we take the basis to be ordered decreasingly: · · · , |2 , |1 , |0 . When s ∈ 1 2 Z ≥0 the infinite-dimensional highest-weight representation contains a finite-dimensional submodule. Indeed, the entry of S + for n = 2s vanishes (cf. Fig. 2), so the subspace labelled by 0 ≤ n ≤ 2s is preserved by all of (A.9). Thus we can truncate to a representation of dimension 2s + 1 with generators [2s − n] q |n + 1 n|, This representation is equivalent to (A.8) by a gauge transformation (conjugation). More importantly, when η = iπ 1 / 2 there exists another truncation yielding an 2dimensional representation, illustrated in Fig. 3. This is sometimes referred to as a nilpotent representation, due to the fact that (S ± ) 2 = 0 in this case. The generators act on the subspace with basis |n for 0 ≤ n ≤ 2 − 1 by [2s − n] q |n + 1 n|, [n + 1] q |n n + 1|, There is one more truncated 2 -dimensional representation that we will use at root of unity: the semicyclic representation. It is similar to the truncated highest-weight representation (A.11) with an additional entry in S + : [2s − n] q |n + 1 n|, [n + 1] q |n n + 1|,

B.1 Twist operator
We define the twist operator E a (φ) for the auxiliary space. Each of the U q (sl 2 ) representations on the auxiliary space from Appendix A.2 is expressed in terms of an orthonormal basis {|d − 1 , · · · , |1 , |0 } with d the dimension of the representation. Here d = 2 s + 1 for the unitary spin-s representation with s ∈ 1 2 Z, d = ∞ for the highest-weight representation with s ∈ C, and d = 2 for the truncation at root of unity. We consider diagonal twist operator E a (φ) given by In view of our ordering of the basis this yields the twist from (3.9) for s = 1/2 (d = 2). In particular, the complex spin-s representation yields monodromy matrix M hw resulting in the transfer matrix T hw s (u, φ) = tr s M hw s (u, φ). When |q| ≤ 1 the diagonal matrix elements of M hw s can be bounded by A N |q N e iφ | n for some constant A N , and so the trace is convergent if |e iφ | < |q| N . At this point we do not know if T hw s can be analytically continued outside this disc of convergence.
The truncation at root of unity η = iπ 1 / 2 likewise has and transfer matrixT s (u, φ) = tr sMs (u, φ). In this case the trace is well defined for any value of of the twist φ.

B.2 Twisted momentum and magnons
Let us compute the first few charges (3.8) generated by the (six-vertex) transfer matrix with auxiliary space V a ∼ = C 2 of spin 1 2 and twist E(φ) = diag(e iφ , 1) as in (3.9)-(3.10). In the periodic case (φ = 0) we get the cyclic (right) translation operator G(0), T(η/2, 0) = sinh N η G(0), G(0) = P 12···N = P 12 P 23 · · · P N −1,N , (B.4) so the charge I (1) = −i log G is the usual momentum operator. In the quasiperiodic case the twist deforms the translation operator to In analogy with the periodic case write its eigenvalue as e ip ; we will get back to the 'twisted momentum' p soon. The next charge is the XXZ Hamiltonian (2.1) up to some constants: where the prime denotes the derivative with respect to the first argument. The twist spoils homogeneity in the traditional sense: for φ = 0 (B.6) does not commute with (B.4). However, has G (φ)-eigenvalue e i p where p = 2πk/N , 0 ≤ k ≤ N − 1, is quantised. When φ is real G (φ) is unitary so the twisted magnons (B.7) are linearly independent and form a basis for the M = 1 sector, so they are also eigenvectors of (B.5) and (B.6). In terms of p := (2πk − φ)/N the (twisted) momentum is p = p + φ/N , the dispersion is ε = cos( p ) − ∆, and the right-hand side of (B.7) looks just like an ordinary magnon, j e −i pj σ − j |Ω . (The sign in the exponential is because we work with the right translation operator.) Let us finally make contact with the algebraic Bethe ansatz: for M = 1 (3.11) gives where the quasimomentum solves the Bethe equation e iN p = e −iφ . The (twisted) momentum and energy given above match the result (2.8)-(2.9) obtained from (3.12) using (B.5) and (B.6). The difference between quasimomentum and (twisted) momentum for M = 1 can be avoided by taking G := e −iφ G to be the twisted translation.

C Bethe roots C.1 Numerical recipe for finding Bethe roots
Here we review a numerical recipe, called McCoy's method, for solving the functional TQ relation. It appears to have been published first in [91]. We follow the description of Haldane [92]. A similar method was also described by Baxter [20]. The idea is that rather than solve the (coupled, nonlinear) Bethe equations one can obtain the Bethe roots by solving a few sets of linear equations. This is done by exploiting the known form of the eigenvalues of the (fundamental) transfer matrix and the Q operator, whose zeroes are the Bethe roots (see Section 5.5). The recipe goes as follows: i. Construct the transfer matrix T 1/2 (u) at some u ∈ C (almost any value will do), and numerically diagonalise it; this is much more efficient than diagonalisation when u is kept free. One obtains 2 N eigenstates that span the physical Hilbert space. From the Bethe ansatz we know that the eigenvectors are independent of the spectral parameter, so these will be eigenvectors for T 1/2 (u) for any u.
ii. The eigenvalues, however, depend on u. Pick one of the eigenvectors. Its eigenvalue is found by acting with T 1/2 (u). This may again be done numerically by writing the eigenvalue as a Laurent polynomial in t = e u of order N , with zeroes τ n that can be fixed by acting with T 1/2 (u n ) for N distinct values u n ∈ C.
iii. The corresponding Bethe roots are the zeroes of the Q operator, found by solving the functional TQ relation (5.19), i.e.
Here T 0 (u) = sinh N (u) and the eigenvalues are of the form where M is the number of down spins of the eigenvector under consideration. The zeroes t m can once more be found numerically by taking t = e u equal to the zeroes τ n of T 1/2 and solving the linear problem.
The zeroes give the Bethe roots u m = log t m . One needs to be careful to interpret the result correctly in the presence of Bethe roots at infinity: u m = ±∞ corresponds to t m ∈ {0, ∞} so the corresponding factor in (C.3) collapses to t ±1 , yielding (5.21): The numerical recipe works very well for the XXZ model away from root of unity, as well as for the XXX model (∆ = ±1). However, one cannot find all the Bethe roots for the XXZ spin chain at root of unity, precisely due to the existence of degenerate eigenstates of the transfer matrix T 1/2 . In that case, our construction for the Q operator still works and gives the correct results, sometimes even analytically.

C.2 Relation between Bethe roots for anisotropies ∆ and −∆
In the gapless regime (−1 < ∆ < 1) there is a simple relation between Bethe roots of all the physical solutions at anisotropy ∆ and those at anisotropy −∆, even though the corresponding eigenstates are different since the B-operators in the algebraic Bethe ansatz differ. We will denote the parameters of the second spin chain by primes: ∆ = −∆ and η = arccosh(∆) ∈ i R, η = arccosh(∆ ) = i π − η.
In terms of these parameters (C.6) reads This precisely of the form (C.6) with η = iπ − η and twist φ chosen such that e −iφ = (−1) N e −iφ . This shows that for each solution {u m } M m=1 at anisotropy ∆ there is a corresponding solution {u m } M m=1 at ∆ = −∆ provided the twist is modified to The two eigenstates are related by the unitary gauge transformation (Removing the prefactor in the expression on the right yields a simpler transformation that also does the job and is its own inverse.) It is easy to check that this transformation changes the sign of ∆ in the Hamiltonian (2.1): Moreover, the eigenstates are related by

C.3 Relation between eigenstates with opposite twist
Recall from Section 5.5 that an M -particle Bethe state |{u m } M m=1 for the XXZ model obeys with eigenvalue Q(u) and P (u) of the form 'beyond the equator' with opposite twist.
Consider the transfer matrices T s (u) with s ∈ 1 2 Z ≥0 . Under global spin inversion these operators simply behave as To see this note that for the unitary spin-s representation (A.8) the Lax operator (3.2) is invariant under total spin reversal, which acts by conjugation by σ x j in the physical space and by the antidiagonal matrix U := 2s n=0 |2s − n n| in the auxiliary space. Spin reversal in the physical space is thus equivalent to spin reversal in the auxiliary space. This property is inherited by the monodromy matrix. In the periodic case it follows that spin flip in the physical space does not affect the transfer matrix, as the trace is invariant under reordering the basis. In the twisted case we only have to correct for the twist (B.1), which is reversed: U E(φ) U −1 = e 2siφ E(−φ). This proves (C.15). Now consider the TQ equation (5.10) with φ inverted to −φ. By conjugating both sides with the global spin-flip operator N j=1 σ x j , using (C.15) and multiplying both sides by e iφ we see thatQ(u, φ) := N j=1 σ x j Q(u, −φ) N j=1 σ x j precisely obeys the TP equation (5.11). Moreover, comparing the eigenvalues in (C.14) shows that the eigenvalues ofQ(u, φ) on Mparticle Bethe vectors are trigonometric polynomials of degree N − M , just as for the P operator. It follows that the eigenvalues ofQ(u, φ) are proportional to those of the P operator; in particular they have the same zeroes: (C. 16) SinceQ(u, φ) and Q(u, −φ) have the same characteristic polynomial this shows that the Mparticle eigenvalues of the P operator are the same as the (N − M )-particle eigenvalues of Q(u, −φ). But we know that the latter can be interpreted as the Bethe roots. Therefore the zeroes of the P operator can be interpreted as the Bethe roots of the spin-reversed Bethe vector beyond the equator with opposite twist. Finally notice that the Bethe vectors (3.11) are constructed using the B-operator, which is independent of the twist, see (3.9). This implies that the result of reversing all spins on an off-shell Bethe vector (for the Hamiltonian with original twist φ) is In this appendix we give another proof of (6.13), i.e. Consider 2s ∈ Z ≥0 with s < 2 2 −1. Proceeding as in Section 5.1 we obtain a decomposition like in (5.5):T s (u, φ) = T s (u, φ) + e i(2s+1)φ T * s (u, φ).

(D.2)
Here T * s denotes the result of restricting the trace over the auxiliary space to the subspace V * a spanned by |n with 2s + 1 ≤ n ≤ 2 − 1.
The restricted Lax matrix L * sj (u) coincides with the transpose (in both auxiliary and physical space) L −s−1,j (u) T of the Lax operator L −s−1,j (u), whose auxiliary space has dimension 2 − 2s − 1. Therefore T * s (u, φ) = T −s−1 (u, φ), and we get Now consider the transfer matrixT * s obtained likeT s truncating the trace to the basis ofṼ a translated (raised) by 2s + 1, i.e. to |n + 2s + 1 with 0 ≤ n ≤ 2 . Again, since 2s + 1 < 2 and S + | 2 = 0, the trace can be split into two transfer matrices. The first one coincides with T * s , while the second one is obtained as T s in the basis shifted by 2 , i.e. |n + 2 with 0 ≤ n ≤ 2s + 1, and therefore picks up a factor ε N compared to T s times the twist contribution. Proceeding as before, one can verify that the second transfer matrix in the decomposition ofT * s andT −s−1 derive from similar Lax matrices and are both equal to e i( 2 −2s−1)φ ε N T s . Thus one has Subtracting e i( 2 −2s−1)φ times (D.3) from (D.4) we obtain (6.13).

E Deforming FM strings E.1 Tuning a small twist
We study the Bethe roots for states that include FM strings at φ = 0 numerically in the presence of a small twist. The results illustrate that FM strings are formed by combining the last two string types from the string hypothesis. Consider the case N = 6 and ∆ = 1 2 , like the example in Section 2.3. At zero twist φ = 0 the primitive state |↑↑↑↑↑↑ has, by the results of Section 7.5, corresponding state beyond the equator with Bethe roots (2.21). The descendant tower further includes the states with energy eigenvalues E = 0. Now we turn on a very small twist and study what happens to the corresponding states by computing the eigenstates with M = 3 and eigenvalues for the transfer matrices T s that are very close to those of (E.1). The resulting Bethe roots are collected in Tables 2-3, where we have expressed the analytic initial values (2.21) numerically to facilitate the comparison. We observe that the Bethe roots are string deviations of two (2, +) strings, namely u 1 , u 2 and v 1 , v 2 , with string deviations that decrease as φ approaches zero. In the limit φ → 0 an FM string forms when two strings, here of type (2, +) and (1, −) respectively, have coinciding real part of the string centres. This motivates our conjecture in Section 10.2 about the relation between the last two string types of the string hypothesis and the string centres of FM strings, even for systems with finite system sizes.

E.2 Tuning the anisotropy
Next we study the behaviour of Bethe roots for states including FM strings at φ = 0 when η is slightly deformed. We give two examples, both with ∆ ≈ 1/2. First we revisit the example with N = 6 from Section 2.3 and Appendix E.1. When ∆ = 1/2 the states (E.1) with Bethe roots (2.21), i.e. numerical values given by the first rows of Tables 2-3, are degenerate for the fundamental transfer matrix, with eigenvalue (recall that t = e u ) where we give the numerical values for later convenience. In particular their energy is E = 0. Now we vary the anisotropy a little, taking η = i ( π 3 ± ) for small > 0, to move away from root of unity. The numerically obtained values of the Bethe roots of the two states at M = 3 with energy E ≈ 0 are given in Table 4, including the values that one would obtain as limits at = 0. (These values obey Eq. (1.11) from [31].) Note that the limiting values at root of unity are rather different from the Bethe roots (2.21) obtained from the zeroes of the Q operator. Here we use the notation  Table 4: The numerical Bethe roots of the two states at M = 3 with energy close to zero as η changes away from iπ/3, including the values 0, ∓iπ/3 and iπ/2, ∓iπ/6 (indicated with a question mark) that one would extrapolate as φ → 0.
eigenvectors of H with E = 0. However, these vectors are not eigenvectors of the twoparameter transfer matrix or the Q operator at ∆ = 1/2. We find that nontrivial linear combinations of them (namely: a multiple of their sum and difference) equal |{u m } 3 m=1 and |{v m } 3 m=1 , which are eigenvectors of Q operator, respectively. Despite the jump in the Bethe roots {u m } m and {v m } m as → 0 we stress that physical quantities are continuous. For instance, fundamental transfer matrix has eigenvalue where t = e u . The behaviour of the Bethe roots under a small change of η is given in Table 5. We see that the Bethe roots change drastically, jumping to form a 2-string and two 1-strings with odd parity. The same is true for the eigenvalues of the Q operator. On the other hand, one can check that the eigenvalues of T s change smoothly. For instance, the eigenvalue of T 1/2 for the corresponding state at η = i(π/3 + 10 −5 ) is very close to (E.7):

F Examples of last two string types of TBA at non-principal root of unity
We present several examples on the last two string types of the TBA at root of unity using Takahashi's notation [78,79] at non-principal root of unity to illustrate our conjecture (10.8) for the string centre of FM strings. We start with η = 2iπ 3 , which is a non-principal root of unity. The allowed string types are (1, −), (2, +), (1, +). Here we have underlined the last two string types, which are of the form α 1 = α − iπ 3 , α 2 = α + iπ 3 , α ∈ R, α 3 = α , α ∈ R. (F.1) According to the conjecture (10.8) the FM string with η = 2iπ 3 should be expressed as Clearly, (F.1) and (F.2) describe the same FM string when the real parts of the string centres coincide. This results are confirmed by all the examples in Section 8 with finite-size calculations.