Ground state solutions of inhomogeneous Bethe equations

The distribution of Bethe roots, solution of the inhomogeneous Bethe equations, which characterize the ground state of the periodic XXX Heisenberg spin-$\frac{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 [1,2] has shown its interest by allowing on to investigate the thermodynamic limit of spin chains [2,[4][5][6][7][8][9] 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 [11] which, in the case of the periodic XXX Heisenberg spin− 1 2 chain with N spins is given by whereq(λ) is the Q operator which commutes with the transfer matrix and whose eigenvalues have the polynomial form 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 [16], 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 [15], Modified algebraic Bethe ansatz [12,24], Separation of variables [20][21][22]. Moreover, it was observed in [15] 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 1 whereQ(λ) is the inhomogeneous Q operator whose eigenvalues have the polynomial form 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 [15,[17][18][19] 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 [2][3][4][5] (while excited states can also have complex roots that behave as strings [1]). Considering that the Hamiltonian of this model can be characterized by both the homogeneous Baxter T-Q equation (1) and the inhomogeneous Baxter T-Q equation (3), 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 2 we recall the Bethe equations, their logarithmic form, and their relation to the homogeneous Baxter T-Q equation (1). The ground state's real Bethe roots are then calculated up to N = 300. In section 3.1 we reconstruct the Q polynomial from the inhomogeneous Baxter T-Q equation (3) 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 4.

Homogeneous ground state solution
The roots of the Q baxter polynomial q(λ) (2), {λ 1 . . . λ M } satisfy the Bethe equations 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 (5), it is much more stable for numerical computation to use the logarithmic form 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 (6) 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 The resulting solutions are presented in Figure 1 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(λ) (4), {u 1 . . . u N } satisfy the Bethe equations with k = {1, .., N }. These equations follow from the Baxter TQ equation (3) 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 (3) and find it roots without having to explicitly solve  (7). 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 (3) 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 (4), since the leading coefficient is know to be 1, only N points are necessary to define exactly 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 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 (8), 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 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 ) ∈R 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 (3) at ± i 2 , one has: Considering as well that t(−i/2) = t(i/2) and Q(i/2) = Q(−i/2), we simply find: and therefore to Q(i/2) = Q(−i/2) = 0 since equation (1) 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: whereQ(z) is a symmetric polynomial (Q(z) =Q(−z)) of order N − 2 with leading coefficient 1 which obeys the T −Q equation: The polynomialQ(z) can be written, in a Lagrange form, only in terms of the grid points in the upper complex plane as: where we used the antisymmetry relation:˜ 2n+1 2 i = −˜ − 2n+1 2 i . In the end, the different values {Q (3i/2) ,Q (5i/2) ,Q (7i/2) . . .Q ((N − 1)i/2)} needed to reconstruct the polynomialQ(z) (and therefore Q(z)) are found as the solution to the linear system of equations: t(3i/2)Q(3i/2) + 3Q(5i/2) + 2 N +1 = 0 The last equation, completing the system, is simply obtained from (14) evaluated at λ = 0.
It is the only one of these equations which involves everyQ (z i ) when explicitly written in terms of {Q (3i/2) . . .Q ((N − 1)i/2)} using the Lagrange basis representation (15) ofQ(0) andQ(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.

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.  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 with which allows one to rewrite the inhomogeneous T-Q equation (3) as The roots belonging to the complex arcs are such that |u c b | → ∞ when N → ∞ and therefore, in the same limit we have 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: 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 and thus It follows that, for a string length N I → ∞, one has: 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: 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 (3) 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 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: We can also extend this hypothesis to the case of the XXX spin chain with an arbitrary twist (as considered in [24]). In that case one also has both homogeneous and inhomogeneous TQ equations. Namely: where the two parameters κ i satisfy κ 1 + κ 2 = κ +κ, κ 1 κ 2 = κκ − κ + κ − , and with ρ 2 − ρ(κ +κ) + κ + κ − = 0. In that case one could still expect a similar threefold structure to emerge, now characterised by: 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 [23]). 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.