Thermodynamics of Inozemtsev's Elliptic Spin Chain

We study the thermodynamic behaviour of Inozemtsev's long-range elliptic spin chain using the Bethe ansatz equations describing the spectrum of the model in the infinite-length limit. We classify all solutions of these equations in that limit and argue which of these solutions determine the spectrum in the thermodynamic limit. Interestingly, some of the solutions are not selfconjugate, which puts the model in sharp contrast to one of the model's limiting cases, the Heisenberg xxx spin chain. Invoking the string hypothesis we derive the thermodynamic Bethe ansatz equations (TBA-equations) from which we determine the Helmholtz free energy in thermodynamic equilibrium and derive the associated Y-system. We corroborate our results by comparing numerical solutions of the TBA-equations to a direct computation of the free energy for the finite-length hamiltonian. In addition we confirm numerically the interesting conjecture put forward by Finkel and Gonz\'alez-L\'opez that the original and supersymmetric versions of Inozemtsev's elliptic spin chain are equivalent in the thermodynamic limit.


Introduction
The Bethe ansatz has been one of the most powerful tools in the field of integrability in the past eighty years. Its origin dates back to Bethe's solution of the Heisenberg model for the ferromagnetic interaction of electrons from 1931 [1]. Since then, analysis of numerous models other than spin chains benefited greatly from this ansatz, including the one-dimensional Bose gas [2], twodimensional lattice models such as the six-vertex model [3] and even N = 4 super Yang-Mills theory [4,5]. Moreover, many extensions of Bethe ansatz have been found, including the thermodynamic Bethe ansatz [6,7], nested Bethe ansatz [8] and asymptotic Bethe ansatz [9,10,11].
Heisenberg's spin-1/2 xxx spin chain is still ubiquitous in the research field centred around the Bethe ansatz. In an effort to generalize this spin chain, Inozemtsev proposed an elliptic spin chain characterized by the hamiltonian where L is the number of sites of the spin chain, J is the interaction parameter and ℘ L is the Weierstraß elliptic function with periods (L, iπ/κ) (for κ > 0) (see Appendix A) and σ is the usual vector of Pauli spin-1/2 operators [12]. Amazingly, this spin chain not only generalizes the Heisenberg xxx spin chain, which is recovered by taking κ → ∞, but actually interpolates smoothly between the (nearest-neighbour) xxx spin chain and the long-range Haldane-Shastry spin chain, obtained in the limit κ → 0. The Haldane-Shastry spin chain is solvable by exploiting its Yangian symmetry already present at finite length [13,14]. Therefore, investigating Inozemtsev's elliptic spin chain may shed light on the relation between these two methods for finding exact solutions. In particular, the integrability of both the Heisenberg xxx spin chain and the Haldane-Shastry spin chain suggest that Inozemtsev's elliptic spin chain might also be integrable. Although a definite proof remains absent to date, research into this question has culminated in a proposed set of L conserved quantities [15] and a description of eigenstates at finite and infinite L, which were found using an extended version of Bethe ansatz [16]. Another piece of evidence interestingly comes from the analysis of the level density of the spectrum of the spin chain, which agrees to great accuracy with some existing conjectures about chaotic versus integrable behaviour of quantum systems [17]. In fact, the spectrum of Inozemtsev's elliptic spin chain has been studied before. Dittrich and Inozemtsev probed the spectrum of Inozemtsev's infinite-length spin chain by classifying its twoparticle bound states [18]. Later, this spin chain was also used in a completely different context to calculate the first corrections to the dilatation operator in N = 4 super Yang-Mills theory. To this end asymptotic Bethe ansatz for Inozemtsev's spin chain was used to calculate corrections to the spectrum of the Heisenberg xxx spin chain as a truncated power series in κ [5], thereby providing some perturbative results on the spectrum of Inozemtsev's elliptic spin chain in the large volume limit.
Finally, the related supersymmetric su(1|1) version of Inozemtsev's elliptic spin chain was studied in [19] and shown to be integrable. It interpolates smoothly between the supersymmetrizations of the Heisenberg xxx spin chain (the xx spin chain at critical strength of the magnetic field) and of the Haldane-Shastry spin chain [20,21]. The thermodynamic limit of the su(1|1) elliptic spin chain was studied and shown to correctly reproduce the behaviour of the aforementioned models in the appropriate limits. In addition, the Heisenberg xxx and Haldane-Shastry spin chain turn out to be equivalent to their supersymmetrizations in the thermodynamic limit and it has been hypothesized that this equivalence also carries over to the elliptic spin chain.
In this work, we aim to gain additional information about the spectrum in the thermodynamic limit by invoking the string hypothesis [1], i.e. by assuming that the solutions of the Bethe ansatz equations in the infinite-length limit completely describe the thermodynamic behaviour of the model 1 . After characterizing all the solutions, usually called strings, we find integral equations describing the system in the thermodynamic limit. This method is quite standard for integrable models [7,23,24] and can be viewed as an extension of the method brought forward by Yang and Yang in [6]. 2 In Section 2 we will recall the relevant models and point out some of their important properties. In Section 3 we review the method of finding strings from a general perspective and apply it to the case of Inozemtsev's spin chain. In Section 4 we give arguments why not all the found strings can be used to parametrize the spectrum and present a set of strings that should describe the thermodynamics of Inozemtsev's elliptic spin chain. In Section 5 we apply the thermodynamic Bethe ansatz to these solutions to derive a set of integral equations that yield the free energy per site. In Section 6 we compare a numerical solution of these equations to a direct computation of the free energy from the hamiltonian. Particular attention is paid to the relation to the Heisenberg xxx spin chain. We conclude in Section 7 by summarizing our results. The appendices cover some basics on Weierstraß elliptic functions (Appendix A), a thorough analysis of the important function φ (defined in equation (12)) (Appendix B) and finally an analysis of convergence of solutions to the Bethe ansatz equations (Appendix C).

