Quantum Stability for the Heisenberg Ferromagnet

Highly spinning classical strings on RxS^3 are described by the Landau-Lifshitz model or equivalently by the Heisenberg ferromagnet in the thermodynamic limit. The spectrum of this model can be given in terms of spectral curves. However, it is a priori not clear whether any given admissible spectral curve can actually be realised as a solution to the discrete Bethe equations, a property which can be referred to as stability. In order to study the issue of stability, we find and explore the general two-cut solution or elliptic curve. It turns out that the moduli space of this elliptic curve shows a surprisingly rich structure. We present the various cases with illustrations and thus gain some insight into the features of multi-cut solutions. It appears that all admissible spectral curves are indeed stable if the branch cuts are positioned in a suitable, non-trivial fashion.


Introduction
The Heisenberg magnet [1] is one of the very first quantum mechanical models, and it serves as the prototypical spin chain. Although it was set up almost 80 years ago it still remains a fascinating subject with many features left to be understood. Only three years after the discovery of the model and owing to its integrability, Bethe was able to write a set of equations (1.1) which determine the complete spectrum [2]. The terms "Bethe equations" or "Bethe ansatz" later became synonymous for the exact solution for generic integrable spin chain models.
Although the Bethe equations describe the complete and exact spectrum, it is virtually impossible (and perhaps not very enlightening) to solve them concretely for generic states somewhere in the middle of the spectrum. Nevertheless some corners of the spectrum are accessible (and interesting). In particular these are the low-energy and highenergy states for very long chains (the thermodynamic limit), where the Bethe equations are approximated by integral equations. Most studies have focused on the low-energy spectrum of the antiferromagnet and this regime is well-understood, see [3] for a review. For example, the antiferromagnetic state is a solution to an integral equation [4] and its excitation quanta are called spinons [5].
Conversely, the low-energy regime of the ferromagnet (or equivalently the high-energy regime of the antiferromagnet) is much less explored. The ferromagnetic ground state coincides with the vacuum of the Bethe ansatz, it is trivial. The excitation quanta are called magnons, and each magnon corresponds to a single Bethe root u describing the rapidity of the magnon. Magnons can form bound states which are usually called Bethe strings. In a k-string centred at rapidity u there are k Bethe roots u 1,...,k arranged in a regular pattern u j = u + i 2 (k + 1 − 2j), i.e. the separation of adjacent constituent magnon rapidities is i, see Fig. 1a. Of course this distribution pattern of Bethe roots is not exact when the length L of the chain is finite and in fact large deviations are observed. Nevertheless the string hypothesis can be used to perform a counting of all states which gives the expected exact result even at finite L [2] (see also [6] for a recent account including references).
An interesting type of "Bethe string" which has not found much attention until recently is one where the rapidity u scales as the length L of the chain. Here one distinguishes between short and long strings where the number of constituent magnons is either of O(1) or of O(L). For these strings the regular pattern is violated strongly, the distance between adjacent Bethe roots deviates much from i. Long strings were first investigated by Sutherland in [7]. This particular string consists of a condensate core  and two tails, see Fig. 2b. Similar long strings were later considered in [8]. The distance of Bethe roots in the core is very close to i while in the curved tails the density of roots decreases to 0 at the ends. Short strings were considered in [9,10]. Unlike for standard strings the distance of Bethe roots is of O( √ L), see Fig. 1b. Short strings can as well be considered as very short versions of long strings [11]. The distinction of the two types is nevertheless useful because of their different energy scales in the limit of large L. Furthermore, long strings can be considered as smooth classical objects for which the discrete magnon constituents play a minor role. Short strings on the other hand are quantum objects and the number of constituent magnons must be a positive integer. In fact a short string is better viewed as a coherent state of bosonic magnons. In the strict thermodynamic limit the magnons do not interact with each other, their rapidities are merely influenced by the long strings which do interact non-trivially among themselves.
Interest in the ferromagnetic regime was recently sparked by the AdS/CFT correspondence. Minahan and Zarembo showed the equivalence [10] of the spectrum of planar one-loop anomalous dimensions in the su (2) sector of N = 4 supersymmetric gauge theory to the spectrum of the Heisenberg magnet with zero overall momentum. Consequently it was observed that the energies of certain long ultra-relativistic string configurations in R × S 3 [12] agree precisely with the energies of certain long strings of the Heisenberg ferromagnet in the thermodynamic limit [8], see [13] for reviews. This extended the earlier observation of [14] that the spectrum of magnons or short string excitations above the ferromagnetic vacuum agrees with the spectrum of quantum excitations of a short ultra-relativistic string orbiting the equator of S 3 . Some time later Kruczenski related both the ultra-relativistic limit of strings on R × S 3 and the thermodynamic limit of the Heisenberg ferromagnet to one and the same Landau-Lifshitz model on S 2 [15], see also [16].
The general solution for this classical model was constructed in [17,18] in terms of spectral curves, see [19]. It was furthermore shown that the moduli space of spectral curves of any given genus has the expected dimension. The features of the spectral curves are in one-to-one correspondence to the integral equations obtained from the above thermodynamic limit of the Bethe equations. It is therefore clear that any state of the Heisenberg ferromagnet in the thermodynamic limit is approximated by a spectral curve. Less obvious is the question if every admissible spectral curve also has a corresponding solution to the Bethe equations, and if not, what are the stability criteria? In [20] one such criterion was derived from the requirement that the mode number of a long string can be defined unambiguously and self-consistently. The statement is that the directed density ρ of Bethe roots should not encircle the points in, n ∈ Z \ {0} when moving along the contour of a string. This is for example achieved if the density is bounded by |ρ(x)| < 1 or | Im ρ(x)| < 1, i.e. the distribution of Bethe roots should not be denser than the ideal Bethe string. In particular we expect |ρ(x)| < 1 where the cut crosses the real axis. It is in fact a natural bound for the quantum model: The values ρ = in (along the imaginary axis) are distinguished because the pairwise interactions of Bethe roots in (1.1) become singular at that point. For the strictly classical model on the other hand the quantity ρ does not have a meaning, and solutions with density bounded by one are not at all different from solutions where the density exceeds one. A notion of stability exists in the semiclassical theory [21]; it demands the absence of tachyonic fluctuations. This condition does not coincide with |ρ| ≤ 1, it turns out that the latter is a stronger requirement in the case of a single cut [22].
So the question how to properly define stability remains. In this paper we would like to gain further insight into this issue by studying the general one-and two-cut solutions. The moduli space of two-cut solutions has two discrete parameters and two continuous ones, and it is sufficiently rich to investigate stability. In particular, we want to understand the role of condensate cores in the context of stability. In order to compare our analysis of the one-and two-cut spectral curves to actual solutions of the discrete Bethe equations (1.1), we further develop a numerical method for the construction of such solutions with large numbers of Bethe roots. Up to now, only few solutions that can be compared to corresponding spectral curves have been constructed [8], and they only have a comparatively small number of Bethe roots. The numerical solutions with finite numbers of Bethe roots represent quantum states on a chain of finite length and hence can also be used to examine finite-size corrections to the thermodynamic limit. For onecut configurations these corrections were studied analytically in [20,23,22], some more configurations of roots were analyzed analytically and numerically in [24]. In this work we examine numerically the leading order finite-size effects explicitly for configurations beyond the one-cut case.
This paper is organised as follows. We start in Sec. 2 with a review of the general solution in the thermodynamic limit by means of spectral curves. In order to illustrate our procedure on a simple example we will first reconsider and discuss the general one-cut solution in Sec. 3. The corresponding construction within the Landau-Lifshitz model is given in App. A. Sec. 4 contains the derivation of the general two-cut solution and its properties. In Sec. 5 we will continue by applying the two-cut solution to study the issue of stability. To that end we have to determine the physical shape of the branch cuts in various regions of the parameter space. Finally in Sec. 7 we test our predictions against solutions of the Bethe equations with large but finite length, constructed numerically, in order to substantiate our claims. We conclude in Sec. 8.

Spectral Curves for the Heisenberg Ferromagnet
We start by reviewing the spectral curves which describe solutions of the Bethe equations (1.1) in the thermodynamic limit [18].

Baxter Equation
To derive the properties of the spectral curves it seems convenient to consider the Baxter equation as advertised in [25] T (u)Q(u) = (u + i 2 ) L Q(u − i) + (u − i 2 ) L Q(u + i).
In terms of the Q-function the momentum and energy read (unless L < 2) .

(2.3)
Note that for gauge theory and AdS/CFT the states are required to be cyclic, i.e. the net momentum must be zero, e iP = 1. Here we will however not restrict to cyclic states but consider states with arbitrary momentum.

Thermodynamic Limit
Before we take the thermodynamic limit, we make the following substitutions: Define and t(x) = T (xL) (xL) L , (2.5) where p(x) is called the quasi-momentum and t(x) is the properly rescaled transfer matrix eigenvalue. It is clear that the function p(x) has logarithmic singularities with opposite prefactors at x = (u k ± i 2 )/L and a pole with residue i/2 at x = 0. Furthermore, we fix the ambiguity of p(x) by shifts of 2π by setting p(∞) = 0. Finally, t(x) is a polynomial in 1/x of degree L, i.e. it has an L-fold pole at x = 0 and is analytic everywhere else. In these variables the Baxter equation reads exp +ip(x − i/2L) It is now easy to take the thermodynamic limit L → ∞. In this limit the magnon number is assumed to scale like M = Lα with fixed total filling α. The magnon rapidities scale like u k ∼ L and are assumed to distribute smoothly along certain contours in the complex plane with distance O(L 0 ). The distance between adjacent magnon rapidities defines the density ρ(x) of Bethe roots Note that the contours typically arrange vertically in the complex plane (but not necessarily strictly along the imaginary direction). Therefore the density function ρ(x) is in general complex and its phase is inversely related to the direction of the contour at the point x. For definiteness, let us decide that cuts generally go upwards in the complex plane, Im u k+1 > Im u k . That means that the density will typically have a negative imaginary part, Im ρ(x) < 0. The Baxter equation in the thermodynamic limit becomes simply The logarithmic singularities in p(x) move together to form linear discontinuities. These discontinuities are interpreted as branch cuts connecting various Riemann sheets of p(x). Apart from these and the pole at x = 0 with residue i/2 the quasi-momentum p(x) is analytic everywhere. The function t(x) is also analytic except for an exponential singularity at x = 0 originating from the pole of degree L → ∞. Solutions p(x), t(x) to the equation (2.8) with the above properties describe the spectrum in the thermodynamic limit. Let us now study the properties of the solutions in more detail.

Branch Cuts
The function p(x) has discontinuities along certain branch cuts whose union shall be denoted by C.
Density. The discontinuity of a cut is proportional to the density ρ(x) of Bethe roots (2.7) in the vicinity of the point x We shall take ε to be an infinitesimally small positive number, i.e. p(x + ε) denotes the limiting value of p at x ∈ C towards the right of the cut. Note that the combination dx ρ(x) must be real and positive which determines the direction of physical branch cuts. The integrated density along a connected cut C k will be denoted by the filling α k Standard Cut. The function t(x) must remain analytic across a branch cut of p(x). Equation (2.8) tells us that p(x) can change sign and shift by a multiple of 2π without causing a discontinuity in t(x). If the sign changes across a branch cut we thus have where n k is the (constant) mode number associated to the connected branch cut C k . Using (2.9) we can relate the density to the quasi-momentum This type of cut can end in a square-root singularity x * of p(x) where it takes the value Note that this singularity of p(x) is indeed compatible with analyticity of t(x) at the singularity, cf. (2.8). At x * the branch cut can be oriented in three different directions: Consider the radial coordinates x = x * + re iϕ . The combination dxρ(x) ∼ e 3iϕ/2 √ rdr must be real and therefore the three possible orientations are separated by 240 • (where the full rotation is by 720 • due to the square root singularity).
Condensate Cut. Conversely, if the sign does not change across a branch cut, the quasi-momentum must obey This type of cut would be required to end on a logarithmic singularity which is not compatible with analyticity of t(x) in (2.8). Nevertheless such a can exist if it ends on other cuts: Indeed, we can view a logarithmic cut (2.14) as the union of two parallel standard cuts (2.11). Therefore a logarithmic cut C can split up into two parallel cuts C L (left) and C R (right) at some point, see Fig. 3. Evidence for this splitting is provided by some numerical solutions to the Bethe equations for small L [9,8]. Compatibility of (2.14) with (2.11) requires n = n L − n R . The integer n k determines the density of Bethe roots (2.9) to ρ(x) = −in k and thus a logarithmic cut has constant integral density and extends strictly along the imaginary direction.

