Ground state solutions of inhomogeneous Bethe equations
S. Belliard A. Faribault


























shapes





    ϵłλδ̣αβ̱γσc̅^
            ⟨#|1⟨#1||#⟩1|#1⟩⟨#|1⟩⟨#1⟩#1<#1|#1|#1>
            
            
            
            Ground state solutions of inhomogeneous Bethe equations
            
            
            
            
            
            
            S. Belliard1 2  †
            A. Faribault3 *
            
            
            
             1 Sorbonne Université, CNRS, Laboratoire de Physique Théorique et Hautes Energies,  LPTHE,  F-75005, Paris, France
            
             2  Institut de Physique Théorique, DSM, CEA, URA2306 CNRS Saclay, F-91191, Gif-sur-Yvette, France
            
             3 Université de Lorraine, CNRS, LPCT, F-54000 Nancy, France.
            
            
            † samuel.belliard@gmail.com
            
             alexandre.faribault@univ-lorraine.fr
            
            
            December 30, 2023
            
            
            
            § ABSTRACT
             
            The distribution of Bethe roots, solution of the inhomogeneous Bethe equations, which characterize the ground state of the periodic XXX Heisenberg spin-1/2 chain is investigated. Numerical calculations shows that, for this state, the new inhomogeneous term does not contribute to the Baxter T-Q equation in the thermodynamic limit.  Different families of Bethe roots are identified and their large N behaviour are conjectured and validated.
            
            
            
            
            
            
            
            
            
            
            
            
            § INTRODUCTION
            
            
            The Bethe-Hulthén ansatz approach  <cit.> has shown its interest by allowing on to investigate the thermodynamic limit of spin chains <cit.> and predict exact physical quantities.
            The spectrum of the Hamiltonians of integrable spin chains follows from the one of the transfer matrix t(λ), the generating function of conserved quantities of the 
            model and can be characterized by the Baxter T-Q equation <cit.> which, in the case of the periodic XXX Heisenberg 
            spin-1/2 chain with N spins is given by
             
            t(λ)q̂(λ)= (λ + i/2)^N q̂(λ - i)+ (λ - i/2)^N q̂(λ + i),  
             
            where q̂(λ) is the Q operator which commutes with the transfer matrix and whose eigenvalues have the polynomial form
            
            q(λ)=∏_k=1^M(λ-λ_k),
              
            with M≤N/2.
            
            The recent development of the Bethe ansatz for the models without U(1) symmetry, where the usual Bethe-Hulthén ansatz approach fails, leads to the discovery of a new family of Baxter T-Q equation <cit.>, called  inhomogeneous Baxter T-Q equation or  modified Baxter T-Q equation. 
            This new family differs from the usual one in two ways:  it involves a new term and the order of the associated Q polynomial is fixed. 
            The implementation of the Bethe ansatz changes and the techniques to obtain this equation and the associated states lies at the intersection of many different methods: Off diagonal Bethe ansatz <cit.>, Modified algebraic Bethe ansatz <cit.>, Separation of variables <cit.>. Moreover, it was observed in <cit.> that this  inhomogeneous Baxter T-Q equation can also characterize 
            finite models with U(1) symmetry and, in the case of the periodic XXX Heisenberg 
            spin-1/2 chain with N spins, is given by [A more general characterization depending on an arbitrary parameter can be constructed <cit.>, here for numerical convenience we fix it to its simplest form which also corresponds to the diagonal limit of the T-Q equation for the twisted XXX Heisenberg 
            spin-1/2 chain in the parametrization of <cit.>.]
            t(λ) Q̂(λ) =- (λ+i/2)^N  Q̂(λ-i)-(λ-i/2)^N  Q̂(λ+i) +4(λ^2+1/4)^N,
            
            where Q̂(λ) is the inhomogeneous Q operator whose eigenvalues have the polynomial form
            
            Q(λ)=∏_i=1^N(λ-u_i).
            
             It is worth mentioning that the existence of such an operatorial T-Q equation is a conjecture contrary to its functional form.
            
            The understanding of the thermodynamic limit for models with a spectrum characterized by such modified Baxter T-Q equation 
            is an important challenge. For models with a boundary term, it was observed numerically in <cit.> that the contribution to the ground state energy for the inhomogeneous term decreases as N^-1, where N is the length of the chain. 
            
            
            In this work, we investigate the evolution of the ground state's Bethe roots as a function of the chain length N and classify the types of roots which appear in the solution of the inhomogeneous Bethe equations. In order to perform the thermodynamic limit a key ingredient is, indeed, knowledge of the structure of the Bethe equation's solution. For example, the periodic XXX Heisenberg spin-1/2 chain, the antiferromagnetic ground state belongs to the class of solution where q(λ) have only real roots <cit.> (while excited states can also have complex roots that behave as strings <cit.>). 
            Considering that the Hamiltonian of this model can be characterized by both the homogeneous Baxter T-Q equation (<ref>) and the inhomogeneous Baxter T-Q equation (<ref>), one can retrieve, using the well-know homogeneous case, the Bethe roots of the inhomogeneous Baxter Q polynomial, without directly solving the inhomogeneous Bethe equations. Here, we will exclusively consider the distribution of Bethe roots which characterizes the antiferromagnetic ground state.
            
            The paper is organized as follows, in the section <ref> we recall the Bethe equations, their logarithmic form, and their relation to the homogeneous Baxter T-Q equation (<ref>). The ground state's real Bethe roots are then calculated up to N=300. In section <ref> we reconstruct the Q polynomial from the inhomogeneous Baxter T-Q equation  (<ref>) and the homogeneous solution giving the transfer matrix in order to numerically find the solution of the  inhomogeneous Bethe equations characterizing the antiferromagnetic ground state of the model. Three different families of Bethe roots are then identified and the large N behavior is discussed. Concluding remarks are given in section  <ref>.  
            
            
            
            
            § HOMOGENEOUS GROUND STATE SOLUTION
             
            
            
            The roots of the Q baxter polynomial q(λ) (<ref>),  {λ_1 …λ_M} satisfy the Bethe equations
            (λ_k + i/2)^N ∏_j k^M (λ_k-λ_j - i) =
            (λ_k - i/2)^N ∏_j k^M (λ_k-λ_j + i),
             
            with k ∈{1,...,M} and M ={0,1,...,N/2} (we always consider N even).
            These equations follow from the Baxter TQ equation (<ref>) evaluated at the zeros of q(λ). 
            The knowledge of the Bethe roots  {λ_1 …λ_M} allows from equation  (<ref>)  a direct numerical computation of the eigenvalue t(λ) at any point λ. This formulation would lead to numerical difficulties in evaluating t(λ) in the immediate vicinity of a Bethe root, as we should divide by q(λ), a problem which is explicitly avoided by the numerical approach used in the following.
            Instead of trying to directly solve the Bethe equations in the form (<ref>), it is much more stable for numerical computation to use the logarithmic form
            
            N ln(1 - 2 i λ_j/1+ 2 i λ_j)= i (N+M+1+J_j) π
            + ∑_k j^M [ln(1 -  i (λ_j -λ_k)/1 +  i (λ_j -λ_k))],
            
            
            
            
            
             where J_j is an integer which specifies the various possible logarithmic branches. 
            For a size, N=4n, such that both N and M=N/2 are even, the ground state solution is found for the values J_j = 2 (j - 1-N) with j=1, 2 …,N/2. 
            
            A simple implementation of the Newton-Raphson method 
            is then sufficient for equation (<ref>) to rapidly converge to the desired ground state solution starting from almost any initial approximation for the position of the Bethe roots. 
            One can exploit the symmetry of the ground state root distribution which come in pairs of ±λ_i, in order to reduce the size of the system by half. Indeed, one can solve for only the positive roots: λ_1 …λ_N/4 by rewriting the N/4 first equations in terms of these variables having replaced the remaining ones by {λ_N/4 +1…λ_N/2}={-λ_1… -λ_N/4}. 
            The resulting solutions are presented in Figure <ref> and uniquely define the corresponding q(λ) polynomial, i.e. the N/2 order polynomial with zeros placed at each of these roots.
            
            
            
            § INHOMOGENEOUS  T-Q GROUND STATE SOLUTION
             
            
            The roots of the Q baxter polynomial Q(λ) (<ref>),  {u_1 … u_N} satisfy the Bethe equations
            (u_k - i/2)^N ∏_j k^N (u_k-u_j + i)-(u_k + i/2)^N ∏_j k^N (u_k-u_j - i) =4 (u_k^2 + 1/4)^N,  
             with k={1,..,N}.
            These equations follow from the Baxter TQ equation (<ref>) evaluated at the zeros of Q(λ). 
            
            Having computed the homogeneous q polynomial through its N/2 zeros in the previous section and therefore the corresponding t(λ) polynomial which stays the same in inhomogeneous case and the homogeneous one, we can retrieve the inhomogeneous Baxter Q operator from the inhomogeneous T-Q equation (<ref>) and find it roots without having to explicitly solve the inhomogeneous Bethe equations (<ref>). This greatly simplifies the calculation since, due to the inhomogeneous term, they cannot be simply brought into logarithmic form so that their explicit numerical solving is a much more arduous task than in the homogeneous case. 
            
            
            
             §.§ Reduction of the inhomogeneous T-Q equation
            
            
            Considering that for any λ, equation (<ref>) couples exclusively the polynomial Q at three points λ,λ± i, it becomes simple to define a system of sparse linear equations. Indeed, through a Lagrange polynomials representation, one can define any polynomial P(λ) of order N, in terms of the values P(z_i) it takes at a set of N+1 arbitrary grid points {z_1 … z_N+1}. 
             In the case at hand (<ref>), since the leading coefficient is know to be 1, only N points are necessary to define exactly
             
             
             Q(λ) = ℓ(λ) [1+ ∑_n=1^N Q(z_n)/ℓ'(z_n)1/λ-z_n]     with   ℓ(λ) ≡∏_n=1^N (λ- z_n).
             
            
            
            
            Considering that the T-Q equation only couples  Q(λ) and Q(λ± i), one can choose to use a grid  of points {z_1 … z_N}, which are separated by i allowing us to fully define Q(λ) in terms of the N values Q(z_i). These Q(z_i) will simply be the solution to a set of (mostly) tridiagonal linear equations.
            Indeed, on a grid such that z_n+1 = z_n + i , at every point except for z_1 and z_N,  equation (<ref>) reads
            
            
            t(z_n) Q(z_n)  + (z_n+i/2)^N Q(z_n-1)+(z_n-i/2)^N Q(z_n+1) = 4 (z_n^2+1/4)^N,
            
            leading to a core set of N-2 linear equations coupling only three of the variables Q(z_n). 
            The two extreme points z_1 and z_N still lead to full equations with non-zero coefficients in front of each variables {Q(z_1), …, Q(z_N+1)} due to the presence of Q(z_1 - i) or Q(z_N+1+i) which, using equation (<ref>), could be expressed in terms of a full sum over the values {Q(z_1), …, Q(z_N)}. 
            
            When considering ground state solutions for even N/2 (N being an integer multiple of 4) we know, from the previous section, that the t(λ) polynomial can be explicitely deduced from the homogeneous Baxter T-Q equation with a q(λ) polynomial whose zeros are all real and come in symmetric pairs: ±λ_i. This fact allows us to conclude that
            
            
            t(z) = t(-z) = t^*(z^*)         Q(z)= Q(-z) = [Q(z^*)]^*.
            
            It therefore leads to further simplifications to choose the (i-separated) grid points z_n to lie on the imaginary axis and to be symmetric around the real axis, namely: z_n= (- N+1/2 + n) i for n = 1,2 … N. By doing so, one insures that Q(z_i) ∈ℝ and that grid points come in pairs, i.e. there exists z_n = -z_n' allowing us to reduce the system of linear equations to only N/2 equations using variables Q(z_i) which will take real values.
            
            Additionally, it is easily shown that, for any N integer multiple of 4,  a pair of Bethe roots is always found at ±i/2. Indeed looking at the specific equations obtained from (<ref>) at ±i/2, one has:
            
            
            t(i/2) Q(i/2)  + Q(-i/2) = 0,     
            t(-i/2) Q(-i/2)  +  Q(i/2) = 0.
            
            Considering as well that  t(-i/2)=t(i/2) and  Q(i/2)  =Q(-i/2), we simply find: 
            
            [ t(i/2)+1 ] Q(i/2)  = 0,
            
            and therefore to  Q(i/2) =  Q(-i/2) =0 since equation (<ref>) gives t(i/2) = q(-i/2)/q(i/2) = 1.
            
            Having shown that Q(z) has a pair of roots at ±i/2 and that our choice of a purely imaginary symmetric set of grid points allows us to reduce the expression for Q(z) by half, one can define:
            
            
            Q(z) = (z-i/2) (z+i/2) Q̃(z),
            
            where Q̃(z) is a symmetric polynomial (Q̃(z)=Q̃(-z)) of order N-2 with leading coefficient 1 which obeys the T-Q̃ equation:
            
            
            t(λ) Q̃(λ)    =   - (λ+i/2)^N-1(λ-3i/2)Q̃(λ-i)-(λ-i/2)^N-1(λ+3i/2)Q̃(λ+i) 
                      +4 (λ^2+1/4)^N-1.
            
            The polynomial Q̃(z) can be written, in a Lagrange form, only in terms of the grid points in the upper complex plane as:
            
            Q̃(z) = ℓ̃(z) ∑_n=1^N/2-1Q̃( 2n+1/2i )/ℓ̃'( 2n+1/2i )[1/z-( 2n+1/2)i - 1/z+( 2n+1/2)i], 
            
            with  ℓ̃(z) = ℓ(z)/(z-i/2)(z+i/2)=∏_n=1^N/2-1(z-( 2n+1/2)i)(z+( 2n+1/2)i),
             
            where we used the antisymmetry relation: ℓ̃'( 2n+1/2i ) = - ℓ̃'(-  2n+1/2i ).
            
            In the end, the different values {Q̃(3i/2), Q̃(5i/2), Q̃(7i/2) …Q̃((N-1)i/2)} needed to reconstruct the polynomial Q̃(z) (and therefore Q(z)) are found as the solution to the linear system of equations:
            
                   t(3 i/2) Q̃(3 i/2) + 3 Q̃(5 i/2) + 2^N + 1 = 0
            
            
                  t(k i + i/2)Q̃(k i + i/2)
              + (k + 1)^N - 1 (k - 1) Q̃( (k-1) i + i/2 ) 
                      + (k + 2) (k)^N - 1Q̃((k+1) i + i/2)  =- 4 (k(k + 1))^N - 1   ∀  k= 2,3 … N/2-2 
            
            
                   t(0) 2^N Q̃(0) - 6 Q̃(i)- 2^4-N =0.
            
            The last equation, completing the system, is simply obtained from
            (<ref>) evaluated at λ=0. It is the only one of these equations which involves every Q̃(z_i) when explicitly written in terms of  {Q̃(3i/2) …Q̃((N-1)i/2)} using the Lagrange basis representation (<ref>) of Q̃(0) and Q̃(i).
            While the resulting linear system requires more than machine precision to be solved numerically even for modest chain lengths N, its sparsity allows one to go to higher precision calculations in what remains a reasonable computation time, even on a standard personal computer. 
            
            Having then found the set of {Q̃(3i/2), Q̃(5i/2), Q̃(7i/2) …Q̃((N-1)i/2)}, one can use any standard solver, at high precision, to find the roots of the resulting Q̃ polynomial, by simply numerically finding the points at which equation (<ref>) equals 0.
            
            
            
            
             §.§ Numerical solution and roots families
            
            
            The general structure of the ground state solution obtained numerically, see figure 2, is threefold. First, it systematically contains a series of N/2 real roots. Secondly, there is a string of purely imaginary roots and finally a series of roots with both non-zero real and imaginary parts which create four symmetric arcs which emerge from both ends of the imaginary string. 
            
            
            
            The general result is best visualised when plotting each of these substructures individually, and therefore, the real-valued roots are first plotted in figure 3.
            
            
            Beyond the realization that the number of such real roots is indeed systematically the same as the number of roots in the homogeneous case, namely N/2, one easily infers from the plot, that, in the large N limit, this set of real roots will correspond exactly to the homogeneous solution. It will therefore become dense in the thermodynamic limit. 
            
            In figure 4 we exclusively plot the purely imaginary roots as well as the number of such roots in figure 5 for system sizes N which are integer multiples of four. 
            
            
            
            
            As can be seen from large enough system sizes, these imaginary roots tend to arrange into of perfect symmetric string of i-separated Bethe roots. The only visible deviations from this ideal string structure appear at the very top and bottom ends of the string.
            
            
            
            
            Within the studied range, the number N_I of such imaginary roots stays bounded by N/8≤ N_I ≤N/8+9/2. The upper bound corresponds exactly to N_I for N=44, N=76, N=108 and N=140, while the lower bound does so for N=256 and N=288. While this does not allow one to conjecture that these bound hold true for arbitrary system sizes, it seems safe to conjecture that, as system size grows, the length of this imaginary string will also keep growing. Indeed, although when going from N to N+4, it is possible to see N_I go down (by 2) therefore giving a shorter string, comparing any N to the corresponding N+8 it is found that N_I is never reduced.
            
            The last subset of roots for which both real and imaginary parts take non-zero values form the arc-like structures and are explicitly shown in figure 6 for the upper right complex plane.
            
            
            
             Higher values of N always lead to arcs whose both end points are further from the real axis, so that, in the figure, they naturally appear as ordered by chain length N. Contrarily to the real roots which densify as N becomes large, the roots on these arcs actually get further apart from one another as the system size gets larger. Consequently, one cannot expect, in the thermodynamic limit N→∞ to be in a position to define a continuous root density to properly describe these structures. 
            
            Nonetheless, these numerical results support strongly the conjecture that, as N→∞, these complex roots will all be such that their norm is also going to infinity.
            
            
            
            The conjectures made from these finite size results can now be explicitly verified to give the correct solution to the inhomogeneous T-Q equation in the thermodynamic limit N→∞. Indeed, we have this threefold structure for which:
            
            
            1- N/2 real Bethe roots u^r_a will tend to the solution of the homogeneous problem.
            
            
            2- N_I imaginary Bethe roots u^i_a that will form an infinitely long string of i-separated roots, with N_I going to infinity.
            
            
            3- The remaining N/2 - N_I  roots u^c_a form discrete points on an arc infinitely far from the dense distribution of real roots |u^c_a|>>|u^r_b|.
            
            
            Under these assumptions, the Q polynomial can be decomposed in three parts
            
            Q(λ) = q^R(λ)q^I(λ)q^C(λ)
            
            with
            
            q^R(λ)=∏_a=1^N/2 (λ -u^r_a),  
            q^I(λ)=∏_a=1^N_I (λ - u^i_a),  
            q^C(λ)=∏_a=1^N/2-N_I (λ - u^c_a)
            
            which allows one to rewrite the inhomogeneous T-Q equation (<ref>) as
            
            
            t(λ) q^R(λ)    =    - (λ+i/2)^N q^R(λ-i)q^I(λ-i)q^C(λ-i)/q^I(λ)q^C(λ)
                  -(λ-i/2)^N q^R(λ+i) q^I(λ+i)q^C(λ+i)/q^I(λ)q^C(λ)
            + 4 (λ^2+1/4)^N/q^I(λ)q^C(λ).
            
            The roots belonging to the complex arcs are such that | u^c_b|→∞ when N→∞ and therefore, in the same limit we have
            q^C(λ± i)/q^C(λ)→1.
            
            
            For arbitrary values of λ such that |λ| << |u^c_a| for every a = 1,2 … N/2-N_I, the inhomogeneous term becomes negligible compared to the three other ones and therefore, in this whole small |λ| region: 
            (λ^2+1/4)^N/q^I(λ)q^C(λ)→ 0.
            
            
            
            On the other hand, the contribution from the string of imaginary roots can be simplified when N_I goes to infinity. Indeed for the pure i-separated string, we have 
            λ^I_b=-iN_I/2+(b-1)  i
            
            and thus
            
            q^I(λ-i)=∏_b=1^N_I (λ-i - λ^I_b)=∏_b=1^N_I (λ+iN_I/2-(b-2)  i)=q^I(λ)λ-iN_I/2+i/λ+iN_I/2+i
            
            q^I(λ+i)=∏_b=1^N_I (λ+i - λ^I_b)=∏_b=1^N_I (λ+iN_I/2-b  i)=q^I(λ)λ+iN_I/2-i/λ-iN_I/2-i
            It follows that, for a string length N_I→∞, one has:
            q^I(λ± i)/q^I(λ)→-1.
            
            
            
            Combining those two limits of q^I and q^C for N→∞ we finally fall back (for small |λ|), on the homogeneous Baxter T-Q equation: 
            
            t(λ) q^R(λ) =  (λ+i/2)^N q^R(λ-i)
            +(λ-i/2)^N q^R(λ+i).
            
            
            Since a polynomial equation which is verified for a whole sector (small |λ|) is, by extension, verified everywhere, the remaining roots u^r_a are then indeed the solution of the homogeneous T-Q equation. The large N limit conjectures made above therefore provide, in the thermodynamic limit, the solution to the inhomogeneous problem. Therefore, for the ground state, the inhomogeneous term does not contribute to the Baxter T-Q equation in the thermodynamic limit.
            
            
            
            § CONCLUSION
            
            
            In this paper we investigated the distribution of the Bethe roots which are solution to the inhomogeneous Bethe equation and characterise the antiferromagnetic ground state of the Heisenberg XXX spin-1/2 chain. We observe that for this state, under some simple conjectures about the behavior of the Bethe roots as function of N, the inhomogeneous
            Baxter T-Q equation reduces to the homogeneous one. This  results from two facts: firstly a large set of complex roots go to infinity as N→∞ while secondly an infinite discrete purely imaginary string of roots changes the sign in front of the two first terms of the inhomogeneous
             Baxter T-Q equation. 
             
            
            As already mentioned, the inhomogeneous Baxter T-Q equation (<ref>) is not the unique parametrization (see footnote p 2). In the general case we have an arbitrary parameter α and the  inhomogeneous Baxter T-Q equation reads
              
            t(λ) Q(λ) =α(λ+i/2)^N Q(λ-i)+α^-1(λ-i/2)^N Q(λ+i) +(2-α-α^-1)(λ^2+1/4)^N.
            
            In this case, it is much more complex to solve numerically for large N due to the additional parameter and the asymmetric form of the equation. However, one could still expect a similar threefold structure to emerge, now characterized by:    
            q^I(λ± i)/q^I(λ)→α^± 1.
            
            
            We can also extend this hypothesis to the case of the XXX spin chain with an arbitrary twist (as considered in <cit.>). In that case one also has both homogeneous and inhomogeneous TQ equations. Namely:
            
            t(λ) q(λ) =κ_1 (λ+i/2)^N q(λ-i)+κ_2 (λ-i/2)^N q(λ+i),
              
            where the two parameters κ_i satisfy κ_1+κ_2=κ+κ̃, κ_1κ_2=κκ̃-κ^+κ^-,
            and
            
            t(λ) Q(λ) =(κ̃-ρ) (λ+i/2)^N Q(λ-i)+( κ -ρ) (λ-i/2)^N Q(λ+i) +2 ρ(λ^2+1/4)^N,
              
            with ρ^2-ρ(κ+κ̃)+κ^+κ^-=0. In that case one could still expect a similar threefold structure to emerge, now characterised by:    
            q^I(λ- i)/q^I(λ)→κ_1/κ̃-ρ,   q^I(λ+ i)/q^I(λ)→κ_2/κ-ρ,
            where we also have κ_1/κ̃-ρ=κ-ρ/κ_2.
            
            
            
            It would be of strong interest to extend this numerical analysis to excited states of the periodic Heisenberg XXX spin-1/2 chain as well as for models where the U(1) symmetry is broken (as the Heisenberg XXX spin-1/2 chain with boundaries). In the later case, 
            the trick to use the homogeneous solutions would fail and one should therefore try to solve directly the inhomogeneous Baxter equation or use an explicit diagonalisation of the transfer matrix to retrieve the roots of the Baxter Q polynomial.
             
            Having identified the Bethe roots distribution for the inhomogeneous Baxter TQ equation, we should be in  a position to consider possible analytical results along the lines of Yang and Yang's works, and many others (see for recent development <cit.>). This question should be discussed elsewhere. 
             
             
            One of the important questions would be to see if it is possible or not to find a model (or a specific state) where the inhomogeneous term can contribute to the Baxter T-Q equation in the thermodynamic limit. Actually, the limited number of cases considered so far in the literature have all shown that it is not the case.  
             
            
            
            
              
            Acknowledgements:
            We thank V. Pasquier, B. Vallet, P. Baseilhac, N. Crampé, D. Serban, N. Slavnov, R. Pimenta, K.K. Kozlowski  for discussions and their interest. S. B. would like to thank the LPCT of Université de Lorraine for its hospitality and the LMPT of Tours University for the invitation during which a part of this work was done. A. F. would also like to thank the LMPT of Tours University for the invitation during which parts of this research was carried out. 
            
            
            
            
            
              
            Funding information:
            S.B. was supported by a public grant as part of the Investissement d'avenir project reference ANR-11-LABX-0056-LMH, LabEx LMH. 
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
             
            	
            99Bet31 H. A. Bethe,  Zur Theorie der Metalle. i. Eigenwerte und Eigenfunktionen der linearen Atomkette, Zeit. für Phys.  71, 205 (1931), 10.1007%2FBF01341708Hul38 L. Hulthén,  Uber das Austauschproblem eines Kristalles, Arkiv Mat. Astrom. Fysik  26 A, 1-106 (1938)
            
            
            CP62 J. Des Cloizeaux and J. J. Pearson,  Spin-wave spectrum of the antiferromagnetic linear chain, Phys. Rev.,  128, 2131-2135 (1962) 10.1103/PhysRev.128.2131CP63 R. B. Griffiths,  Magnetization Curve at Zero Temperature for the Antiferromagnetic Heisenberg Linear Chain, Phys. Rev.  133, 768-775 (1963) 10.1103/PhysRev.133.A768YY66a C. N. Yang and C. P. Yang,  One dimensional chain of anisotropic spin-spin interaction. I. Proof's of Bethe's hypothesis for ground state in finite system, Phys. Rev.  150,  321-327 (1966), 10.1103/PhysRev.150.321YY66b C. N. Yang and C. P. Yang,  One dimensional chain of anisotropic spin-spin interaction. II. Properties of the ground state energy per lattice site in finite system, Phys. Rev.  150, 327-339 (1966), 10.1103/PhysRev.150.327Gau71 M. Gaudin,  Thermodynamics of the Heisenberg-Ising Ring for Δ≥ 1^*, Phys. Rev. Lett.   26 1301 (1971),  10.1103/PhysRevLett.26.1301Tak71 M. Takahasi,  One dimensional Heisenberg model at finite temperature, Prog. Theor. Phys.  46 401-215 (1971), 10.1143/PTP.46.401FT81a L. D. Faddeev and L. A. Takhtadzhyan,  Spectrum and scattering of excitations in the one-dimensional isotropic Heisenberg model, J. Sov. Math  24, 241-267 (1984), 10.1007/BF01087245FT81b L. D. Faddeev and L. A. Takhtadzhyan,  What is the spin of a spin wave?, Phys. Let. A  85, 375-377 (1981), 10.1016/0375-9601(81)90335-2Bax72
            Baxter R.J. ,Exactly solved models in statistical mechanics, Academic 
            Press, (1982), 10.1142/9789814415255_0002 Faribault A. Faribault, H. Tschirhart,  Common framework and quadratic Bethe equations for rational Gaudin magnets in arbitrarily oriented magnetic fields, SciPost Phys.  3, 009 (2017), 10.21468/SciPostPhys.3.2.009gaudinbook
            M. Gaudin,  La fonction d'onde de
            Bethe, Paris: Masson, (1983)
            korepin
            V. E. Korepin, N. M.  Bogoliubov and A. G.  Izergin, 
            Quantum Inverse Scattering Method and Correlation Functions, Cambridge: Cambridge
            University Press, (1993), 10.1017/CBO9780511628832ODBA
            Y. Wang, W.-L. Yang, J Cao and K. Shi,  Off-Diagonal Bethe Ansatz for Exactly Solvable Models
            , Berlin Heidelberg: Springer-Verlag, (2015), 10.1007/978-3-662-46756-5ODBA2
            J. Cao, W.-L. Yang, K. Shi and Y. Wang,  Off-diagonal Bethe ansatz and exact solution of a topological spin ring, Phys. Rev. Lett.  111, (2013),  137201, 10.1103/PhysRevLett.111.137201ODBA3
            Y-Y. Lie, J. Cao, W.-L. Yang, K. Shi and Y. Wang,  Thermodynamic limit and surface energy of the XXZ spin chain with arbitrary boundary fields, Nuc. Phys. B  884, 17-27  (2014), 10.1016/j.nuclphysb.2014.04.010ODBA4
            F. Wen, T. Yang, Z. Yang,J. Cao, K. Hao and W.-L. Yang,  Thermodynamic limit and boundary energy of the su(3) spin chain with non-diagonal boundary fields, Nuc. Phys. B  915, 119-134  (2017), 10.1016/j.nuclphysb.2016.12.003NC13
            R. I. Nepomechie and C. Wang,  Boundary energy of the open XXX chain with a non-diagonal boundary term, J. Phys. A  47, 032001 (2014), 10.1088/1751-8113/47/3/032001sklyanintrends
            E. K. Sklyanin,  Separation of Variables: New Trends, Prog. Theor. Phys.  118, 35 (1995) 10.1143/PTPS.118.35Niccoliap
            G. Niccoli and V. Terras,  Antiperiodic XXZ Chains with Arbitrary Spins: Complete Eigenstate Construction by Functional  Equations in Separation of Variables, Lett. Math. Phys.  105, 989
            (2015), 10.1007/s11005-015-0759-9Niccoliff
            G. Niccoli,  Form factors and complete spectrum of XXX antiperiodic higher spin chains by quantum separation of variables, J. Math. Phys.  54, 053516 (2013),  10.1063/1.4807078Kozl K. K. Kozlowski,  On condensation properties of Bethe roots associated with the XXZ chain,
            Commun. Math. Phys.  357 1009 (2018), 0.1007/s00220-017-3066-8belliardnou1
            S. Belliard and R. A. Pimenta,  Slavnov and Gaudin-Korepin Formulas for Models without U(1) Symmetry: the Twisted XXX Chain, SIGMA   11 099 (2015), 10.3842/SIGMA.2015.099