Inozemtsev's elliptic spin chain
Inozemtsev's elliptic spin chain with spin 1/2 as defined by the Hamiltonian in equation (1) has been studied extensively (see e.g. [12,15,26]). It is expected to be integrable, although this has not been completely proven since there is no proof that the found conserved quantities actually commute. There does exist a set of exact eigenfunctions in the form of a generalized Bethe ansatz and transcendental equations that determine the quasi-momenta. Various models can be reached starting from the elliptic spin chain by varying either the parameter κ and/or the length L of the chain. All of these spin-1/2 models are characterized by hamiltonians of the form where the potential V can depend on the length L of the chain, which is possibly infinite. In this work we will focus solely on the ferromagnetic case J > 0. Following [19,27], we accommodate these 1 In the paper [22] a study of the thermodynamics of Inozemtsev's spin chain was announced, but to the author's best knowledge this study has never been published. 2 An application of the method by Yang and Yang to Inozemtsev's spin chain can be found in the author's [25].
In that unpublished work one can also find an account of part of the results discussed in the present work. limits by redefining the hamiltonian (1) by rescaling and shifting the potential by site-independent factors: from now on we take Inozemtsev's elliptic spin chain to be defined by the hamiltonian where where ζ L is the Weierstraß ζ-function with quasi-periods (L, iπ/κ). If one sends κ to infinity we reach the Heisenberg xxx spin chain (see [12] or Appendix A of [19]) with potential If one sends κ to zero, one obtains the hamiltonian of the Haldane-Shastry (HS) spin chain with potential [13,14] V (L) On the other hand, if we keep κ fixed and send L → ∞ we reach Inozemtsev's infinite-length spin chain with potential which was treated extensively in [16]. All limits are summarized in Fig. 1.
The exact solution of Inozemtsev's spin chain of infinite length is based on the su(2)-invariance of its local hamiltonians allowing for the M -particle ansatz where W denotes the set of all m ∈ Z M such that 0 ≤ m i ≤ M − 1 for all 1 ≤ i ≤ M and S M is the symmetric group of M symbols. The coefficients c m 1 ···m M (p) can be found solving the set of equations where Z n,n = {k ∈ Z | max(−n, n − M + 1) ≤ k ≤ min(M − 1 − n, n )}. These eigenfunctions are closely related to the eigenfunctions of the continuous Calogero-Moser-Sutherland model with 1/ sinh 2 -interaction [16]. The associated eigenvalues are additive, the energy of an M -magnon state being given by with (see also Fig. 2) where the Weierstraß functions ℘ = ℘ 1 and ζ = ζ 1 are defined on the lattice with periods (1, iπ/κ). Note that, unlike the finite-length case, solving the eigenvalue problem with this ansatz does not lead to any restrictions on the quasi-momenta and one needs to resort to other methods to find the spectrum of the model. A way to introduce Bethe equations is to follow the asymptotic Bethe ansatz scheme (ABA), which can be summarized as imposing periodic boundary conditions on the asymptotic form of the eigenfunctions (8) [27,28]. This leads to Bethe equations (BE) (see [28]) where M denotes the total number of magnons and the meromorphic function φ is given by Solving these equations at L → ∞ yields sets of quasi-momenta that are good candidates for parametrizing the spectrum of Inozemtsev's infinite-length spin chain, but one needs to to verify this by different means since usually the relation between quasi-momenta and physical states is not one-to-one. The solutions to (11) might also be used to study the thermodynamic limit (M, L → ∞ with M/L fixed) of Inozemtsev's elliptic spin chain, since at very large L the eigenfunctions of the elliptic spin chain can be approximated by those of the infinite spin chain, as was shown equivalently in the su(1|1) case in [19]. 3 The system of equations (11) is the usual form of BE, where for example the Bethe equations of the homogeneous Heisenberg xxx spin chain (BE xxx ) are of this form with the function φ replaced by although in that particular case p can be replaced by φ xxx in the BE xxx altogether due to the form of φ xxx . Note that lim κ→∞ φ = φ xxx , implying that the BE xxx can be found from equation (11) by taking this limit. It is therefore natural to expect that the solutions to the BE (11) with (12) are closely related to the known results for the BE xxx .