Observables
According to (2.3) the momentum P and energyẼ = EL appear in the expansion of the quasi-momentum at x = 0 (2.15) Figure 3: Two standard cuts C L and C R with mode numbers n L and n R whose physical contours (as determined by the condition that ρ(x) dx be real) meet at a common point. The standard cut contours start and end at square-root singularities, which are indicated by crosses. Beyond their common point, the cut contours can be placed on top of each other -the result is a condensate cut as shown in the right picture. The density ρ c on the condensate cut can be determined by placing the individual contours infinitesimally close and parallel to each other, as shown in the left picture. It turns out that the density on the condensate cut is constant and purely imaginary: ρ c (x) = i(n R − n L ), see (2.14). Consequently, the contour of the condensate cut must be vertical.
The relation between t(x) and the momentum and energy is slightly trickier: The limit of the expansion (2.2) gives two asymptotic expansions at the exponential singularity of t(x) at where the terms O(u L ) in (2.2) are interpreted as exponential singularities. Putting the two expansions together and comparing with (2.8) is in agreement with (2.15). The total filling α = M/L is easily obtained from the definition (2.5) of the quasimomentum The expansion of t(x) around x = ∞ in (2.4) yields a compatible expression, which however does not fix the above sign (2.18) Figure 4: Cycles of the spectral curve dp. Shown are the two sheets dp 1 , dp 2 of the Riemann surface of the spectral curve dp. The two sheets are connected by cuts C i . To each cut C i is associated an A-cycle A i and a B-cycle B i . The A-cycle encircles the cut contour, staying on the same sheet, while the B-cycle starts at x = ∞ on one sheet, passes through the cut and ends at x = ∞ on the other sheet.
One can also derive two useful relations between the total and partial fillings 19) and the overall momentum and the mode numbers (2.20)

Spectral Curve
For the differential dp(x) the analytic structure is somewhat simpler: The shifts of p(x) by constants which arise when passing through a branch cut are not seen in dp(x). The differential dp(x) thus has merely two Riemann sheets with opposite signs and it is called the spectral curve. In particular, the condensate cuts drop out in the spectral curve. Only the standard branch cuts are seen in dp(x); they change the sheet or equivalently flip the sign. The above parameters of the solution can be read off from the spectral curve. First we need to introduce A-cycles and B-cycles for the cuts. A-cycles wind around the cuts while B-cycles extend from x = ∞ through a cut and back to x = ∞ but on the other Riemann sheet, see Fig. 4. 2 For a normal cut the A-cycle does not intersect with any branch cut of the quasi-momentum. The A-period of dp(x) is therefore zero A k dp(x) = 0. (2.21) The B-period yields the mode number n k of the cut B k dp(x) = 2πn k . (2.22) To understand this we split up the contour into the two parts before and after crossing the cut at point x. The first part yields p(x + ε) − p(∞) = p(x + ε) due to our normalisation p(∞) = 0. On the other sheet the sign of the quasi-momentum is flipped and consequently the second part yields −p(∞) + p(x − ε) = p(x − ε). Altogether we obtain p(x + ε) + p(x − ε) and by (2.11) this equals 2πn k . Note that a vanishing A-period of dp enables us to be rather unspecific about how the B-period returns to x = ∞ on the second sheet. Finally, the partial filling (2.10) of a cut can be expressed through the A-period of x dp(x) This follows by substitution of (2.9) and partial integration. Note that the partial integration assumes that the A-period of dp is zero. The above discussion shows that the cut contours and corresponding cycles can be deformed continuously without affecting the parameters of the solution. However, special care has to be taken when a cycle moves into a cut or passes through a singularity. In particular, the reality condition dρ(x) ≥ 0 appears to play no important role in the purely classical approximation.
In the case of many cuts the assignment of cycles is not necessarily unique anymore. Let us discuss what happens when two standard cuts join to form a condensate cut with four tails. Now there are various ways in which the four branch points could be connected, see Fig. 5. The standard method to set up the cuts is to connect complex conjugate pairs of branch points more or less directly, as in Fig. 5a. Then the above assignments of parameters works well. All cuts are of the standard kind and condensate cuts arise from two standard cuts running parallel for a while, see Fig. 3.
Another option would be to connect the branch points horizontally with a condensate cut forming between the two standard cuts (Fig. 5b). This option leads to non-zero Aperiods Figure 5: Possible cut contours (blue/red,thick) and corresponding cycles for a core with four tails. In the standard configuration (a), complex conjugate pairs of branch points are connected more or less directly. In this case, the A-periods vanish and the B-periods are non-ambiguous, even when the central cut segments join and form a condensate core as in Fig. 3. When the branch points are connected horizontally (b), the mode numbers (B-periods) of the cuts are ambiguous, hence the configuration requires a condensate (black,vertical) between the two cuts. This implies that the A-periods do not vanish. Branch cuts that wind around each other (c) also require a condensate cut, hence complex conjugate branch points carry different mode numbers.
where n j is the density of the condensate cut being intersected by the A-cycles. The nontrivial A-cycles then lead to ambiguities in the definition of mode numbers n k and fillings α k . Also non-technically it is unclear how to associate suitable n k and α k to these two individual cuts. The situation becomes even worse if the cuts wind around the branch points (Fig. 5c). Therefore the standard straight connection of complex conjugate pairs of branch points appears best to describe the parameters.

Finite Gap Solutions
The simplest types of spectral curves with the desired properties are the finite-gap solutions. These have finitely many branch cuts ("gaps") and therefore finitely many branch points. The general ansatz for p (x) is given by where g(x) and h(x) are polynomials of degree K and 2K, respectively The roots x k are the square root branch points which must come in (complex conjugate) pairs. A pair is typically connected by a branch cut and therefore K is the number of standard branch cuts. The genus of the algebraic curve p (x) equals g = K − 1.

Stability
Stability addresses the question which classical spectral curves with real and positive densities on the cuts can be realised approximately by solutions of the Bethe equations with large but finite L. An important stability criterion was derived in [20] The criterion implies that the branch of the above logarithm can be uniquely defined on an isolated standard branch cut C k . The argument of the logarithm approximates the scattering term in the Bethe equation (1.1). In a logarithmic form of the Bethe equations the condition implies that mode numbers can be unambiguously associated to the cuts, even in the full quantum theory. Unlike the classical quantities, here the density appears independently of the differential dx and therefore it is crucial to use the physical contour of the cut C k . The most interesting case is n = 1; the conditions for n > 1 appear to be less constraining and perhaps redundant. It is quite clear that if the absolute density is bounded by unity everywhere, |ρ(x)| < 1, (2.29) the stability criterion will be satisfied. However this condition is too restrictive in general. Nevertheless, for a stand-alone cut the density is typically highest at the centre x 0 where the cut crosses the real axis and where the density becomes purely imaginary. It is then necessary to obey |ρ(x 0 )| < 1, (2.30) in order to satisfy the stability condition. Otherwise the argument of the logarithm in (2.27) will cross the negative real axis at x 0 and wind around the origin once.
What remains obscure at this point is how to interpret the stability condition for two standard cuts joined by a condensate. Also it is not clear if (2.27) alone can ensure that a given configuration of cuts can be realised by a concrete solution of the Bethe equations. In the following we shall study the general one-cut and two-cut solutions in order to shed some light on the various classically allowed configurations of cuts and when the stability condition is satisfied.

One-Cut Solution
In this section we review the simplest type of spectral curve with one cut. It has two parameters, the mode number n and the filling α. Its genus is zero and thus we will only encounter algebraic and trigonometric functions. We also review the corresponding construction for the Landau-Lifshitz model in App. A.

Solution
The general one-cut solution was obtained in [18]. It has one mode number n, one filling α and it takes the form (3.1) The derivative of the quasi-momentum defines the spectral curve and matches with (2.25,2.26). The quasi-momentum is in agreement with the expansions at x = 0, ∞ (2.15,2.17) and with the condition for the branch points (2.13). The higher terms in the expansion at x = 0 yield the total momentum and energy If we choose to restrict to cyclic states as required for AdS/CFT we have to set P = 2πm with integer m. Then the total filling must be a rational number α = m/n. However, all these solutions are unstable. As we shall see later this is related to the fact that the total momentum leaves the first Brillouin zone, |P | > π.

Cut Contour
The physical contour of the cut is not a simple function. In particular, it is not given by the natural branch cut of the square root in p(x); it lies somewhat closer to the origin. In order to find the contour, we shall make use of the identity (2.12) which relates the density to the quasi-momentum. The integrated density must be real and positive; therefore we will need the integral of the quasi-momentum On the cut the combination Λ(x) − πnx must be real. At the branch points it takes the values Λ(x ± )−πnx ± = ± i 2 πα and thus the contour is defined by Λ(x)−πnx ∈ i 2 πα[−1, 1]. In Fig. 6 we have plotted the contour of a sample cut. In fact as discussed below (2.13), a branch cut with positive density can, in principle, originate from a branch point in three different directions. The shortest path will turn out to be the only physical choice.
The longer circular path which encircles the origin is equivalent to the one-cut solution with opposite mode number n = −n and conjugate filling α = 1−α. This can be seen as follows: We deform the shorter cut continuously to the longer cut by rotating it towards the right by 240 • . 3 At some point the cut has to pass the point x = ∞. This has the following two effects. Firstly, the A-cycle describing the filling (2.23) intersects x = ∞ twice. This adds to α the residues 1 2 − α and 1 2 − α (the opposite sheet with opposite circulation) so we obtain α = α + 2( 1 2 − α) = 1 − α. Secondly, the point x = ∞ now resides on a different sheet. This implies p(∞) = 2πn and we have to subtract the constant 2πn to normalise to p(∞) = 0. At the branch point that leads to p(x * ) = −πn, i.e. the new mode number is n = −n. Note that this configuration has a large filling α ≥ 1 2 because the shorter cut has α ≤ 1 2 . It is therefore unphysical. For the third choice we rotate the contour towards the left by 240 • . This contour escapes towards ±i∞ on both sides. Moreover the density ρ on the contour approaches a finite value at infinity. Therefore the total filling on the contour is infinite, and the configuration is inconsistent.

Stability
We are now ready to consider the issue of stability. The stability condition for the onecut solutions implies that the density must be bounded by unity where the cut crosses the real line. Assume the cut crosses the real axis at x 0 which must be the solution of the equation Λ(x 0 ) = πnx 0 . This equation is transcendental and we can only solve it numerically for given values of n and α. Once we have the value x 0 the absolute density is given by The density for solutions with n = 1 is plotted in Fig. 7. We are interested in the filling α cond where the density reaches the maximum allowed value, |ρ(x 0 )| = 1 or p(x 0 ) = π(n + 1). Let us assume for simplicity that n > 0. Because of the way p(x) and Λ(x) scale with x and n, the intersection point of a cut with mode number n and filling a and the density at that point are directly related to the intersection point and the density of the cut with mode number 1 and the same filling: x 0 (n, α) = x 0 (1, α) n , ρ n (x 0 (n, α)) = nρ 1 (x 0 (1, α)) . Hence it is sufficient to obtain the fillings α at which the absolute density at the centre of the cut with mode number 1 equals 1/n. These fillings are exactly the maximal fillings of the cuts with mode number n that are allowed by the stability criterion; they are indicated in Fig. 7.  Table 1: Maximal fillings α cond and the corresponding energiesẼ cond for onecut solutions for the first few mode numbers n. Above these fillings, solutions require a condensate for stability.
Alternatively, one can solve the equation p(x 0 ) = π(n + 1) for x 0 in terms of α and n. The solution is substituted in Λ(x 0 ) = πnx 0 which yields the equation Here q is an auxiliary variable that encodes the maximal filling α cond and the intersection point x 0 We list the first few values of α cond in Tab. 1. Note that α cond < 1/2n, which implies that stable solutions have a momentum P < π. Hence, cyclic one-cut solutions cannot be stable. In fact one can solve the above equation perturbatively using that q ∼ 1/n.

Fluctuations
Fluctuations are very small cuts which can exist in a background the one long cut [8]. For the one-cut solution they were discussed in [8,26,22]. Their position x * is determined by the long cut, p(x * ) = πn * , (3.11) but they are not strong enough to back-react on the position of the long cut substantially. We shall again assume that n > 0. The solution to the above equation reads The momentum and energy of a fluctuation mode quantum with δα = 1/L are given by [26,22] (3.13) If the long cut is rather short, the fluctuations will reside close to their vacuum positions x * ≈ 1/2πn * . For fluctuations of the same mode number as the cut, n * = n, the condition p(x * ) = πn * is solved by the branch points, which means that this particular fluctuation will increase the size of the long cut. The longer the long cut gets the more will it attract fluctuations with nearby mode numbers n * ≈ n towards it [22]. At some value of α, the fluctuation with n * = n + 1 collides with the long cut (from the left). This happens precisely when the density reaches unity at the real axis, cf. (3.6), and so the above discussion applies. Note that at this point the fluctuation with n * = n − 1 is still at some distance from the cut.
In conclusion, we can infer that for all stable one-cut solutions the fluctuations are well-separated from the long cut. When the fluctuation collides with the cut we can actually argue independently for an instability: Now the long cut can be filled not only from both ends, but also from the middle. In practice we expect that a condensate cut (with unit density) will form in the middle of the long cut.