Solving the Bethe equations asymptotically
We are interested in solving the system of M equations (11) for an M ∈ N in the limit L → ∞ for sets of noncoinciding 4 complex momenta {p j } ∈ D, with We can restrict −π ≤ Re(p j ) < π for all 1 ≤ j ≤ M due to the translation invariance of the spin chain. The total momentum and energy of these sets should be real, that is The behaviour of the terms in (11) in the limit L → ∞ is quite simple: The left-hand side only depends on s j = sign(Im(p j )), Depending on the exact form of φ, this implies that different sets of momenta might correspond to a single set of minimal data {(θ j , s j )} j≤M with a sign s j ∈ {+, 0, −} to indicate the sign of the imaginary part of the associated momenta. To analyze the possible solutions as structured as possible, we will first characterize the allowed sets of minimal data by the usual analysis for string solutions [1,23].
Consider a case in which s 1 = +. We see from the BE for j = 1 that there must be an n ≤ M (and we can take n = 2 without loss of generality) such that θ 1 − θ 2 → −i as L → ∞, such that also the right-hand side of the BE converges to zero. We will say colloquially that (θ 2 , s 2 ) helps (θ 1 , s 1 ) to satisfy its BE. It means that in the limit the real parts of θ 1 and θ 2 coincide and that Im(θ 2 ) = Im(θ 1 ) + 1. There are three options for s 2 . If s 2 = +, the reasoning continues along the same line until we find an s j with either s j = − or s j = 0. We will not treat cases with s j = 0 here for simplicity, since they can be derived from our results without much work, but show a possible configuration on the right in Fig. 3 nevertheless. If s 2 = −, however, we see that the Bethe equation for j = 2 is already satisfied due to the presence of (θ 1 , s 1 ). Therefore, we do not need to add more tuples (θ n , s n ) to the set to make it consistent (provided that the reality condition on momenta and energy are satisfied). By carrying out a similar reasoning for the case s 1 = −, we see that the basic structure of a set of minimal data is a string of pluses and minuses as in Fig. 3. From this analysis, we see that an allowed set of points {θ j } j≤M should be a subset of for a certain m ≤ M and certain θ R,I ∈ R. Note that we allow several θ to occupy the same point in θ-space, as seems allowed by the above analysis: as long as there is another tuple (θ n , s n ) that provides the correct limiting behaviour as L → ∞, we can include any tuple we want. In particular, if (θ 2 , s 2 ) is such that it helps (θ 1 , s 1 ), it will also help any (θ n , s n ) that satisfies θ n = θ 1 and s n = s 1 . We will see in Section 3.3.1 that as long as the basic structure is present, we can almost freely associate as many tuples as we want to a single point in θ-space. These solutions cannot as easily be depicted as in Fig. 3. Of course, we are at the moment ignoring possible issues with convergence, which we will address in Section 4.3. Also note that actual solutions to the BE should in the end have real momentum and energy. However, given a set of momentum associated to a set of minimal data {(θ j , s j )} j≤M as derived above, we can always add the complex conjugates of these momenta to the set to make sure that both total momentum and energy are real, as long as {φ(p j )} ∩ {φ(p j )} is either empty or consists of one real element (see the middle and right configuration in Fig. 3 respectively). This is possible due to the meromorphicity of φ and the one-particle energy .
The real question now is whether this very general analysis (and in fact more general than usually considered) is even necessary in the present case. Before going into details about this question, let us first make the connection with the known results for the xxx spin chain and see why we do not need this general approach in that case. R ⊂ C is indicated by the dashed line and the arrows indicate the structure that solves the BE: an arrow from sign s m to sign s n indicates that the BE with j = n are satisfied because of the presence of (θ m , s m ) on its right-hand side (so (θ m , s m ) helps (θ n , s n )). The left configuration is the standard string solution as the ones occurring for the xxx model. The middle configuration is a new feature of Inozemtsev's BE and consists of two connected components. For odd M , the allowed sets of minimal data look like the right configuration, with a real momentum in the middle, indicated by the 0.

Solutions for the xxx spin chain
The Bethe equations for the Heisenberg xxx spin chain are where φ xxx (p) = 1 2 cot p 2 , which is in this case usually called the rapidity function. The structure of the solutions to these equations is very simple. For each M , there exists a one-parameter family of string solutions of length M , which can be most conveniently parametrized in terms of the rapidities λ j = 1/2 cot (p j /2) and is given by The reason for this simple structure is the bijectivity of φ xxx as a function on D. Following the reasoning introduced in the previous section, we want to solve equations of the form φ xxx (p 2 ) = φ xxx (p 1 ) ± i (where p 2 is the unknown). By the bijectivity of φ xxx , these equations have unique solutions, which leads to a unique set of momenta as soon as p 1 is fixed. Additionally, the sum of momenta must be real to ensure that the energy of the solution is real, which imposes that the rapidities have the prescribed imaginary parts given in (19). So due to the bijectivity of φ xxx , all the asymptotic solutions to the BE xxx are of the form as in (19) and usually called string solutions. This is no longer the case if the bijectivity of φ is lost, which turns out to be the case for Inozemtsev's spin chain. The function φ appearing in Inozemtsev's BE is odd and quasiperiodic, satisfying

Behaviour of
which means that its behaviour on the region completely determines its behaviour on D. One can prove using the argument principle that φ : D ≤κ → C is almost bijective. 56 φ is certainly surjective, but it attains twice those θ ∈ C for which Im(θ) = ±1/2 and −θ crit < |Re(θ)| < θ crit , where θ crit > 0 depends on the parameter κ and is defined by where p crit is the unique solution on [0, π] to the equation The preimages of these θ's lie on the top and bottom boundary of D ≤κ , i.e. where Im(p) = ±κ. This behaviour is illustrated in Figure 5. The quasi-periodicity and almost bijectivity of φ when restricted to D ≤κ inspires to introduce a partition of D into regions such that φ is bijective when restricted to such a region: the fundamental region D f is defined as and the region D n is defined as the region obtained by shifting D f by 2κin, that is D n = D f +2κin, as can be seen in Fig. 4. The partition {D n } n∈Z of D is such that the restrictions φ n to D n are bijective functions onto C. This will make it easier to categorize the momentum sets that belong to a certain set of minimal data {(θ j , s j )}. Finally, note that there exists exactly one other partition consisting of connected sets that differs from ours, which can be created by mirroring this partition in the real line.