Condensate Formation
Let us add a vertical condensate cut ending on the existing branch cut, see Fig. 8. This is achieved by shifting the quasi-momentum p(x) by −2π in the region enclosed by the two contours, cf. Fig. 9. Clearly the new curve satisfies all conditions for the classical spectral curve. In effect, the shift decreases the mode number on the inner part of the original branch cut by one unit according to (2.11). Furthermore, the density function is changed according to (2.9,2.12) and thus the inner part of the branch cut has to be moved to obtain a real density (3.14) The condensate cut can end on any point of the branch cut, and we show the various potential configurations in Fig. 8. Any of these configurations with a condensate appear to be possible from the point of view of spectral curves. However, the Bethe equations will single out one particular configuration as the distribution of Bethe roots in the thermodynamic limit. To understand the physical distribution we have to distinguish three qualitatively different cases, see Fig. 8.
In the first case the filling is below the threshold for condensate formation, α < α cond , as discussed in the previous subsections. Here all potential deformations leave the branch cut almost vertically and then circle around the origin. Due to the residue at the origin the filling of the configuration is altered and becomes larger than 1 2 . Inserting a condensate cut therefore leads to an unphysical configuration when α < α cond .
At α = α cond the fluctuation point with mode number n * = n + 1 crosses the branch cut and effectively acquires the mode number n * = n − 1. We now have two fluctuation points with the same mode number n − 1. For definiteness, let us call the new point n * = (n − 1) . When increasing α further, it will eventually collide with the other fluctuation point. This happens at α = α crit with   : Possible deformation of a branch cut: For large enough filling α > α cond , a closed loop cut (thick/green) can be added to the branch cut (thin/red). In the region enclosed by the loop cut, the quasi-momentum p(x) is changed to −p(x) + 2πn L , where n L is the mode number of the loop cut. As shown in Fig. 3, the density ρ C on the common part of the two contours equals i(n L − n B ), where n B is the mode number of the branch cut. Hence the resulting configuration has a condensate core, as shown in the right Figure. Consider now the case α cond < α < α crit . The deformations originating from the branch cut can now have two qualitatively different shapes. The contours originating from near the ends of the branch cut behave like in the first case above and are thus unphysical. The contours originating from near the centre however are just small deformations of the original cut and consequently they lead to the same filling. All of these configurations are in principle okay, however, only one can be physical. Indeed one of the configurations is special: the limiting shape, i.e. the largest possible small deformation. It is distinguished by a cusp at the fluctuation point for n * = (n − 1) , see Fig. 8b. We can argue that this is the physical configuration: When we add macroscopically many Bethe roots to the fluctuation point, it will split up and form two square root singularities. The deformed cut will split at the cusp and form a condensate with four tails as in Fig. 10b. This configuration is a genuine two-cut solution as we shall see in the next section. When we take the new Bethe roots away we should return to the one-cut solution. This is possible only if the deformed contour meets the fluctuation point with mode number n * = (n − 1) , as in Fig. 8b,10a. In conclusion the one-cut solution for α cond < α < α crit is a degenerate case of a two-cut solution.
The final case is α crit < α. Here the two fluctuation points with mode numbers n * = (n − 1) and n * = (n − 1) have joined and branched off into the complex plane. The configurations are similar to the above case of α cond < α < α crit , but now the limiting shape meets both fluctuation points and thus has two cusps joined by an approximately vertical contour. In this case we cannot add Bethe roots to any of the two fluctuation points individually because it would violate the reality condition (reflection symmetry of the configuration about the real axis). They can only be excited in pairs to create four new square root singularities and a three-cut solution as in Fig. 10d. In other words, the configuration is a degenerate case of a three-cut solution. The third cut is the line segment joining the two fluctuation points. By removing Bethe roots from this cut one lowers the energy and we might therefore call this solution unstable. Nevertheless we expect it to be a perfectly well-behaved solution of the Bethe equations in the thermodynamic limit albeit with non-minimal energy. A local minimum is obtained by shifting all Bethe roots from the third cut to the second, as in Fig. 10e. This is a true two-cut solution which will be discussed in generality in the next section.
Finally we remark that for n = 1 the third case does not exist because α crit = 1 2 . This is because the fluctuation point with n * = 0 is always at x * = ∞ and cannot join with the one for n * = 0 . Therefore the solution with n = 1 represents a local minimum of the energy for all 0 < α < 1/2.

The General Two-Cut Solution
As described in Sec. 2, a solution to the Bethe equations (1.1) in the thermodynamic limit is given by a quasi-momentum p, whose domain is a multi-sheeted cover ofC. The corresponding spectral curve dp has B-periods in 2πZ and vanishing A-periods, while the quasi-momentum p(x) has a simple pole of residue 1/2 at x = 0. In this section, the general spectral curve with two cuts will be constructed and its properties will be investigated. Spectral curves with two cuts have a Riemann surface of genus one and thus can be described in terms of complete elliptic integrals of the first, second and third kind. These are denoted by K, E and Π and are defined as The two-cut solution can be obtained by direct integration of the general ansatz (2.25). Here a different approach will be followed which makes use of the result obtained for the symmetric case in [27].

The Symmetric Two-Cut Solution
We will obtain the general two-cut solution by generalising the symmetric solution [28,27] p 0 (z) := − ∆n a 0 z where q = 1 − a 2 0 /b 2 0 . As a function of z, the elliptic integral of the third kind Π(z| q) has branch points at z = 1 and z = ∞. Its discontinuity across the branch cut between these two points (from the lower to the upper half plane) equals 2πi times the residue of the integrand in (4.1) at t = 1/ √ z: .
The functionΠ(z) has two branch cuts that lie symmetrically around z = 0. In principle, a 0 and b 0 can be arbitrary complex numbers. However, physical branch cuts must be symmetric about the real axis, therefore only the case b 0 =ā 0 is physical. The discontinuity across the cuts can be inferred from (4.3) and equals .

(4.4)
In (4.2), the factor in front ofΠ(z) makes the discontinuity constant, but also switches sign across the cut. As a result, the sum of the limiting values of p 0 on either side of a cut is This means that for ∆n ∈ 2Z, the Bethe equations (1.1) are satisfied and the cuts (a 0 , b 0 ) and (−a 0 , −b 0 ) have mode numbers ∆n/2 and −∆n/2 respectively. 4 SinceΠ(0) = K(q), p 0 (z) has a simple pole at z = 0 with residue For b 0 =ā 0 , q = 1 − exp(4i arg a 0 ) lies on a unit circle centred at 1 5 and hence depends only on the argument of a 0 , not on its modulus. As a result, one real degree of freedom remains after imposing the condition res 0 p 0 (z) = 1/2, which corresponds to the filling of the two symmetric cuts.

Construction of the General Two-Cut Solution
The solution (4.2) can be generalised to the non-symmetric case by forming the compo- It turns out that transforming the symmetric solution this way preserves enough of its structure while providing sufficient freedom for constructing arbitrary two-cut solutions. By suitably choosing the Möbius transformation µ and the modulus q, one can map any four points {a, b, c, d} ∈C onto the branch points {a 0 , b 0 , −a 0 , −b 0 }, i.e. one can construct a function with any four branch points {a, b, c, d} by solving the system of equations for a 0 , b 0 and the parameters of the transformation µ. In fact, the equations (4.8) do not completely fix a 0 , b 0 and µ. However, only solutions whose cuts are symmetric under reflection about the real axis are physical. Hence it is useful to set b 0 =ā 0 and restrict oneself to real Möbius transformations (i.e. s, t, u, r ∈ R), which map complex conjugate pairs to complex conjugate pairs. Applying the transformation to p 0 directly would move the pole from z = 0 to x = µ −1 (0) = u/t. In order to have the pole at x = 0 in the transformed function, it has to be moved to z = µ(0) = u/s in the original function p 0 . This can be achieved by adding a term to p 0 The prefactor is necessary for retaining the structure of cuts of the original function without introducing additional poles at z = ±a 0 . The function to be transformed now readsp , (4.10) and the candidate for the generalised solution is The constant C must be chosen such that p(∞) = 0, therefore it has to be defined as The function p has branch points {a, b, c, d} = µ −1 ({a 0 , b 0 , −a 0 , −b 0 }) and branch cuts C 1 := (a, b), C 2 := (c, d). As can be verified directly by computing the derivative of (4.11), dp is of the form (2.25). The following sections will show that the function p indeed meets all requirements for being a valid quasi-momentum while providing maximal freedom in the choice of physical parameters.