Solutions for Inozemtsev's spin chain
The fact that φ is so far from being injective has great consequences for the solutions of the BE of Inozemtsev's spin chain. Equations of the form 5 With almost bijective we mean that there exists a restriction of φ to a domain differing from D≤κ by a set of measure zero that is bijective. 6 This is shown in Appendix B.
for a given θ ∈ C, have an countably infinite set of solutions, parametrized by the region index in which each of the solutions lies. In particular, the equation (24) has solutions for p with positive and with negative imaginary parts. This makes it possible for a string solution to consist of two parts, as in the middle of Figure 3: each part in itself forms a consistent solution to the BE, but only the sum of the parts has real energy and momentum. These new solutions also have more degrees of freedom than the usual string solutions: whereas the usual string solutions (such as the left configuration in Fig. 3) have no freedom in choosing the imaginary parts of the θ's, the new solutions can be shifted in the imaginary direction as long as the two parts remain complex conjugate and distinct. More precisely, for m distinct θ j we can choose θ I parametrizing the imaginary part of the θ j as in (17) to be anything from the set As an example, consider the solution consisting of the four momenta for the case where κ = 1.26. It consists of two connected components and has m = M (i.e. noncoinciding θ j ), θ I = 1.8 and θ R = 0.6. This is just one of the countably infinite number of solutions specified by these θ R , θ I : there are infinitely many D n from which p 1 can be chosen and the same is true for p 2 . We see that solutions with m = M are not a one-parameter family (as was the case for the xxx spin chain); there are 2 continuous parameters and M/2 discrete ones needed to specify an M -momenta solution of this type.

Solutions with coinciding θ j
Although the sets of minimal data considered in the previous section are already an extension to the usual string analysis, Inozemtsev's BE allow even more general sets. The fact that we are only interested sets of non-coinciding momenta does not mean that also the set of θ j should be non-coinciding. The non-injectivity of φ precisely allows us to associate any number of momenta to any particular value θ. Moreover, in many cases we can associate momenta with both positive and negative imaginary part to each value. In order to be able to characterize these sets of minimal data more easily, we will no longer allow the θ j to be coinciding, but instead associate multiple signs s j,i j to a single θ j , that is we rewrite minimal data where now θ j = θ n implies j = n. We can depict these sets of minimal data, which we will call coincident minimal data, by placing the s j,i j belonging to the same θ j on a horizontal line (the level ). In this way, we can depict a set of coincident minimal data as has been done in Figure 6, where level j contains P j pluses and M j minuses and the total number of levels is m ≤ M . These numbers satisfy Thus in Figure 6, there are M 1 momenta with negative imaginary part associated to the image point θ R +θ I i with biggest imaginary part, M 2 momenta with negative imaginary part to the image point θ R + (θ I − 1)i and P 2 momenta with positive imaginary part, etc. In this configuration, a sign s j,i j = + on level j receives help from all the s j+1,i j+1 , whereas a sign s j,i j = − receives help from all the s j−1,i j−1 . Sets of coincident minimal data can be parametrized as follows: we fix an integer M and an m ≤ M and choose an allowed sign configuration conform Figure 6. 7 Then we choose θ R ∈ R and θ I ∈ R M . An actual solution to the BE, a coinciding solution, also requires the choice of a region D n for each of the M signs in our configuration. Generically, there is an infinite number of allowed regions. It is clear that these solutions enjoy even more freedom than the ones considered in the previous sections.
However, one might argue that for coinciding solutions we can no longer follow the naive construction of string solutions, since in this case taking the limit L → ∞ becomes more problematic; indeed this is the case, because many different terms might converge or diverge on the right-hand side of a given BE for one of the momenta of a coinciding solution. In Appendix C we address this question more carefully, but we will see in the rest of our analysis that the question whether or not one can still take the limit is irrelevant. with energy E 8 = 0.211. As shown in Appendix C, this solution does have a defect: there is no consistent way to consider the limit L → ∞ for this solution, implying it is not a good candidate to parametrize the spectrum.

Pruning the solution set
The solutions presented in the previous sections obey the rules that are usually obeyed by string solutions of Bethe equations, such as the ones for the xxx model. Some of the features of the new solutions do raise questions: the solutions have too many degrees of freedom to execute the usual string hypothesis program, i.e. assume that the string solutions accurately describe the thermodynamic behaviour of the model and derive thermodynamic Bethe ansatz equations. In particular, excluding momenta in favour of θ's would mean having to introduce an uncountably infinite number of types of particles. We might expect, however, that these issues are due to the fact that only a subset of our set of solutions contains information about the spectrum of Inozemtsev's spin chain and most of the solutions to the BE are actually non-physical: they have some sort of defect that forces us to discard them as physical solutions. This is indeed the case, there are four main types of defects to be found in our set of solutions: 1. the associated wavefunction does not vanish at infinity, that is the momenta do not parametrize a bound state.
2. they do not correspond to a string solution to the Heisenberg xxx spin chain in the limit κ → ∞.
3. there is no consistent way to consider the limit L → ∞.

the associated wavefunction is identically zero.
In the following sections we will consider these defects and discard the solutions that suffer from these defects.
Unfortunately, the complicated form of the wavefunctions for M > 2 makes it difficult to prove a similar statement for bound states consisting of more than 2 particles, although numerical analysis of the wavefunctions shows that the statement seems to be true at least up to M = 6.