A-periods and Mode Numbers
For being a valid quasi-momentum, p must have vanishing A-periods and integral mode numbers (B-periods), as was established in Sec. 2.5. Because the function (4.11) is single-valued by construction, the integrals of dp over A-cycles vanish, as long as one chooses the cuts C 1 and C 2 of p between the branch points ( in a way that they do not cross each other. Since p vanishes at infinity, its mode numbers n 1 and n 2 are (cf. (2.22)) and similarly (4.14) For obtaining physical solutions, the parameters of the function p must be chosen such that n 1 and n 2 are integer numbers.

Energy and Momentum
The residue of the pole as well as the total energy and the total momentum of a solution p are encoded in the expansion of p(x) at x = 0. Performing this expansion, comparing with (2.15) and demanding that the residue of the pole be 1/2 yields the equations

Solving for Physical Parameters
The parameters of the solution (4.11) are a 0 , b 0 , ∆n and the three parameters of the Möbius transformation µ. The solution p is only physical if these parameters are adjusted such that all physicality conditions are satisfied. Namely, the mode numbers (4.13) and (4.14) must be integral, the filling fractions (4.19) must be real and (4.15) must hold. The correct behaviour at infinity (2.17) is provided by construction. Setting b 0 =ā 0 and restricting to real Möbius transformations results in complex conjugate branch points, which guarantees that the filling fractions are real. Further, equation (4.15) can be solved for t (or r, alternatively): Imposing these constraints leaves the complex a 0 , the real ∆n and two real parameters of µ as free parameters. Taking into account that the curve p is invariant under the rescaling {a 0 , t, u} → λ {a 0 , t, u}, λ ∈ R, only four real parameters remain. They correspond to the freedom in the choice of the mode numbers n 1 , n 2 and the filling fractions α 1 , α 2 . This is consistent with the result of [18], according to which the algebraic curve with k cuts has 2k parameters: k discrete mode numbers and k continuous fillings.
The expressions for the mode numbers (4.13) and (4.14) can be easily solved for ∆n, yielding ∆n = n 1 − n 2 . (4.21) The physicality constraints cannot be solved analytically for the remaining parameters, hence one has to resort to numerical methods for finding explicit solutions p for given mode numbers and filling fractions (or other sets of physical quantities, such as energỹ E and total filling α).

Finding the Cut Contours and the Density
For a given solution p, the physical contours of the two branch cuts are determined by the condition that ρ(x)dx be real. Together with the relation (2.12) between the density and the quasi-momentum, this condition yields a first-order differential equation for the cut contours. The physical cut contours can be obtained by numerically integrating this equation. Alternatively, one can obtain an explicit expression for the integrated density, similar to (3.4) in the one-cut case. The integral of the density reads where Λ(x) is the integral of the two-cut quasi-momentum p(x) and n j is the mode number of the respective cut. The integral Λ can be calculated analytically and is given in App. D. On the physical contour, the expression (4.22) must be real when x * is one of the branch points. Note that the function Λ(x) can also be used for obtaining initial positions of Bethe roots for the computation of discrete Bethe root distributions. Namely, when one wants to approximate a given two-cut spectral curve p with a solution to the discrete Bethe equations on a chain of length L, then this solution has to consist of M = αL roots. M j = α j L of these roots have mode number n j and must lie on or close to cut C j , j = 1, 2. Taking the solutions x k to the equations as initial positions for the Bethe roots reproduces the density function on the respective classical cut well and thus gives a good chance of numerically finding an approximating solution of the Bethe equations. Of course, the same method can be employed for the one-cut solution, using the respective integral (3.4). This way, the solutions shown in Fig. 1b,2 were obtained. In Sec. 7.1, we describe another simple way for obtaining discrete distributions of Bethe roots which become exact in the limit of large length L and fixed number of roots M . They can be easily constructed for any number of cuts and we use them to obtain numerical solutions of the Bethe equations in Sec. 7.2. However (4.23) gives a better approximation when L is not sufficiently large.

Moduli Space for Consecutive Mode Numbers and Stability
In the last section, the general two-cut spectral curve has been constructed. For given mode numbers n 1 and n 2 , the moduli space of this solution is two-dimensional. It can be parametrised by the two fillings α 1 and α 2 , by the total filling α and the energyẼ or by any other pair of physical quantities. In this section, the moduli space of configurations with two consecutive mode numbers is investigated. It turns out that this space has a very rich structure and even is globally connected, that is the spaces of different pairs of mode numbers are interconnected. Some features of the moduli space for non-consecutive mode numbers are presented in the following section.

General Features of the Moduli Space
As discussed in Sec. 3, a single cut with mode number n c and a small filling is centred around 1/(2πn c ). When the filling of the cut increases, the neighbouring excitation points with mode numbers n c ± 1 get attracted by the cut, until the excitation n c + 1 (for positive n c ) collides with the cut when the absolute density at its centre reaches unity. The same happens when the excitation n c + 1 has a small but finite filling, i.e. when it is replaced by a small cut. As will be shown in more detail below, larger cuts in general attract smaller cuts when their filling is increased.
To begin with, the moduli spaces of configurations with two consecutive mode numbers can be divided into several regions, as is shown in Fig. 11,12,13,14,15 for the cases of n 1 = 1, n 2 = 2 and n 1 = 2, n 2 = 3. Solutions in each region and on each line separating the regions share certain characteristics. In the following sections, the different regions, lines and points of the moduli space will be discussed in detail and examples will be given. For simplicity, only non-negative mode numbers n i ≥ 0 are considered, which correspond to cuts that lie to the right of the imaginary line. The discussion for negative mode numbers is completely analogous. The two cuts are named "cut one" and "cut two". Unless otherwise stated, the mode number n 1 of cut one is smaller than the one of cut two, n 2 = n 1 + 1.
In all figures that show cut contours or densities, cut one has blue colour (thick lines), while cut two has red colour (thin lines). Plots of cut contours in the complex plane always have an aspect ratio of 1 : 1. In all density plots, the density is shown versus the distance from the cut's centre, measured along the cut's contour.

Regions
The moduli space of configurations with two given consecutive mode numbers can be divided into three regions. Configurations in each region share certain characteristics: In region I, the two cuts are disjoint; in region II, the cuts join with each other and form a condensate with four tails; in region III, the two cuts are disjoint but require a third, closed cut and a condensate for stability. In the following paragraphs, the characteristics of the different regions are described in more detail and example solutions are shown.     Fig. 15. Structurally, the moduli spaces of higher pairs of mode numbers like this one are very similar to the case n 1 = 1, n 2 = 2 shown in Fig. 11. One can already see that the part of the space shown here is connected to the part shown in Fig. 11.
Region I n,n+1 : Two Disjoint Cuts. In this region, the contours of the two cuts with mode numbers n and n + 1 are disjoint, and the cut with mode number n + 1 lies further left (closer to the origin) than the one with mode number n. The absolute density is bounded above by 1 everywhere on both cuts. An example is shown in Fig. 16. When the filling of one of the cut increases, this cut attracts the other one, as can be seen in Fig. 17,18.
Region II n,n+1 : Condensate with Four Tails. In this region, the standard contours 6 of the two cuts bend towards and cross each other. An example is shown in Fig. 19. As explained below (2.14), the sum of the densities at the two common points of the standard contours equals −i (since the difference of the mode numbers is n 1 − n 2 = −1) and solutions in this region allow for straight condensate cuts with constant density ρ c (x) = −i between the two common points of the cuts. The condensate-cut version of the configuration in Fig. 19 is shown in Fig. 20. Beyond the two common points, the remaining parts of the standard cut contours form four tails that are attached to the ends of the condensate. The tails of cut two (with mode number n + 1) lie further left than the ones of cut one (with mode number n).
Whenever a particular solution p allows for a condensate cut, it will be assumed that the configuration with the condensate cut is the physical version of the solution. This assumption may not seem very motivated at the moment, but will become more justified during the subsequent discussion. A first sign of its verity is the fact that there is no well-defined direction on standard cut contours that cross each other. This can be seen as follows: According to (2.12), the density ρ on a standard cut contour is proportional to the difference between the limiting values of p(x) on either side of the cut contour: When a standard cut contour crosses another standard cut contour, the value of p(x ± ) on either side of the contour changes to −p(x ± ) + 2πn, where n is the mode number of the cut that is traversed; this follows from the Bethe equations (1.1). Hence the density ρ on the contour switches its sign, and in order that the corresponding root density ρ(x) dx stays positive, the direction on the contour must change as well. However a unique direction on each contour was assumed in the derivation of the spectral curve and quasi-momentum in Sec. 2.3. This argument makes it seem very unlikely that configurations with crossing standard cut contours can be realised as limits of Bethe root distributions.
As indicated by the dashed lines in Fig. 11,12,13,14,15, region II n,n+1 can be further subdivided. Between the lines D n,n+1 , G n,n+1 and C n , the standard cut contours have the usual form, as in Fig. 19: Neither do they encircle the origin, nor do they extend to infinity. In contrast, beyond line G n,n+1 the standard contour of cut one encircles the origin, as in the example in Fig. 21. Moving a contour with filling α 1 past the pole at the origin changes its filling to α 1 = α 1 + 1; this follows from the expression of the filling fraction as a contour integral (2.23), the residue theorem and the fact that the residue at x = 0 switches its sign as the contour moves past it. As a result, the absolute density on the standard contour typically exceeds unity at its centre, as in Fig. 21. This and  (b) Absolute densities along the cuts vs. distance from the cuts' centres (measured along the contour) Figure 16: Cut contours and absolute densities of an example solution in region I 12 (see Fig. 11). The fillings are α 1 = 0.016 (blue/thick, n 1 = 1) and α 2 = 0.03 (red/thin, n 2 = 2).    (b) Absolute densities along the cuts vs. distance from the cuts' centres (measured along the contour) Figure 19: Standard cut contours and absolute densities of an example configuration in region II 12 (see Fig. 11). The fillings fractions are α 1 = 0.11 (blue/thick, n 1 = 1) and α 2 = 0.05 (red/thin, n 2 = 2). This particular configuration lies below line G 12 in the (α,Ẽ) plane. The physical version of the configuration has a condensate and is shown in Fig. 20. : An example configuration in region II 12 beyond line G 12 (see Fig. 11). The filling fractions are α 1 = 0.032 (blue/thick, n 1 = 1) and α 2 = 0.07 (red/thin, n 2 = 2). The physical version of this solution has a condensate cut and is shown in Fig. 22.
the fact that the pole at the origin has the wrong sign in this configuration are strong indications that the physical realisation of such a solution is the one with a condensate cut, as in Fig. 22.
Physically, there is no qualitative difference between configurations on either side of line G n,n+1 , as their physical versions always contain a condensate cut and thus do not see the central part of the standard cut contours. Note that below line G n,n+1 , as in the example in Fig. 19, the absolute density does not exceed unity on the standard cut contours. Hence the stability criterion (2.29) is obeyed by the standard contours and no condensate is necessary to fulfil (2.27). Yet a condensate version of the configuration is possible. Beyond line G n,n+1 a condensate is necessary for stability, and for continuity reasons one would assume that also below line G n,n+1 the physical version is the one with a condensate cut. This is another sign in favour of the assumption that if a condensate cut is possible, it is also realised physically.
The other dashed lines H n,n+1 and J n,n+1 that further subdivide region II n,n+1 are discussed in Sec. 5.3 below.
Region III n,n+1 : Two Disjoint Cuts, One with a Condensate. In this region, the contours of the two cuts are disjoint as in region I n,n+1 , but the absolute density exceeds unity at the centre of cut two, i.e. the stability condition (2.29) is violated. An example configuration is shown in Fig. 23. Configurations in this region are still physical, because the same discussion as in Sec. 3.5 is applicable: A third, closed cut can be added in order to form a condensate at the centre of cut two. Since the absolute density at the centre of cut two exceeds unity, the excitation point with mode number n + 2 has   Figure 23: An example solution in region III 12 (see Fig. 11). The fillings are α 1 = 0.0038 (blue/thick, n 1 = 1) and α 2 = 0.07 (red/thin, n 2 = 2). The absolute density exceeds unity at the centre of cut two, hence the stability criterion (2.29) is violated. The physical version of the configuration involves a condensate cut and is shown in Fig. 24. already passed through the contour of cut two and changed its mode number to n. The closed cut originates and ends at this excitation point. The resulting physical version of such a configuration is shown in Fig. 24.

Lines
In the following paragraphs, the features of configurations on the various lines shown in Fig. 11 Line A n : Only One Cut. This line represents configurations where only the cut with mode number n is present and has a filling α < α cond (cf. Tab. 1), such that the absolute density at the cut's centre always stays below unity. Among all configurations with mode numbers n and n + 1 and given total filling, the state on this line is the one with the lowest energy. Adding a second cut with higher mode number while keeping the total filling constant obviously increases the energy. Conversely, adding a second cut with lower mode number while keeping the total filling constant decreases the total energy.
Line B n : Only One Cut with a Condensate. This line is the continuation of line A n . Solutions on this line have only one cut with mode number n which has a filling α cond < α < α crit and were discussed in Sec. 3. The absolute density at the centre of the cut's standard contour exceeds unity; hence a second, closed cut that starts and ends at the excitation (n − 1) and forms a condensate with the central part of the cut is required for stability, as in the second example in Fig. 8.
Line C n : Condensate with Two Tails. Line C n separates region II n,n+1 from region II n−1,n , which both consist of solutions with a condensate cut. Solutions on line C n consist of a larger cut (cut one) with mode number n and a smaller cut (cut two), whose end points lie exactly on the standard contour of the larger cut. Therefore, the physical version of such a configuration consists of a condensate with one tail at each of its ends; the smaller cut is completely hidden in the condensate. Since the density of cut two at the cut's end points vanishes, the density of cut one at these points must equal −i, as explained below (2.14). Consequently, also the direction of cut one at these points is purely imaginary, i.e. vertical. Hence, the cut contour of the condensate version of this configuration and the density along it are smooth. Solutions on this line were investigated by Sutherland in [7], refer to Figure 2 in that reference for a picture of the absolute density on the condensate and the tails. A series of configurations in which the branch points of cut two pass the contour of cut one is shown in Fig. 25. The sequence passes line C 1 from left to right in Fig. 11. During the series, the filling fraction of cut two stays constant, while the filling fraction of cut one increases. At the beginning, the configuration is in region II 12 : The branch points of cut two lie to the left of cut one, and the mode number of cut two is n 2 = 2, while the mode number of cut one is n 1 = 1 (Fig. 25a,25b,25c). While the filling fraction of cut one increases, cut two gets attracted by it and moves further right. During this process, the cuts increasingly bend towards each other and both become longer. At a certain value of the large cut's filling, the end points of the small cut lie exactly on the large cut's contour, the configuration has arrived at line C 1 (Fig. 25d). One can see that the contour of cut one is indeed vertical at the end points of cut two. The filling of cut two stays constant during the series, but its length increases, hence the absolute density along its contour decreases.
When the filling of cut one increases further, the small cut's branch points move on to its other side and the configuration reaches region II n−1,n (Fig. 25f,25g). While the end points of cut two pass the large cut, multiple quantities change: According to (4.13), the mode number of a cut is n i = πp(x * ), where x * is a branch point of the cut. Across a cut with mode number n 1 , p switches sign and shifts by 2πn 1 (2.11). Therefore, a cut with mode number n 2 whose branch points pass another cut with mode number n 1 changes its mode number to n 2 = −n 2 + 2n 1 . Specifically, if n 2 = n 1 + 1, as in the series in Fig. 25, then n 2 = n 1 − 1. Further, while cut two passes through the contour of cut one, the sign of the density (2.12) along cut two changes, as shown above in the discussion of region II n,n+1 . This does not affect the density of roots along the cut, but it switches the sign of the filling: α 2 → α 2 = −α 2 . Thirdly, the contour of cut one and its filling change by a finite amount when cut two passes. The filling changes from α 1 to α 1 = α 1 + 2α 2 , which is consistent with the fact that the total filling α 1 +α 2 , as a coefficient of the expansion of p at infinity (2.17), must vary continuously during the series. Also the total momentum m = n 1 α 1 + n 2 α 2 varies continuously during the transition from region II n,n+1 to region II n−1,n .
For a qualitative picture of the densities on the cuts in the series, refer to Fig. 19b,26b. As a consequence of the change in the contour and the filling of cut one, in all solutions beyond line C n (right of C n in the (α,Ẽ) plane) the absolute density exceeds unity at the centre of cut one's standard contour. The standard contour of cut two is disjoint from cut one, as in Fig. 26. In this form, the configuration does not allow for a condensate cut and is therefore unstable. However, there still is a stable version of such solutions: As discussed below (2.13), there are three possible starting directions for a cut at each branch point. Choosing a different direction for cut two, this cut's contour can be made to cross the contour of the first cut, and the unphysical density at the first cut's centre becomes hidden in a condensate. For the solution in Fig. 26, his type of configuration is shown in Fig. 27. These condensate-cut versions of the solutions beyond line C n not only cure the instability at the centre of cut one, they also allow for a smooth transition from configurations on line C n to solutions beyond line C n . This can be seen in Fig. 28, where the physical version of the series of Fig. 25 is displayed.
The changes that occur when the branch points of cut two pass the contour of cut one imply that configurations on line C n have an ambiguity in the mode number of cut two and in the standard contour of cut one. This can be seen in the example in Fig. 25d,25e. In the limit where cut two approaches cut one from the left, the mode number of cut two is n 2 = n 1 +1, and the contour of cut one is bent to the left. In the other limit, where cut two approaches cut one from the right, the mode number of cut two is n 2 = n 1 − 1, and cut one is bent to the right. The fillings of the cuts behave accordingly, as described in the discussion of the series above. On line C n , the contour of cut one is unique between (b) Absolute densities along the cuts vs. distance from the cuts' centres (measured along the contour) Figure 26: Standard cut contours and absolute densities of an example solution in region II n−1,n , beyond line C 1 (see Fig. 11). The fillings are α 1 = 0.25 (blue/thick, n 1 = 1) and α 2 = 0.04 (red/thin, n 2 = 0). The absolute density exceeds unity at the centre of cut one, hence this configuration is unstable. The physical, stable form of this solution is shown in Fig. 27.  Fig. 26. Cut two (red/thin) has a different starting direction than in Fig. 26 and hence a different contour, which allows for the formation of a condensate cut that hides the central part of cut one (blue/thick) and renders the configuration stable. its end points and the branch points of cut two. The fact that the contour is ambiguous between the two branch points of cut two is consistent with the fact that the expansion of p has a square root term at the branch points x * (cf. (2.13)): because this results in the expansion of the density (as follows from (2.12)) which implies that the solution of the differential equation ρ 1 (x)dx ∈ R has an ambiguous second derivative at the branch points x ± and hence there are two choices for the continuation of the contour beyond the branch points x * . However, this ambiguity has no physical significance, as the physical configuration with the condensate is non-ambiguous.
Line D n,n+1 : Two Tangential Cuts. This line separates region I n,n+1 from region II n,n+1 . The two cut contours of configurations on this line bend towards each other and touch each other in one point on the real axis. At the common point of the contours, the sum of the two densities equals −i. An example is shown in Fig. 29. Beyond this line, the formation of condensate cuts begins.
Line E n,n+1 : Cut Two Violating Stability. This line can be reached from region I n,n+1 by increasing the filling of cut two while keeping the filling of cut one small, such that the two cuts stay disjoint. Line E n,n+1 is reached when the absolute density at the centre of cut one reaches unity. At this point, the excitation with mode number n + 2 collides with the contour of cut two from the left. Further increasing the filling of cut two makes the excitation pass to the right of the contour and change its mode number to n. The absolute density on the standard contour of cut two now exceeds unity and a third, closed cut is required for stability; the configuration has reached region III n,n+1 (cf. Fig. 24).
Line F n,n+1 : Two Disjoint Cuts and a Cusp on Cut One. This line separates region III n,n+1 , where the standard cut contours of the two cuts are disjoint, from solutions in region II n,n+1 , in which cut one encircles the origin and consequently crosses the contour of cut two. A sequence of configurations that passes line F 12 from region III 12 to region II 12 is shown in Fig. 30. The transition happens when the excitation n , which carries a closed loop that forms a condensate with cut two, reaches the contour of cut one. At this point, the contour of the closed loop merges with cut one, whose standard contour encircles the origin afterwards. The transition point is marked by line F n,n+1 . On this line, the density (2.9) on cut one decreases to zero at the merging point x 0 , as p(x 0 ) = πn. Around this point, the density expands to ρ 1 (x 0 + ε)dε ∼ ε 2 , which implies that the contour has a cusp with opening angle 90 • .
Line G n,n+1 : Condensate with Tails and a Virtual Cusp. This line in region II n,n+1 separates solutions in which the standard contour of cut one closes on the right side of the origin from solutions in which it encircles the origin (cf. Fig. 19,21). The qualitative change that happens to cut one's contour when crossing this line is the same as on line F n,n+1 : The excitation point n collides with the standard contour, making the density decay to zero at its centre. Unlike on line F n,n+1 , here the process is less (c) 0.0705 Figure 30: A sequence of configurations that passes line F 12 . Cut one (blue/thick) has mode number n 1 = 1 and constant filling α 1 = 0.005, cut two (red/thin) has mode number n 2 = 2 and its filling α 2 is indicated below each subfigure. The excitation with mode number 1 is shown and carries a closed cut (green/thick) which forms a condensate cut at the centre of cut two. Configuration (a) lies in region III 12 , configuration (b) lies on line F 12 and configuration (c) lies in region II 12 . The dashed contours all close around the origin. The energy increases during the sequence.
significant because the central part of cut one is hidden in the condensate. An example configuration that lies on line G 12 is displayed in Fig. 31.
Line H n : One Cut with Vanishing Virtual Filling. On this line, the filling α 1 of cut one (with mode number n − 1) vanishes. However, as shown above in the discussion of region II n,n+1 , the standard contour of cut one encircles the origin and the filling on this contour equals one. Since most of cut one is hidden in the condensate cut, the vanishing filling of the physical cut contour does not have any further implications for the physical configuration. An example is given in Fig. 32. Because the filling of cut one vanishes, there is a one-cut solution with the same mode numbers and fillings for each solution on this line. These one-cut solutions are of the type shown in Fig. 8c and have a higher total energy than the corresponding solutions on line H n . Refer to the discussion of point Y n below for more details on the relation between these two types of solutions.
Line J n : Excitation n Collides with the Condensate. These lines mark all configurations where the excitation point with mode number n is located on the condensate cut. When line J n gets traversed, the excitation point with mode number n collides with the condensate from the left and passes to the right (for positive mode numbers), thereby changing its mode number to n − 2. But unlike an excitation point passing a standard   The fillings are α 1 = 0 (blue/thick, n 1 = 1) and α 2 = 0.085 (red/thin, n 2 = 2). The parts of the standard contours that are hidden in the condensate cut are shown as dashed lines. Since the standard contour of cut one encircles the origin, the integral of the absolute density along it is not zero, but unity.
cut, this process does not create an additional, closed cut (as on line E n,n+1 ); in fact, the passage of the excitation point has no effect on the configuration at all.

Special Points
As a result of the discussion of the various regions and lines in the previous sections, three types of special points emerge, labelled by X n , Y n and Z n,n+1 in the Fig. 11,12,13,14,15. These points will be discussed in the following.
Point X n : Excitation Point Meets Single Cut. The points X n mark the one-cut solutions with mode number n whose absolute density equals unity at the cut's centre, which means that the excitation point with mode number n + 1 is located on the cut's contour. The filling fractions of these solutions for n ≤ 10 can be found in Tab. 1. The point X n connects line A n , on which the absolute density on the single cut's centre falls below unity, with line B n , on which the absolute density at the centre of the standard contour exceeds unity and a condensate begins to form with the help of a second, closed cut. Also the lines C n , D n,n+1 and E n−1,n begin at the point X n (the last one only if n > 1). Among all configurations on line A n , the one at this point has the highest energy (cf. Tab. 1).

Point Y n : The Two Excitation Points with the Same Mode Number Collide.
Point Y n marks the end of line B n and thus represents a one-cut solution with a condensate. On line B n , the two excitation points (n−1) and (n−1) both lie on the real line, as in Fig. 8b. At the point Y n , these two excitation points collide. Analytically continuing line B n beyond point Y n lets the branch points diverge into the complex plane and leads to one-cut solutions of the type shown in Fig. 8c. This analytic continuation of line B n shall be called line B n . As discussed in Sec. 3.5, solutions on this line are limiting cases of three-cut solutions. Hence line B n is not part of the moduli space of two-cut solutions (or of its boundary). The only way to continue beyond point Y n while staying in the regime of consecutive mode numbers is to enter either region II n−1,n or region III n−1,n . In both cases, the excitation points (n − 1) and (n − 1) separate again. The former stays on the real line, while the latter splits into two square root branch points that diverge into the complex plane. Upon entering region III n−1,n , the excitation point (n − 1) continues to carry a closed cut as on line B n while the new branch points form the ends of a second standard cut, as in Fig. 30a. When the configuration enters region II n−1,n , the afore closed cut opens up at its cusp and becomes a standard cut that ends on the new pair of branch points, as in Fig. 30c. Either way, a second standard cut emerges, which can have a positive, negative or vanishing filling. Configurations with vanishing filling belong to line H n and lie in region II n−1,n . Since the closed cut in the solutions on line B n also has zero filling, line H n seems to be the most natural continuation of line B n .
As noted at the end of Sec. 3, solutions on line B n can be transformed into the corresponding solutions on line H n by separating the two segments of the closed cut in Fig. 8c. In doing so, the two excitation points in the complex plane split into two square-root branch points each, which results in a three cut solution. The filling of the third, vertical cut can then be decreased to zero while the filling of the second cut increases to zero, such that the total filling fraction α stays constant. During this process, the energy decreases. Comparing the energy of solutions on line H n to the energy of the corresponding one-cut solutions on line B n (Fig. 33) shows that a second order phase transition occurs as a configuration follows line B n and continues on line H n . For example, the deviation of the energy on line H n from the energy on line B n has the following expansion in the total filling at point Y n : The prefactor of the quadratic deviation, including its dependency on n, has been guessed on the basis of numerical data which, for the first few values of n, matches this form with high accuracy. On a technical level, this phase transition is related to the one discussed by Douglas and Kazakov in [29]: It involves a transition between the same classes of algebraic curves and it is also related to a density threshold, see App. B for a technical review. The authors study the large-N limit of the partition function of continuous pure Yang-Mills theory on the two-sphere and find a phase transition with respect to the area of the sphere. 7 Point Z n,n+1 : Cusp on Cut One Meets Cut Two. At this point, the lines F n,n+1 and G n,n+1 meet: Cut one has a cusp at the centre of its contour whose tip lies on the contour of cut two. Since the tip of the cusp marks the excitation point with mode number n + 2, this implies that the absolute density of cut two at this point equals unity and that lines E n,n+1 and J n+2 also end on point Z n,n+1 . The two cuts touch each other only in one point, thus also line D n,n+1 ends here.

Global Structure of the Moduli Space
As a consequence of the discussion about the lines C n in Sec. 5.3, all moduli spaces of two consecutive mode numbers are analytically connected through the lines C n . Pictures of the resulting global space are shown in Fig. 34,35. In the latter, also the lines of unstable one-cut solutions B n are shown for comparison, while in the former these lines are the same as the lines H n . In particular, the global connectedness of the space shows that mode numbers lose their intuitive meaning at large filling fractions: A small change in the number of excitations (filling fractions) can alter the mode numbers, hence they can no longer be interpreted as the number of periods of a coherent state on the spin chain. Note that the case n 1 = 0 (cut one), n 2 = 1 (cut two) qualitatively differs from the other cases with consecutive mode numbers. The excitation point n * = 0 always lies at x = ∞ and hence cannot be excited, only the excitation point n * = 0 can. Consequently, line B 1 cannot be crossed, and cut one always has a negative filling and its standard contour winds around the origin. That means that solutions with mode numbers n 1 = 0, n 2 = 1 always have a condensate and regions I 01 and III 01 do not exist.
The one-cut lines A n and B n , n > 1, form branch cuts of the global moduli space; across them, the derivatives of physical quantities such as the energy and the momentum are discontinuous. The points Y n , n > 1 are the corresponding branch points.
As mentioned at the end of Sec. 3.5, at point Y 1 the excitation point n * = 0 reaches x = ∞. At that point, the filling fraction of cut two equals α 2 = 1/2 and cut one extends along the whole imaginary line while the branch points of cut two are x ± = ±i/2π. The absolute density of this solution is shown in Fig. 36. When the filling of cut one gets reduced to finite negative values while α 2 stays at 1/2, the condensate and the tails of the configuration remain on the imaginary line, but their length decreases. Solutions of this type constitute line C 0 ; an example is shown in Fig. 37. All solutions on this line have a momentum of P = 2π(n 1 α 1 + n 2 α 2 ) = π.
As mentioned, a consequence of the fact that the phase space of consecutive mode numbers is connected is that the mode numbers can change by a continuous variation of the fillings. In particular, if one follows a straight horizontal line in Fig. 34 from the left to the right, the branch points cycle around each other while the mode numbers decrease until the configuration reaches line C 0 and both pairs of branch points and the contour lie on the imaginary line. A schematic sketch of this process is shown in Fig. 38. The sequence can even be continued further by following it in reverse, but in the other (left) half of the complex plane: The mode numbers become negative, but all other physical quantities are the same as in the mirror solution with positive mode numbers. This means that there exist two copies of the space in Fig. 34, one for positive and one for negative mode numbers, and the two parts are connected through line C 0 .     Also the lines B n , n > 0 of unstable one-cut solutions (as shown in Fig. 8c) can be followed until the branch points reach the imaginary axis, similar to the sequence of two-cut solutions that approaches line C 0 . At this point, also the fluctuation points n * = (n − 1) and n * = (n − 1) lie on the imaginary line and the filling of the branch cut reaches α = 1/2, as follows from (3.5). But, according to (3.3), this implies that P = πn > π. Since the maximal physical momentum on a lattice is P = π, the physical momentum is P − 2π for P > π. The maximal physical momentum P = π is reached by the configurations on line B n long before the branch points arrive at the imaginary axis. Reaching the regime of such high momenta is prevented by the phase transition discussed at the end of Sec. 5.4, which makes the one-cut solution transform into a two-cut solution at the critical value of the filling α crit = α Yn . Physically, nothing peculiar would happen when the momentum would increase to values larger than π, but it is interesting to see that the two-cut solutions by themselves implement the restriction to momenta |P | ≤ π, as can be seen in the pictures of the phase space in Fig. 34,35.

Stability and Gauge Theory States
The expositions of this section show that the moduli space of two-cut solutions has a very rich structure, even though only solutions with consecutive mode numbers have been considered. Most remarkably, there seem to be no unstable two-cut solutions: For all spectral curves with consecutive mode numbers there is a configuration of standard and condensate cuts that generate the curve and that do not violate the stability condition. Here, by stability condition is meant that the absolute density never exceeds unity at any point on the cut contour whose tangent is vertical: A remark about suitable gauge theory states is appropriate at this point. Only spin chain states with a momentum P ∈ 2πZ correspond to gauge theory states. None of the solutions studied in this chapter satisfies this condition: All two-cut configurations with consecutive mode numbers have a non-zero momentum 0 < P ≤ π. Nevertheless, all qualitative results should generalise to valid gauge theory states. Namely, configurations whose cut structure is point symmetric about the origin manifestly have a vanishing momentum P . 8 It is reasonable to expect that such point symmetric solutions with mode numbers (n, n + 1) and (−n, −n − 1) and with symmetric fillings exhibit the same qualitative features as the solutions investigated here (such as condensate formation and connectedness of the moduli space). Some point-symmetric configurations of two-cut and degenerate multi-cut solutions are studied in the next section.

Moduli Space of Non-Consecutive Mode Numbers
All qualitative features of the general moduli space with arbitrary mode numbers can already be seen in the case of consecutive mode numbers, which was studied in the previous section. In general, cuts with non-consecutive mode numbers lie further away from each other and hence influence each other less. Some examples will be studied in this section.
6.1 Mode Numbers n 1 = 1, n 2 = 3 First, the comparatively simple space of solutions with mode numbers n 1 = 1 (cut one), n 2 = 3 (cut two) shall be exposed. It is shown in the (α, P/2πα) plane in Fig. 39. In region I, the fillings of both cuts are small, their contours are disjoint and their absolute densities are bound by unity, hence there are no condensates. Line D 1 marks the configurations in which the density on cut one equals −i. Beyond this line lies region II 1 , whose configurations have a condensate on cut one and an associated closed loop-cut with mode number 0 = 2 . Similarly, configurations on line D 2 have an absolute density of −i at the centre of cut two; beyond this line lies region II 3 , where cut two has a condensate core and an associated closed loop-cut with mode number 2 = 4 . In region II 13 , both cuts have a condensate core and associated closed loop-cuts. Finally, line C 24 marks configurations in which the excitation points with mode numbers n = 2 and n = 4 have collided. Configurations beyond this line are unstable in the sense that they are degenerate cases of solutions with higher genus (higher number of cuts).
Regions II 1 and II 13 can be further subdivided, which is signified by the dashed and dotted lines in Fig. 39. To the left of these lines (at low total filling α), the two cuts with their condensates and loop-cuts are disjoint. On the green, dashed line, the excitation point with mode number 2 = 4 (which is the cusp of the loop-cut of cut two) collides with the condensate core of cut one. Upon further increasing the total filling, the condensate core of cut two (or its standard contour) collides with the condensate core of cut one; this event is marked by the orange, dotted line. At even higher total filling, the branch points of cut two pass through the condensate core of cut one (orange, dashed line), such that the whole contour of cut two (with its associated loop-cut) is encircled by the loop-cut of cut one. Nothing special happens to the condensate core of cut one during the passage of cut two -this is sensible, as the condensate core is invisible in the spectral curve dp, which contains all the information about the solution. When cut two has completely passed through the condensate core of cut one, its mode number has changed to n 2 = n 2 − 2 = 1.
When the total filling gets increased up to α = 1/2 (line E −1 ), the excitation point with mode number 2 = 0 , which marks the cusp of the closed loop-cut that is connected to cut one, reaches x = ∞. This results in a configuration similar to the one at the point Y 1 (see Fig. 36): The contour of the closed cut passes through the branch points of cut one, such that cut one is completely hidden in a condensate that has two vertical, infinitely long tails. Consequently, the mode number of cut one is ambiguous: It is either n 1 = 1 or n 1 = −1. In addition, there is cut two, whose mode number is n 2 = 1. It has a standard contour (below line D 3 ), or a standard contour with a condensate core and a closed loop-cut (above line D 3 ). At the upper end of line E −1 lies point V, at which also the excitation point that marks the cusp of cut two's closed loop-cut has reached x = ∞ and the configuration is completely symmetric about the origin (Fig. 40a). Consequently, the momentum vanishes at this point, P = 2π, and the solution represents a valid gaugetheory state. Fig. 40b shows another cut configuration that realises the solution at point V. It can be constructed by taking another possible standard contour for each cut and replacing their common centre with a double condensate, that is a condensate with root density |ρ| = 2. Bethe root distributions that show this type of pattern have been found before [9,8]. Behind point V lies line E 1 . Configurations on this line are "mirror solutions" of the ones on line E −1 : Cut three has mode number 1/ − 1 and consists of a condensate with two vertical tails that extend to infinity, while cut one has mode number −1 and consists of a standard contour with or without a condensate core and an associated closed loop-cut.  Fig. 39. This solution is symmetric about the origin and hence represents a valid gauge-theory state (The momentum equals P = 2π). Figure (a) shows the configuration that is reached when coming from region II 13 or from line E −1 in Fig. 39. There is a dual configuration of cuts, shown in (b), that generates the same solution. It has a "double-condensate" with root density |ρ| = 2 at it's centre.

Mode Numbers
The moduli space of solutions with mode numbers n 1 = 1 (cut one), n 2 = −1 (cut two) is shown in Fig. 41,42 in the (α, P/2πα) and in the (α, E) plane. Regions, lines and points with the same labels but subscripts 1 and −1 fall together in Fig. 42, since their configurations only differ by a reflection about the imaginary line and hence have the same energy.
As in previous examples, the lines A ±1 and B ±1 represent solutions with a single cut with mode number n, solutions on the lines A ±1 have a standard contour without a condensate core, while solutions on the lines B ±1 have a condensate core and a closed loop-cut. In the regions I ±1 and II ±1 , both cuts with mode numbers n 1 = 1 and n 2 = −1 are present. Line D ±1 separates solutions that only have a standard contour (region I ±1 ) from solutions that have a condensate core with a closed loop-cut (region II ±1 ). The closed loop-cut has mode number n = 0 for both cuts. Hence on line D ±1 , the absolute density at the centre of the cut with mode number ±1 equals |ρ| = 1.
In the regions I 1 and II 1 , the filling fraction of cut one is larger than the filling fraction of cut two, α 1 > α 2 . In the regions I −1 and II −1 , it is the other way round, α 1 < α 2 . These regions are separated by the lines S and S c , which consist of symmetric solutions with α 1 = α 2 . Such solutions have vanishing momentum P and hence represent valid gauge theory states, they were studied in [27]. Obviously, solutions on line S have plain standard contours, while solutions on line S c have condensate cores in both contours. Point X s marks the symmetric solution with maximal filling that does not have condensates. The lines E ±1 and point V are the same as the ones discussed above in Sec. 6.1, hence the   Fig. 39,41 are connected. The points Y ±1 actually both represent the same solution, namely the one shown in Fig. 36 whose cut configuration lies entirely on the imaginary line.

Numerical Solution of the Bethe Equation
In this section we develop a method for obtaining numerically exact distributions of Bethe roots for the Heisenberg spin chain. The main goal of this section is to test our predictions from the analysis of the one-cut and two-cut spectral curves in Sec. 3,5,6. Thus we are especially interested in distributions which can be compared with the spectral curve, i.e. solutions with large numbers of roots that align on a set of contours.
We examine in detail some subtle points of the previous analysis, such as stability (Sec. 2.7), condensate formation (Sec. 3.5) and phase transitions between spectral curves with different topologies (Sec. 5.4). We will see that condensate cores with roots separated by ∆u ≈ i indeed appear in all cases considered. We will also confirm that the phase transition from one-cut to two-cut solutions indeed happens. Furthermore, in Sec. 7.3 we will analyse small deviations from the thermodynamic limit. In particular we will numerically analyse the leading finite-size corrections and show that they do not exhibit any non-analyticity due to the formation of the condensate.
The numerical configurations shown in this section are approximate solutions to the discrete Bethe equations (1.1) with finite length L. It is astonishing to observe how well the predictions from the thermodynamic limit describe the actual distributions of Bethe roots.

Description of the Numerical Method
Here we describe the numerical method used to solve the Bethe equations (1.1) in detail. In order to find roots of a set of equations numerically, it is important to specify a good set of starting points that are sufficiently close to the exact solution. We will focus on this problem first and then describe a simple programme for finding numerical solutions.
Starting Points. In Sec. 2.3 we introduced the mode number associated to a branch cut in the thermodynamic limit. We will need to generalise the notion of mode numbers to the quantum case. The mode numbers naturally appear in the quantum theory when one takes the logarithm of (1.1): The integers n k represent the 2πi-ambiguity of the logarithm, in the thermodynamic limit they exactly become the mode numbers introduced in equation (2.11). This can be seen by taking the large-L limit of (7.1), as in [7,8,18]. The solution to (7.1) is uniquely specified by a set of integers n k (at least for large L and small M ). In order to find the solution numerically by Newton's method, one has to specify starting points for the search. For generic parameters it is very difficult to give suitable starting points that lead to convergent iterations. However, as we will see it is pretty simple to find good starting points in the limit of small filling fractions L M where the Bethe equations (7.1) can be solved approximately [9]: To the leading order we can drop the r.h.s. of (7.1) and, assuming that u k ∼ L (which is true in the case we are interested in), we get that u k L/2πn k . We conclude that Bethe roots with the same mode number n k are close to each other and separated by ∆u ∼ L from other groups of roots with different mode numbers. For resolving the positions of the Bethe roots with the same mode number n, we can neglect terms in the r.h.s. of (7.1) containing roots with different mode numbers and use the ansatz Assuming that n ∼ 1 and z k ∼ 1, we expand (7.1) and get [9] 1 , and one clearly sees that the relative mismatch is indeed smaller for larger L. It is possible to continue the 1/ √ L expansion in a systematic way to get a better description of Bethe roots for large L and mode number n fixed. For our goal the leading approximation of the Bethe roots by zeros of Hermite polynomials provides us with a good set of starting points for sufficiently large L.
Iteration Procedure. Having starting points at hand it is almost immediate to construct a numerical solution of the Bethe equations using Newton's method. However, we know the starting points only for large L. To overcome this difficulty we can view the Bethe root distributions as functions of the length L while keeping the magnon number M fixed. We can solve the Bethe equations numerically, gradually varying L from a sufficiently large value L M to the desired value with L = O(M ). In each step we use the numerical solution obtained in the previous step as new starting values for smaller L. In this way we obtain a series of configurations whose filling α = M/L slowly increases.
Some important examples of such series will be considered in the next subsection.
For various reasons it is favourable to use the Bethe equations in product form (1.1) and not in logarithmic form (7.1) for the Newton iteration. Using this method, the mode numbers enter into the calculation only through the starting points. Since we vary L slowly, at each step we inherit the information about the mode numbers from the previous step, and hence the resulting configuration corresponds to the analytically continued curve from the regime of large L.
Note that for large L and fixed M the roots are separated by ∆u ∼ √ L and thus the stability condition is satisfied, |ρ| ≈ 1/∆u 1. However, when we decrease L, at some point the separation between Bethe roots approaches ∆u = i. The Bethe equations (1.1) are singular at this point which leads to a numerical instability at small L. 9 Therefore the step size for L must be chosen sufficiently small in this regime. For practical purposes one can even choose L to take real intermediate values. We expect the above numerical procedure to lead to a reasonably unambiguous final configuration if L is changed in sufficiently small steps.
Sample Programme. Let us give a basic realisation of above method in Mathematica. For simplicity the programme presented in Tab. 2 is restricted to the one-cut case. To generalise it, it is enough to change the set of starting points u0 according to the discussion in the previous section. We tested this programme on Mathematica versions 5.2 and 6.0, but it should be compatible with earlier versions as well.
This programme can be optimised to work much faster, especially for the case of a large number of Bethe roots: A major improvement can be achieved by replacing the built-in FindRoot function, which is too general and inefficient, by a simple implementation of Newton's method. Furthermore the step size for L can be made adaptive. One can even decrease L in small steps during each Newton iteration and thus avoid an excessively precise calculation which is merely used for new starting values.  Table 2: Mathematica programme for the generation of a finite-length numerical Bethe root distribution corresponding to a one-cut spectral curve.

Some Series of Numerical Solutions
In this section we will exploit numerical solutions of the Bethe equations obtained with our programme to test the predictions from the thermodynamic limit. We begin with a study of the simplest one-cut solution, subsequently proceed to more complex configurations with two cuts and consecutive mode numbers, and finish with a two-cut solution with non-consecutive mode numbers. In the last configuration, one of the two cuts transforms into a two-cut complex in a phase transition, the result being a three-cut solution.
One-Cut Solution, n = 1. In the thermodynamic limit this solution is discussed in Sec. 3. An interesting phenomenon which we expect to arise for this configuration is the appearance of a closed loop-cut for filling fractions above α cond (see Fig. 8b).
The spectral curve of the one-cut solution can be discussed very explicitly, in particular the quasi-momentum is given by a simple algebraic formula (3.1) and one can find the shape of the physical branch cut quite easily. In the first figure of Fig. 44 the contour is shown together with the numerical solution for M = 60 roots and length L = 2925. We almost do not see any deviation of the roots from the asymptotical curve. On the next picture we considerably increased the filling fraction to M/L = 60/315 ≈ 0.19, and one can see that the distance between roots close to the real axis is just slightly above i. This means that the filling α is very close to α cond and the stability condition is almost violated. According to the discussion in Sec. 3.5, this implies that the fluctuation with mode number 2 is very close to the cut. In the figure, we indicate it by a red diamond.
When we continue increasing the filling fraction to values above α cond , the fluctuation point passes through the cut contour and produces a loop cut. At the intersection point of the loop cut with the original standard cut we expect the condensate to be formed. As is demonstrated in Fig. 44, this is indeed the case.
We can almost reach the maximal filling α = 1 2 numerically. As discussed in Sec. 5.5, this would correspond to a configuration with a condensate centred at the origin together with two tails that extend along the imaginary axis towards the fluctuation point, which has diverged to infinity. Finally, we notice that the roots belonging to the two tails on the right of the condensate do not seem to align on the contour predicted by the spectral curve (red/dashed). In fact this is an artefact of the finiteness of the number of roots. One can observe that the situation becomes better when the number of roots M is increased (while the filling fraction M/L is held fixed). The finite-size effects become very relevant at the tip of the loop because the density drops to zero at the fluctuation point forming the tip. Close to a standard branch point the roots align on the contour more accurately, because the density decreases only as the square root of the distance from the branch point, whereas at the cusp of the loop the density is proportional to the distance itself and hence the deviation from the spectral curve is more severe.
One-Cut Solution, n = 2. For larger mode numbers, the one-cut solution behaves similarly to the n = 1 case when the filling fraction is increased. However, we expect that a configuration with n > 1 exhibits a new type of behaviour when its filling fraction reaches α crit given by (3.15). As discussed in Sec. 3.5, we expect that a phase transition from one cut to two cuts occurs.
In Fig. 45, we compare the distributions of numerical Bethe roots with mode number n = 2 for different filling fractions. We see that the condensate forms when the nearest fluctuation from the left crosses the cut contour at α = α cond , which is in full agreement with the prediction of Tab. 1. In contrast to the n = 1 case, there is also a fluctuation point to the right of the cut (with mode number n = 1), even at small fillings α. When   we increase α, the fluctuation point with mode number n = 3, which initially was on the left of the standard cut, passes through the cut and approaches the fluctuation point with mode number n = 1 (second figure in Fig. 45). Beyond the point of collision there are two possibilities -the configuration could continue as a one-cut solution (with a loop cut with two cusps that could naturally be decomposed into two standard cuts), or it could transform into a two-cut solution with lower energy as explained in Sec. 3.5. We found that in our study always the second possibility is realised, our numerical method prefers the two-cut configuration with lower energy. It seems that this configuration is numerically more stable. In particular we see that in the last two pictures of Fig. 45 the roots perfectly follow the two standard contours of the two-cut spectral curve. Accordingly, Fig. 46 shows that the energy of the numerical configuration after the point of collision matches the energy of the two-cut configuration (up to deviations due to finite-size effects), which clearly indicates that the phase transition is indeed taking place.
Two-Cut Solution, n = 1, 2. The spectral curve for two-cut solutions was constructed in Sec. 4 and is less simple than the one-cut curve. The moduli space of two-cut configurations is also much more complicated, for the case of consecutive mode numbers it is described in detail in Sec. 5. Numerically, however, it is not difficult to construct a solution with an arbitrary number of cuts. We mainly want to show here that, in accordance with the prediction of Sec. 5, an intersection of the two cuts leads to the formation of a condensate that connects the two intersection points. In Fig. 47 we compare numerical two-cut solutions with mode numbers n = 1, 2 with the predictions from the spectral curve. In the first two pictures the two cuts are disjoint and hence we are in region I 12 of the phase diagram shown in Fig. 11. The last three configurations belong to region II 12 , and we see that indeed a condensate with roots separated by approximately i connects the intersection points of the two contours. On the third picture, the right (blue) n = 1 cut does not encircle the origin of the complex plane, but it has a cusp at the point of intersection with the real axis. This means that this configuration lies exactly on the line G 12 . In the last configuration shown in Fig. 47, the branch points of the n = 2 cut (red) are almost touching the n = 1 cut (blue) which means that the configuration is very close to the line C 1 .
Two-Cut Solution, n = 1, 3. Two-cut solutions with non-consecutive mode numbers allow for even richer properties. As is explained in Sec. 6.1, for some values of filling fractions both cuts can develop condensate cores before intersecting with each other. At higher fillings the cut with the higher mode number, which initially lies further left, can pass through the condensate core of the other cut. As explained in Sec. 6.1, this should not affect the inner structure of the cuts.
From the distributions of numerical roots for different values of the fillings in Fig. 48, we see indeed that the group of n = 3 Bethe roots (left, red) reaches the condensate core of the n = 1 roots (right, blue) and passes through the condensate. 10 These numerical solutions confirm that the passage of a cut complex through a condensate core does not affect the structure of the individual complexes, i.e. their shape and their root densities.
We did not plot cut contours for all configurations in this case, because similarly as in the previous example, the configuration undergoes a phase transition. Namely, this series of configurations follows the line P/2πα = 1.5, which is slightly above the critical value ∼ 1.4 so that we are crossing the line C 24 of the Fig. 39. When the two fluctuations denoted by the orange diamonds (n = 2, 4) collide, the corresponding spectral curve transforms into a three-cut curve, as is expected from the discussion in Sec. 6.1. Note that we have no means to plot the exact contours for three-cut configurations and thus these lines are missing from the figures.
Conclusions. We have presented numerical solutions for various types of configurations discussed in this paper, and in all cases we observe a perfect match with the analytical predictions of the spectral curve. This confirms the validity of the analysis of the thermodynamic limit in the previous sections.

One-Cut Solution and Finite-Size Effects
In this paper we mainly focus on the strict thermodynamic limit. However, the numerical solutions containing a finite number of Bethe roots allow us to easily analyse the leading correction to the thermodynamic limit. In string theory they correspond to the semiclassical energy shift, see [30,31] for a general discussion of canonical quantisation for the spectral curve.
The energy shift was studied analytically for one-cut configurations of the Heisenberg spin chain in [20,23,22]. The result reads Note that this result equals the sum of zero-mode energies of the fluctuation modes in (3.13) upon zeta-function regularisation of the sum. The result may in principle include exponentially suppressed contributions, see [32]. It is however conceivable as well that these corrections are absent in the Heisenberg/Landau-Lifshitz spectral curve which is considerably simpler than the spectral curve for the principle chiral model/string theory. The above energy shift was obtained assuming that the roots are distributed along a single cut and it is a priori questionable if this result remains true for the case α > α cond when the loop cut appears and the distribution of roots is more sophisticated. We can however use the intuition gained studying the strict thermodynamic limit: We saw that the energy is described by the same analytical expression (3.3) before and after the condensate is formed. Therefore we expect (7.4) to be valid also in the full range of α's.
Here we perform a numerical test of this natural guess. In Tab. 3 we give a series of deviations of numerical values of the energy from their asymptotical values (3.3) for different numbers of roots and for two fixed filling fractions below and above the point of condensation α cond ≈ 0.21. To compute δE ∞ we fit the values of the energies from the table for finite M by a constant and inverse powers of M . The value of the constant is precisely δE ∞ given in the table. For both filling fractions we see a perfect match with the value given by (7.4), which shows that the correction (7.4) is also valid for fillings beyond α cond .
A possible explanation for this could be based on the algebraic curve. The algebraic curve appearing in the thermodynamic limit can be modified to reproduce also the leading  Table 3: Deviations δE M of the numerical energy from its asymptotical value for different numbers M of Bethe roots and for two fixed filling fractions α above and below α cond ≈ 0.21. δE ∞ is the extrapolation to an infinite number of roots. The value E (1) (7.4) is given for comparison. quantum corrections [22,24]. The same type of analysis we have performed for the leading order hopefully can be applied for this corrected curve as well.

Conclusions and Outlook
In this article we have considered the general solution (Sec. 2) for the spectrum of the Heisenberg ferromagnet in the thermodynamic limit by means of spectral curves. Equivalently, we have studied the semiclassical Landau-Lifshitz model on S 2 , see App. A.
In particular, we have constructed the general two-cut (elliptic) solution (Sec. 4) and studied its properties (Sec. 5,6). This allowed us to shed light on the issue of quantum stability of classical solutions (Sec. 3,5) brought up in [20,22]. We have also tested our claims by constructing some numerical solutions (Sec. 7). Let us summarise our findings briefly: The two-cut solution is specified by two mode numbers n 1,2 and two fillings α 1,2 . Its quasi-momentum is an elliptic function which encodes all the relevant classical observables. In order to quantise this solution, one needs to understand the positions of the two branch cuts which correspond to the distribution of Bethe roots. We have studied these in detail and we have found a rich "phase space" (Fig. 34,39) of different configurations depending on the fillings. There are two main configurations, see Fig. 49: When the fillings are sufficiently small, there will be two separate branch cuts, see Fig. 49a. When the fillings are increased, the branch cuts bend towards each other so that they eventually would overlap. Instead of overlapping, a straight vertical condensate of Bethe roots with constant density forms. The configuration looks like a Bethe string but with four curved tails of lower density, see Fig. 49b.
Classically there is no quantitative difference between the two types of solutions, but the condensate serves several purposes for the quantum model: It sets an upper bound to the density of Bethe roots which is an important feature for quantum stability. In some sense, the condensate "hides" some potentially dangerous distributions of Bethe roots on an unphysical sheet of the spectral curve. Moreover, solutions with different mode numbers (but equal n 1 − n 2 ) can be connected analytically via the condensate. For example, there is only a single connected phase space for solutions with consecutive mode numbers (n 2 = n 1 + 1). Finally, we have observed that the naturally continued overall momentum P never crosses the boundary of the first Brillouin zone; it is apparently bounded by |P | ≤ π and the condensate formation plays an important role. At P = ±π the two-cut solution (with consecutive mode numbers) flattens out onto the imaginary axis where it connects smoothly the two regimes with P > 0 and P < 0.
Concerning quantum stability we find that an instability does not actually put a threat to the quantisation of some classical solution; the instability is rather caused by an unsuitable way of looking at the classical solution. This can be observed conveniently Figure 51: The configurations in Fig. 50b,50c should be understood as degenerate cases of two-and three-cut solutions. This can be seen by the addition of small but macroscopic numbers of Bethe roots to the fluctuation points, which "regularises" the configurations. The resulting configurations in (b), (c) are perfectly well-behaved two-and three-cut solutions. Exciting the fluctuation points of the solution in Fig. 50c lowers the energy, which is consistent with the (only seemingly) unstable excitation modes.
by means of the one-cut solution (Fig. 50) Here there are three main classes of cut configurations (which serve as subclasses for the two-cut solution): For small filling, there is just a simple cut, see Fig. 50a, in which case the solution is perfectly stable. For intermediate filling the naive one-cut solution violates the stability criterion but without any observable flaws like unstable fluctuation modes. In this case, the physical cut configuration develops a condensate and a loop with a cusp, see Fig. 50b, which achieves to bypass the stability criterion; it is nevertheless semiclassically equivalent to the one-cut solution. For large filling, the naive one-cut solution acquires an unstable fluctuation mode which appears to render the solution unphysical. The physical cut configuration now has a condensate and a loop with two cusps, see Fig. 50c, which gives a clue of how to understand the instability: The latter two dressed configurations should be interpreted as singular limits of twocut or three-cut solutions, respectively, see Fig. 51. In the naive one-cut configuration, a fluctuation would break apart the loop at a cusp by a tiny amount. This is fine when the cusp is on the real axis (for intermediate filling), but away from the real axis (for large filling) it leads to complex momenta and thus an instability. Conversely, fluctuations of the regularised three-cut configuration would merely change the filling of either the twocut complex or the third cut; in either case the configuration remains symmetrical with respect to the real axis, and there is no physical instability. Therefore the classical onecut state with large filling does have a corresponding quantum state. This state however does not minimise the energy locally in accordance with the apparent instability. The energy is lowered by taking away Bethe roots from the third cut and moving them onto the two-cut complex. At the end of this process the third cut has evaporated and only the two-cut complex remains which is a genuine two-cut solution (Fig. 52c,10e). This leads to the conclusion that the one-cut solution undergoes a phase transition at the point which separates intermediate and large filling: When the one-cut solution is analytically continued from small to large filling one obtains a degenerate three-cut solution. However in the vicinity of the latter there exists a two-cut solution with slightly smaller energy. It is useful to consider this solution the proper physical continuation of the one-cut curve. We find that at the phase transition point the second derivative of the energy is discontinuous. An exact expression for the discontinuity has been found on the basis of numerical data; an analytic confirmation of its form would be desirable. Technically, the phase transition turns out to be equivalent to the one observed in 2D QCD on S 2 [29] or equivalently a particular matrix model. In App. B we summarise the relevant functions and illustrate the phase transition for that model. Note that the apparent instability and the phase transition are merely artefacts of the thermodynamic limit in accordance with general results of statistical physics; they are resolved naturally at finite size, i.e. by the quantum theory: At finite L there exists no precise notion of continuous cuts and the distribution of Bethe roots onto these cuts is often ambiguous. Thus two cuts never actually end on the very same point which was an essential assumption for the instability and the phase transition. In particular we show that numerically exact solutions with finite L are smoothly interpolating between two phases and there is no real instability when the condensate forms.
We hope to have presented a fairly complete picture of the Heisenberg ferromagnet in the thermodynamic limit, its quantum fluctuations and the issue of stability. However there are several further questions which may be addressed in future work. One of them would be how precisely the third cut interacts with the two-cut complex in solutions of the type shown in Fig. 50c,51c. An investigation of this issue would require an explicit expression for the (hyperelliptic) three-cut quasi-momentum. Further, it would be interesting to construct the general elliptic solution to the Landau-Lifshitz model which includes the solutions obtained in [15,33]. It should be equivalent to the two-cut solution presented here and the circular solution (A.8) should transform into the elliptic solution in the same way the two-cut solutions emerges from the one-cut solution at large fillings. Another question would be how the thermodynamic solutions of the Heisenberg magnet behave under the influence of a magnetic field and how this would change the picture of the phase space.
Concerning the AdS/CFT correspondence, a generalisation of the results to the SU(2) principal chiral model describing the dual string theory on R × S 3 would be desirable. As discussed in [18], the solution to this model is given by a spectral curve that reduces to the Heisenberg spectral curve in a certain limit. Because of their strong structural similarity, we expect that stability and loop-formation work in qualitatively the same fashion. The phase space of general two-cut solutions (or elliptic string states) is likely to be similar as well, but new features may appear due to the more elaborate singularity structure on the curve. Therefore it is worth obtaining and discussing the general twocut curve for the principal chiral model. Some special cases of two-cut solutions have appeared in [27,34]. Analogously one could find the corresponding elliptic solution of the equations of motion explicitly, e.g. using the results of [30]. This solution may be related to the solutions studied in [35].
Finally, it remains an open question whether a similar study of the thermodynamic phase space can be carried out for spin chains which have different symmetry groups, quantum-deformed symmetry, or for which the spins transform in different representations. For example, studying the spectrum for symmetry groups of higher rank requires working with spectral curves with more than two sheets, which complicates the analysis substantially. Moreover for more general symmetry groups the algebraic curves are quite sensitive to different boundary conditions of the spin chain [24], thus making the analysis even more complicated (and interesting). cular solution and its fluctuation spectrum which corresponds to the one-cut solutions discussed in Sec. 3.

A.1 Definitions
The Landau-Lifshitz model on S 2 is a classical two-dimensional sigma model with two fields θ(σ, τ ) and φ(σ, τ ) which are the standard angular coordinates on S 2 . We identify the site k of the corresponding Heisenberg chain with the coordinate σ = 2πk/L of the sigma model. The expectation value of the spin projection S z at that site corresponds to the value 1 2 cos θ at σ. In analogy to the closed spin chain we take the world sheet to be closed, i.e. the coordinate σ is identified with σ + 2π: Consequently, the embedding coordinates θ and φ (more precisely the equivalence classes describing points on S 2 ) must be 2π-periodic in σ. This leads to the following periodicity conditions where the constants L/4π and π/2L are adjusted to match with our conventions for the Heisenberg chain. In fact the first term in the Lagrangian originates from a topological Wess-Zumino term in three dimensions and for a fully consistent treatment one should make the replacementφ cos θ → The equations of motion following from the Lagrangian arė Taking the derivative of L with respect toθ andφ yields the conjugate momentum densities The Lagrangian is homogeneous and therefore the total energy and total momentum are conserved. The rescaled total energy and a suitably shifted total momentum read E = EL = L π θθ + π φφ − L dσ = π 2 θ 2 + φ 2 sin 2 θ dσ, P = 1 2 φ dσ + 2π L π θ θ + π φ φ dσ = 1 2 1 − cos θ φ dσ. (A.6) Furthermore the Lagrangian is invariant under shifts φ → φ+ε which yields the conserved total spin S = L( 1 2 − α) with the filling α = 1 2 + 1 L π φ dσ = 1 4π 1 − cos θ dσ. (A.7)

A.3 Fluctuations
We would like to obtain the spectrum of fluctuations around this solution. For the case of α = 1 2 this problem was solved in [37] which coincides with the ultra-relativistic limit of spinning strings discussed in [21]. To this end, fluctuation modes with Fourier mode k and frequency ω k are added to the background solution with the following ansatz δθ = a √ L cos kσ + ω k τ + c , δφ = b √ L sin kσ + ω k τ + d .
Now it turns out that the total energy, total momentum and filling are not yet completely fixed. They still depend on the value of a and curiously also on the second variation of the σ-independent mode 11 δ 2 θ = f L .
(A.14) 11 This feature is related to the fact that this mode is already macroscopically filled in the classical solution. Essentially f corresponds to changing α in the classical solution by a microscopic amount.
We shall parametrise our ignorance by two constants N k and N 0 as follows This leads to the following shift in the observables The excitations enumerated by N 0 and N k carry each one quantum 1/L of filling equivalent to one unit of spin. They have mode numbers n and n + k, respectively. 12 Their energies can be read off from the above expression. It is not hard to see that a quantum of mode n essentially just increases α by 1/L in the classical solution, cf. (A.9). The charges for mode n + k are in complete agreement with the fluctuations in the ferromagnetic Heisenberg magnet (3.13). A careful semiclassical analysis actually shows that the numbers N 0 and N k must be integers. 13 The classical solution (A.8) is unstable if some of the oscillator frequencies ω k are complex. Obviously, this is the case if 2n α(1 − α) > 1 . (A.17) The transition from stable to unstable classical solutions exactly happens at the critical density α crit (3.15), where the one-cut solution of the Heisenberg magnet transforms into a two-cut solution, as discussed in Sec. 5.4.

B The Douglas-Kazakov Transition
Here we shall review the Douglas-Kazakov transition [29] for QCD on S 2 on a technical level. It mimics the transition between our one-cut and two-cut solutions, albeit in a somewhat simplified and more transparent fashion.

B.1 Definitions
The description of the partition function in the large-N limit of QCD on S 2 uses a density function ρ(h) with support on the real axis. It must obey certain integral equations which will not be needed in detail here. In order to extract the relevant observables it is convenient to introduce the continuation u(h) of ρ(h) into the complex plane. The function u(h) has a branch cut on the real axis such that the real/principal part equals the density The function u(h) is the analog of our quasi-momentum p(x) and its derivative u (h) is an algebraic curve. Apart from the branch cut, the only singularity is at h = ∞ with This expansion also defines the area A and the derivative of the free energy F . It is clear that these quantities play similar roles as the total filling α, momentum P and energỹ E.

B.2 Rational Regime
The simplest solution has one branch cut at [−h max , +h max ]. The function u(h) consequently is of semi-circle form with h max = 2/ √ A, see Fig. 53a. The resulting free energy from (B.2) is simply [38]  This function is physically correct for A ≈ 0. However, it was noticed that it does not have the correct behaviour at A → ∞. For instance F → − 1 24 but not F → 0.

B.3 Elliptic Regime
On the technical level, the pathology can be attributed to a density which grows beyond a threshold of ρ max = 1. (B.5) A density with ρ(h) > 1 is unphysical. Indeed for A > A crit with A crit = π 2 (B.6) the maximum density ρ(0) = √ A/π exceeds the maximum allowed density ρ(0) > ρ max . This problem is resolved by inserting a condensate with constant density ρ(h) = 1 on the interval [−b, b] [29]. The density then falls off to zero on the intervals [−a, −b] and [b, a]. This leads to the two-cut solution, cf. Fig. 53b u The expansion (B.2) at h → ∞ leads to in agreement with physical expectations [39].

B.4 Transition
We can now compare the two regimes. The expansion of the derivative of the free energy around the critical point A = A crit = π 2 yields in the two cases implying a third-order transition in the free energy F . In fact the function u(h) deviates only quadratically at a generic point away from the branch points, see Fig. 54b. It can be expected that the same holds whenever an elliptic curve degenerated into a rational curve. This leads to the conclusion that all quantities which can be read off directly from the curve (e.g. in an expansion) will deviate at least quadratically.

C Integration of the Partial Filling Fractions
Analytically integrating the partial filling fractions α 1 , α 2 of the two-cut solution (4.11) is possible because they can be expressed as integrals over A-cycles (2.23) of the derivativẽ p 0 (x), which is an algebraic function: dx . (C.1) In the following, this integral will be evaluated for the case 0 < a 0 < b 0 < ∞. The result can then be analytically continued to the case of complex a 0 and b 0 =ā 0 . For 0 < a 0 < b 0 < ∞ the contours A 1,2 encircle the real intervals I 1 := [a 0 , b 0 ], I 2 := [−b 0 , −a 0 ]. Except for an overall sign, the integrandp 0 (x)µ −1 (x) has the same values above and below the intervals I 1,2 . Therefore, the integrals along the contours A 1,2 can be expressed as integrals over I 1,2 : the square root in the integrand on the right hand side becomes The partial filling fraction integrals can now be written as where α 1 has the positive and α 2 the negative sign. The constants are A = u 2 K(q) − b 2 0 s 2 E(q) rs , A t/r = (a 2 0 rs − tu)(b 2 0 rs − tu) K(q) − b 2 0 (st − ru) 2 E(q) r 2 (st − ru) , A u/s = −(a 2 0 s 2 − u 2 )(b 2 0 s 2 − u 2 ) K(q) s 2 (st − ru) . (C.6) The integral (C.5) can be performed term by term. The first term is the simplest: The second and third term are of the same form, but the results for either contour differ by a sign: For obtaining these two forms, the integrals b 0 c 2Π (c) (C.10) were used. The latter can be found in [40], formula 217.12. Putting it all together and writing the occurring square roots in the form (4.17) in order to select the branch of the square root consistently, one arrives at the result (4.19).

D Integral of the Two-Cut Quasi-Momentum
In order to find the physical cut contours for a given two-cut quasi-momentum p(x), one can use the fact that expression (4.22) must be real on the physical contours. That expression contains the integral of the quasi-momentum, which will be given here. Direct integration yields · F arcsin (a 0 − b 0 )(a 0 + y) (a 0 + b 0 )(a 0 − y) and where F(z| q) is the incomplete elliptic integral of the first kind: Using the integral (D.1), one can obtain a closed expression for the integral Λ(x) of the two-cut quasi-momentum p(x) (4.11), which can be expressed as