Relation to the Heisenberg xxx spin chain
To get more evidence for the fact that the solutions with momenta outside D ≤κ are non-physical, we take a closer look at the relationship between Inozemtsev's infinite spin chain and the xxx spin chain. Let us therefore first consider a general solution to Inozemtsev's BE: it consists of a set of θ j , an assignment of momenta to each of the θ j and the regions where each of these momenta can be found. If we take the limit towards the xxx spin chain (κ → ∞), all the momenta that do not lie in D ≤κ acquire infinite imaginary part, since they lie on the outside of D ≤κ which fills up D entirely in this limit. Comparing this to the allowed string solutions for the Heisenberg xxx chain (19) and stipulating that all solutions should flow to a xxx string in the limit κ → ∞ also suggests we should abandon solutions that have momenta outside of D ≤κ .

Convergence issues
The arguments above potentially reduce the solution set enormously, but their origin lies in applying restrictions that are not related to the solving of the BE itself. Interestingly, there is a subtle issue arising because some solutions have coinciding θ j , which forces us to look more carefully at the procedure of taking the limit L → ∞ when we are looking at the BE. For a standard string solution, there are usually only two terms on the right-hand side of the BE of p j that do not have a finite limiting value, but converge to zero or diverge. One of the two terms is there to make sure that the equation is satisfied in the L → ∞ limit, but the other one actually counteracts this. By associating to each momentum a speed with which its limiting value is reached, one can take the limit in a consistent way. For solutions with coinciding θ j however, the case is more complicated, because more terms influence the limit. In Appendix C, we discuss this matter in detail and show that there is indeed a consistent way to consider the limit for standard string solutions and also give an example of a coinciding solution for which there is not. This excludes some of the coinciding solutions from being in the spectrum, but most of them remain to be candidates. In particular, we cannot exclude more solutions than can already be be excluded using the previous two arguments.

Vanishing wavefunctions
If we restrict the domain of our momenta to D ≤κ , we have almost completely excluded the possibility of coinciding θ j . In fact, the only remaining solutions of this kind must be built up from momenta living on the boundary of D ≤κ , since only on the boundary is φ still non-injective (see Fig. 5). So for |θ R | < θ crit there are 2 solutions of the equation which we name p 1 and p 2 . In this way we can build four two-particle bound states, by pairing the momenta as follows: The bound states on the first line are of the form p 1 = p − iκ, p 2 = p 1 . A simple computation shows that the wavefunctions of these bound states vanish identically. Numerical analysis up to M = 6 suggests that this holds in general for wavefunctions parametrized by a set of momenta for which p n = p m + 2κi for some n, m. If we assume this is indeed true, we can no longer built coincident solutions. However, we can still build two types of particles out of momenta in the region D ≤κ for |φ R | < θ crit and even M : there are two two-particle bound states (see the lower line of equation (30)) that one can use as a basis for building a solution, after which there is no choice for the remaining particles, since they are fixed by the requirement of real energy. However, one can check that these two types of bound states have exactly the same energy and total momentum and in fact parametrize the exact same wavefunction; they parametrize the same bound state. This situation is very similar to the one encountered in [29] and our solution is the same: we only keep one of the two. Although the choice is arbitrary, we can choose by restricting the allowed momenta domain even further, from D ≤κ to D f . Choosing the other bound state amounts to partitioning the domain D in the alternative way as described in Section 3.2.
Note that the two-particle bound states we have found here are actually very peculiar: they are not self-conjugate when viewed in momentum space, which is a novel feature of Inozemtsev's spin chain (see also [18]), but the total momentum and energy of this bound state are real. In fact, this is even more interesting when one realizes that the self-conjugacy of the string solutions in the spectrum of the Heisenberg xxx model can be traced back to the underlying algebraic structure [30]. Despite this difference, we will see in the next section that the equations describing the thermodynamic behaviour of Inozemtsev's spin chain bear a striking resemblance to those for the xxx model.

Remaining solutions
All of the arguments presented above indicate that we should only consider solutions built up out of momenta from the fundamental region D f . As another check to see that we are on the right track we have plotted in where E M is the energy of an M -string lying in the fundamental region as a function of total momentum. We see clearly that the inequality is satisfied for all plotted M and numerical analysis shows this is true at least up to M = 40. This implies that the solutions are indeed bound states for positive J, because the energies of these states are smaller than the sum of one-particle energies. The next step now is to invoke the string hypothesis and assume that the remaining solutions accurately describe the thermodynamic behaviour of Inozemtsev's spin chains. That will allow us to perform the thermodynamic Bethe ansatz program to derive equations that describe the free energy of Inozemtsev's spin chains, which we will do in the next section.

Thermodynamic Bethe ansatz
Since we are interested in the thermodynamic regime of Inozemtsev's elliptic spin chain, we want to send the number of quasi-particles M and the length of the chain L to infinity, while keeping the ratio M/L fixed. The well-known method we will deploy here to take this limit for the BE (11) is called thermodynamic Bethe ansatz (TBA) and the resulting set of equations are usually called TBA-equations.
Let us first summarize which string solutions we consider to be relevant in the thermodynamic limit: since we now have a one-to-one relation between p's and θ's, we can no longer create solutions with coinciding θ j , but also lost our freedom to choose θ I . Our restricted φ is meromorphic and bijective, implying we can no longer make strings like the one portrayed in the middle of Fig. 3, because for all p ∈ D f we have that sign(Im(p)) = −sign(Im(φ(p)). Therefore, the remaining strings are of the form which we will call Q-strings. Their total momentum and energy are given by where is the one-particle energy given in (10) and unfortunately no explicit formula for φ −1 is known. As any momentum set {p j } must be built up from bound states, we can fuse the BE for the composite particles parametrized by our Q-strings: where the N Q ≥ 0 denotes the number of Q-strings and where Taking logarithms in (35) we get c P (θ P,l )L = I P,l , where the I P,l are integer quantum numbers and are the counting functions. We can check numerically that these functions are monotonically increasing. Now we can introduce particle and hole densities ρ Q ,ρ Q that should satisfy In the limit L → ∞ the counting functions get transformed as the summations become integrals: We define the convolution Taking the derivative explicitly, we see that (40) becomes ρ P (θ) +ρ P (θ) = 1 2π where we have used the kernels K P (θ) = 1 π P P 2 + θ 2 for P ≥ 1 and K 0 (θ) = δ(θ) Note that these kernels are exactly the same as those appearing in the derivation of the TBAequations for the xxx spin chain. In fact, our entire derivation differs from that one only because our formulae for p Q and E Q cannot be written in terms of elementary functions. By varying (43), we get where we sum over the repeated indices. Now we can introduce a free energy density and find the point of thermodynamic equilibrium: The entropy density is defined as Varying, substituting equation (45) and changing integration variables in the kernel term we end up with This directly leads to the TBA-equations of Inozemtsev's spin chain, which when we introduce the which are of exactly the same form as those for the Heisenberg xxx model; the only difference sits in the definition of the energies E Q . Moreover, after sending κ → ∞ we recover the TBA-equations for the xxx spin chain.

Free energy
Much of the information about the system in thermal equilibrium can be extracted from the density of the Helmholtz free energy f . We can express the free energy solely in terms of Y -functions as follows: plugging in our definition of Y 's into our definition of f (46) leads to Now, using equation (45) we can replaceρ Q by ρ Q : Using the fact the the Y -functions should obey the TBA-equations (50) we are left with the expression In the next section, we will use this expression for the free energy to compare our equations describing the thermodynamics of the model with a more straightforward method starting from the hamiltonian. However, let us for completeness first mention how one could simplify the TBAequations further.

Y -system
One can find a simpler-looking set of equations using the function s(θ) = 1 4 cosh πθ 2 (54) together with its pseudoinverse s −1 defined by and the property that s (K P −1 + K P +1 ) = K P for P ≥ 1.
These equations, the Y -system for Inozemtsev's spin chain, take the following form: with M > 1 in the second line and where the superscripts ± indicate shifts of ±i and the inversê K 1 −1 means thatK The Y -system of Inozemtsev's spin chain reduces to the Y -system for the xxx spin chain in the by now well known limit κ → ∞ given by but for finite κ the exponential prefactors do not seem to simplify.
The Y -system looks simpler than the TBA-equations, but also has a downside: it admits more solutions than the ones we are interested in alone. One has to complement the Y -system with a set of asymptotics and possibly other analytic data to find the solution describing the thermodynamics of Inozemtsev's spin chain. Therefore we will later use the TBA-equations (50) to find numerical results. One important piece of data one can extract from the Y -system is the value of the Yfunctions as T → ∞: in this limit the exponential prefactors become unity, leaving us with the Y -system (59). From (50) we see that in this limit the Y -functions should become constant functions and plugging this information into the Y -system gives the asymptotic result which one can use to solve the TBA-equations numerically.

Solving the TBA-equations numerically
Obtaining numerical results 8 for Inozemtsev's spin chains using the TBA-equations (50) is interesting for several reasons: firstly, we can use numerical results to check our equations, thereby implicitly corroborating the usage of the string hypothesis, used to the TBA-equations, as well as checking our treatment of the solutions of Inozemtsev's BE in Section 4. Secondly, with numerical results we can check the hypothesis put forward in [19] that the normal and supersymmetric versions of Inozemtsev's spin chains are equivalent in the thermodynamic limit.
To perform these checks we have gathered three types of numerical data: (1) the free energy for the xxx and Inozemtsev's spin chains from TBA-equations 9 , (2) the free energy for the xxx, Haldane-Shastry and Inozemtsev's spin chains computed from the finite-length hamiltonians and (3) the free energy for the supersymmetric versions of these three types of spin chains.

Methods
(1) Solving the TBA-equations can be done quite fast using Fast Fourier Transform to compute the convolutions in a Picard iteration scheme. To be able to this we have cut the number of Yfunctions to not more than 35 and treated the real line as a grid of (typically) 2 8 , 2 9 points. In our particular case the tricky part is computing the energies E Q efficiently, because their definitions contain the inverse of φ, for which no explicit formula exists. We have written a program that is capable of finding inverses numerically, but finding inverses is still a time-consuming task compared to performing the iterations.
The iteration scheme takes a set of Y -functions and iterates using the TBA-equations until stability is reached, that is until the biggest pointwise difference between ingoing and outgoing Yfunctions is smaller than 10 −10 . Then the number of Y -functions is increased and a stable solution of this set of Y -functions is found. For both sets of Y -functions the free energy is calculated and if the relative difference between the found free energies is smaller than 10 −5 J we declare the solution stable and otherwise keep increasing the number of Y -functions. This approach is very similar to the one used in e.g. [32,33] and following these authors we believe that our results should be accurate at least up to a few percent. In particular, we have also explicitly computed the free energy from the TBA-equations for the Heisenberg xxx spin chain.
(2) We have also compared the results from the TBA-equations with another calculation of the free energy per site: given a spin chain hamiltonian H for finite L, one straightforwardly derives that where L is the length of the spin chain and the subscript reminds us how we calculated this free energy. As long as L is not too large (typically L ≤ 15) one can perform these matrix operations explicitly reasonably fast. We have done this for the xxx, HS and Inozemtsev's spin chains, using the finite-length hamiltonians with potentials (5), (6) and (3) respectively for increasing L until the results stabilized (in this case that means that consecutive terms differ by less than 1 − 2%). Extrapolating the relative difference suggests that the results are also accurate up to a 5 percent.
(3) Finally, we have reproduced the results in [19], giving the free energy per site of the supersymmetric versions of the xxx, HS and Inozemtsev's spin chains. Computing the relevant integrals numerically was done with a high degree of precision (up to 50 digits). We see that as κ increases the free energy of Inozemtsev's spin chain converges to the free energy of the xxx spin chain. Also, for decreasing κ the free energy approaches the free energy of the HS spin chain. This shows that our TBA-equations reproduce the thermodynamic behaviour of the two limiting spin chains in the appropriate limits and nicely interpolate between them at finite κ.

Results
In Fig. 10 we have plotted free energies of our three models as calculated from (61) and the TBA-equations when relevant, accompanied by the free energy of the supersymmetric version of these models. We see that all the free energies agree to very high accuracy for T 5J, whereas deviations occur for smaller T . The differences between the different functions scale as J/T , as can be confirmed by repeating the analysis for different values of J. Moreover, these deviations occur for all our models, including the Heisenberg xxx and Haldane-Shastry spin chains for which previous studies have confirmed the correctness of the underlying equations [19,31,34]. The deviations are most likely caused by numerical inaccuracies related to exponentiating large numbers (as happens in equation (61) and in the TBA-equations (50)) and to restricting the real line to a finite interval. We estimate the observed discrepancies to be within the error of these numerical effects. Therefore, we regard the data in Fig. 10 as confirmation that our TBA-equations (50) truly determine the thermodynamic behaviour of Inozemtsev's spin chains. In particular, this validates our usage of the string hypothesis in the derivation of the TBA-equations which is non-trivial in itself. Moreover, the matching with thermodynamic data of the supersymmetric models confirms that the hypothesis brought forward by Finkel and González-López in [19] that the supersymmetric version coincides with the non-supersymmetric model in the thermodynamic limit. It would be interesting to see whether it is possible to derive the defining equation for the free energy of the supersymmetric model from our TBA-equations (50), perhaps providing more insight into why this correspondence between certain models and their supersymmetrization exists.
One can further check the claim that Inozemtsev's finite-length model (3) and infinite-length model (7) coincide in the thermodynamic limit by computing the free energy of the finite-length spin chain given by the hamiltonian (2) with potential using (61). The resulting free energy coincides to such high accuracy with the result obtained using the hamiltonian with potential (3) that they would not be separately discernable in the plots in Fig. 10. This, combined with the fact that the limits to the xxx and HS spin chain behave as expected provide additional evidence that our TBA-equations (50) are correct.

Conclusions
In this paper we have investigated the thermodynamics of Inozemtsev's elliptic spin chain. Starting from the Bethe ansatz equations of Inozemtsev's infinite-length spin chain, we classified all the solutions to these equations that have real energy and total momentum. We then analyzed whether they parametrize bound-state solutions of Inozemtsev's infinite-length spin chain, for example by comparing the results with one of the limiting cases, the infinite-length Heisenberg xxx spin chain. This reduces the number of solutions immensely, leaving a set of solutions that is structurally very similar to the string solutions of the Bethe ansatz equations of the xxx spin chain. One interesting new feature is the presence of solutions with non-selfconjugate momenta. Carrying out the thermodynamic Bethe ansatz program we have derived a set of coupled integral equations and an associated set of finite-difference equations (Y -system) that allows one to compute the free energy of the model at thermal equilibrium. We have solved the integral equations numerically and compared them with the free energy computed directly from the finite-size hamiltonian as well as with the free energy of its limiting models, the Heisenberg xxx and Haldane-Shastry (HS) spin chains. All the results seem to be consistent, corroborating the correctness of our derived integral equations. Moreover, we also compared the free energy of Inozemtsev's spin chain with the free energy of the supersymmetric version of this spin chain obtained by Finkel and González-López and concluded that these models coincide in the thermodynamic limit. Our findings extends the relationship between Inozemtsev's spin chain and the xxx and HS spin chains to the thermodynamic regime, in line with the finding that this is also true for the supersymmetric version of these models [19]. One might wonder whether similar relations exist for other generalizations or deformations of such spin chains. Also, further research could be conducted to see whether one can analytically show the equivalence of the normal and supersymmetric Inozemtsev spin chain in the thermodynamic limit as has been done for their limiting models. It would be interesting to get a better understanding of the relation between our Y -system for Inozemtsev's elliptic spin chain and Y -systems for related su(2)-invariant models, for example because it might lead to an elliptic extension of the kernel identities we used to derive the Y -system. Moreover, it would be interesting to see whether one can simplify the Y -system even further in light of the recent advances in simplifying the Y -system for N = 4 super Yang-Mills theory to what is known as the quantum spectral curve [35].
−κi κi π −π 0 Re(p) Im(p) Figure 11: The contour around which we integrate to find the number of zeroes in D ≤κ for φ.
edge of this contour, the imaginary part of φ is constant. Let θ ∈ C be arbitrary, but such that φ(z) = θ has no solutions when φ is restricted to the contour. Then the functionφ(z) = φ(z)−θ has no zeroes or poles on the contour and we can use the argument principle to state that where N is the number of zeroes and P the number of poles ofφ in the interior of the contour, which is the fundamental region of φ. In this case, we have P = 1. We can calculate the integral on the left-hand side: the contributions from the vertical parts of the contour cancel each other due to the periodicity ofφ.
For the contributions of the top part, we see the following: becauseφ is 2π-periodic in the real direction. Note that we could evaluate the integral using the logarithm, because we know that the imaginary part ofφ is constant along the path, allowing us to find a holomorphic branch for the logarithm on a neighbourhood of the top part of the contour. In a similar fashion, one can show that the contribution from the bottom part vanishes, thus we end up with implying that for all the θ we considered,φ has exactly one zero in the fundamental region, thus φ(z) = θ has exactly one solution in this region.
On the boundary of D ≤κ , the following holds. The restriction x → φ(−π + ix) (with x ∈ [−κ, κ]) has negative derivative everywhere. Moreover, since φ(−π ± κi) = ∓i/2, we can conclude that this restriction maps bijectively onto [−i/2, i/2]. This shows that φ : [−π, π[ ⊕ ]−iκ, iκ[ → A ⊂ C maps bijectively onto its image A. On the top part of the contour we can write x → φ(x + iκ) for the restriction. A plot of this function is shown in Figure 5, which shows that this restriction is not bijective onto its image. In fact, all image values are attained exactly twice. We call the graphs maximum θ crit and by symmetry, its minimum is −θ crit . The value of p for which Re(φ(p+iκ)) = θ crit we call p crit . By symmetry, the minimum is attained at −p crit . The behaviour of the real part of φ along the bottom boundary is exactly the same.
Bethe equation of the momentum associated to this plus sign reads where we have θ j,n j = φ(p j,n j ) and θ k belongs to level k. Note that the terms belonging to other momenta on level j are all 1 and are not written explicitly and that the θ j,n j do not actually depend on n j . As L → ∞, the left-hand side converges to 0. Most of the terms on the right-hand side converge to finite values and are irrelevant for the analysis. The interesting terms are those belonging to level j ± 1. They form the product However, to each momentum we have associated a convergence rate and we can let all the fractions in this product converge to their limiting value with different rates. In the infinite-L limit, the term belonging to p j+1,γ j+1 (with γ = α, β) on the level j + 1 behaves as while the term belonging to p j−1,γ j−1 behaves as θ j,n j − θ j−1 + i θ j,n j − θ j−1 − i ≈ O exp min δ j,n j , δ j−1,γ j−1 L .
From now on, we write (x, y) := min(x, y). In total, the product of terms belonging to level j + 1 converges as δ j,n j , δ j+1,β j+1     and combining this with the similar result for the level j − 1 we see that the right-hand side of the Bethe equation (68) behaves as δ j,n j , δ j+1,β j+1 and therefore goes to zero only when the convergence rates obey δ j,n j , δ j+1,β j+1 δ j,n j , δ j−1,β j−1 < 0.
In a similar fashion, one can derive that the Bethe equation corresponding to a momentum p j,n j with negative imaginary part is satisfied only when δ j,n j , δ j+1,β j+1 δ j,n j , δ j−1,β j−1 > 0.
For a valid solution of the Bethe equations, equation (73) must be satisfied for all plus signs, while equation (74) must be satisfied for all minus signs. Note that these restrictions arise simply because there is more than one term that exhibits vanishing or divergent behaviour and we should include more information to find the behaviour of the product. This problem already exists in many of the previously known cases (such as the Hubbard model or the Heisenberg xxx model), but as far as we know, this has never been addressed. Fortunately, however, the restrictions (73),(74) simplify drastically for the usual simple string solutions occurring in the aforementioned cases and can easily be solved. The system of restrictions for a string solution without a real momentum involved read where the superscripts indicate the sign of the imaginary part of the associated momenta. It is solved by the ordering δ 1 < δ 2 < · · · < δ mp = δ mp+1 > δ mp+2 > · · · > δ M , together with δ j := δ + j = δ − j . However, determining whether the system of equations consisting of (73) and (74) for a general tree solution can be solved is a much more complicated question. In the next section, we treat some cases and include an example from which it follows that not every sign configuration has a consistent set of convergence rates.

C.2 Examples
A coinciding solution consists of at least 2 levels. The 2-level case (illustrated in subfigure (a) in Figure 12) can also be solved in general, because the inequalities are trivially satisfied. However, − − · · · · · · − + + · · · · · · + M 2 > 0 already the 3-level case harbours an example of a configuration that cannot have a consistent set of convergence rates. Consider the example in subfigure (b) in Figure 12. The relevant set of equations is (δ where we omit the superscript (±) when it is not necessary. We first try to deduce which of the δ's should be the smallest one of these four. From the upper equation, we conclude that neither δ (+) 2 nor δ 3 can be the smallest, while the lower equation tells us that neither δ (−) 2 nor δ 1 can be the smallest. Therefore, none of the 4 rates can be the smallest, thus no solution can exist. Note that this example can be extended: if we include P 2 > 0 pluses and M 2 > 0 minuses on level 2, the resulting set of restrictions has the system (77) as a subsystem and cannot be solved. In particular, this shows that example (b) we treated in Section 3.3.1 is not a valid solution to the BE after all, although we could find momenta to match the configuration. Moreover, any sign configuration that contains this 3-level structure cannot be solved. However, all other 3-level configurations do admit a consistent solution as a careful analysis of the cases shows.
We have not been able to find a general algorithm to solve these complex coupled sets of inequalities or prove the existence (or absence) of a solution. The only configurations we found that lead to inconsistent inequalities are of the type described in the previous paragraph. In any case, the structure of the solutions is complicated, but in the present analysis we do not need it.