Semi-classical quantisation of magnetic solitons in the anisotropic Heisenberg quantum chain

Using the algebro-geometric approach, we study the structure of semi-classical eigenstates in a weakly-anisotropic quantum Heisenberg spin chain. We outline how classical nonlinear spin waves governed by the anisotropic Landau-Lifshitz equation arise as coherent macroscopic low-energy fluctuations of the ferromagnetic grounds state. Special emphasis is devoted to the simplest types of solutions, describing precessional motion and elliptic magnetisation waves. To resolve their internal magnon structure, we carry out the semi-classical quantisation of classical spin waves using the Riemann-Hilbert problem approach. We describe overlaps of semi-classical eigenstates employing functional methods and correlation functions of these states with aid of classical phase-space averaging.


Introduction
Describing thermodynamic systems of interacting quantum particles represents one of the most challenging theoretical problems in physics. In spite of inherent complexity in describing the quantum many-body dynamics, there has been a tremendous progress in recent years in understanding its various facets such as thermalisation [1,2], far-from-equilibrium dynamics [3,4], quantum transport [5,6] and entanglement entropies [7]. In this respect, one-dimensional models play a special role, providing a playground for developing and testing analytical approaches and permitting efficient numerical simulations. While non-perturbative closed-form solutions are rarely accessible, there are prominent exceptions such as integrable models. In recent years, the focus of study has largely shifted towards various non-equilibrium phenomena and quantum transport, particularly after the inception of generalised hydrodynamics [8,9]. This paper is devoted to properties of semi-classical eigenstates and emergence of classical dynamics in an integrable quantum lattice model. And the phenomenon remains one of the most fundamental conceptual aspects of quantum many-body systems that has received very little attention thus far. Recently the interest is renewed in the the study of the N = 4 SUSY Yang-Mills theory [10,11], due to the emergence of integrable structures [12,13]. More specifically, semi-classical limit of isotropic Heisenberg model has been scrutinised in Ref. [10,11]. In this work we shall examine the semi-classical spectrum of the anisotropic Heisenberg spin-1/2 chain in the easy-axis regime, partly following the lines of Ref. [11].
The other source of motivation for studying semi-classical eigenstates and the emergent classical magnetisation dynamics comes from recent studies of spin transport in the anisotropic Heisenberg spin-1/2 chain. There are indeed many conspicuous qualitative similarities with spin dynamics in the classical continuous Landau-Lifshitz ferromagnet [14], which points to a particular type of a classical-quantum correspondence [15][16][17]. One piece of evidence for it comes from the exact solution of the classical initial value problem for the initial domainwall configuration, obtained in Ref. [16]. Depending on the value of anisotropy, three types of dynamical regimes are: (i) ballistic spin transport in the easy-plane regime, (ii) absence of transport in the easy-axis regime and (iii) diffusion with a multiplicative logarithmic correction at the isotropic point, independently confirmed to present also in quantum isotropic Heisenberg spin-1/2 chain [15]. Absence of transport in the easy-axis regime can be attributed to the presence of a kink [16]. Hoping to gain better insight into this curious correspondence, we set out to perform the semi-classical quantisation of classical nonlinear spin waves that arise in the weakly-anisotropic ferromagnetic spin chain.
Remarkably, the classical-quantum correspondence of spin transport indeed continues to hold even in thermal equilibrium at finite density. Most prominently, in the isotropic Heisenberg quantum chain, finite-temperature spin transport (in the zero magnetisation sections) is superdiffusive [18][19][20], belonging to the KPZ universality class with characteristic dynamical exponent z = 3/2 [21,22]. The anomalous nature of spin transport can been attributed to thermally dressed interacting 'giant magnons' -these refers to the semi-classical eigenstates residing in the low-energy spectrum of the model which at the classical level manifest themselves as soliton modes (as opposed to genuine quantum magnons are bounds states thereof (Bethe string) are completely suppressed in this regime) [23]. The anomaly nevertheless again disappears upon introducing interaction anisotropy: in the easy-plane regime on has ballistic transport characterized by a finite spin Drude weight [24,25], whereas in the easy-axis regime normal diffusion gets restored [21,26,27]. In our view, that these recent insights call for a careful and systematic study of semi-classical eigenstates in integrable quantum chains.
In this work we employ the Asymptotic Bethe Ansatz approach [28] to extract the semiclassical spectrum of excitations in the Heisenberg spin-1/2 chain in the presence of the uniaxial interaction anisotropy, following the lines of Refs. [10,11]. We proceed to describe the associated Riemann-Hilbert problem and outline the algebro-geometric technique for its solutions [29,30]. This provides complete information about the classical finite-gap spectrum of the anisotropic Landau-Lifshitz ferromagnet. Our attention will be largely devoted to the emergence of special types of semi-classical solutions which physically represent bions and kink configurations, which crucially require anisotropy for their stability. In addition, we shall provide closed expression for the semi-classical norms and overlaps, building on earlier works [31][32][33][34]. Finally, we briefly discuss the structure of the semi-classical limits of static correlation functions.
The article is structured as follows. In Sec. 2 we start with a short introduction to the spectral curve using the example of harmonic oscillator. Then we introduce the asymptotic Bethe ansatz technique for anisotropic Heisenberg chain in Sec. 3. We solve explicitly for the classical finite-gap solutions and provide 2 concrete examples in Sec. 4 and 5. The semiclassical quantisation is carried out in Sec. 6 and calculations of semi-classical overlaps is given in Sec. 7. We compose the conjecture of the classical-quantum correspondence for the correlation functions in Sec. 8. Finally, we summarise the conclusion and outlook in Sec. 9.

Harmonic oscillator: introducing the spectral curve
Before plunging into the realm of many-body systems, we wish to first familiarise the reader with key concepts of algebro-geometric methods. To this end, we wish carry out the semiclassical analysis on a single quantum-mechanical degree of freedom -the good old quantum harmonic oscillator,Ĥ In the second-quantisation formalism, the quantum harmonic oscillator can be expressed in the diagonal formĤ = ω â †â + 1 2 , (2.3) in terms of the bosonic annihilation operator The ground state (Fock vacuum) |0 satisfiesâ|0 = 0. Excited eigenstates |n , obeyinĝ a †â |n = n |n , are obtained from |0 by successive application of the creation operatorâ † , |n = â † n √ n! |0 , (2.5) such thatĤ |n = E n |n = ω n + 1 2 |n . (2.6) Eigenstates ψ(u) = u|n in dimensionless coordinate u = mω x can be parametrised in term of a normalised logarithmic derivative of the wavefunction called quasi-momentum p(u) (see e.g. [35]), (2.7) Here any subsequently we omitted the subscript n of the wavefunction and energy eigenvalue for clarity of notations. The Schrödinger equation then transforms into a Ricatti equation (2.8) Nodes (i.e. zeros) of excited wavefunctions |n have now turned into simple poles. After a short exercise, one can deduce the following form of the solution [35] with poles u j satisfying a system of equations (2.10) Note that every eigenstate |n gets assigned a unique set of poles u j , with j = 1, . . . , n. From this vantage point, Eqs. (2.10) bear direct analogy to Bethe ansatz equations of integrable interacting quantum systems. Furthermore, we can introduce the Q-polynomial Q n (u) = n j=1 (u − u j ), (2.11) which takes analogous role to the Baxter's Q-function in (integrable) quantum spin chains. By integrating the quasi-momentum we can readily retrieve the wavefunction for the nth excited state with energy E n = n + 1 2 ω, ψ n (u) = 1 √ N e −u 2 /2 Q n (u), (2.12) where N is the normalisation. The solution to the Bethe-like equation (2.10) is given by is the Hermite polynomials H n (u) of order n, Q n (u) = H n (u), (2.13) as expected; poles u j are thus zeros (roots) of Hermite polynomials.

Semi-classical limit
Our next task is to analyse the semi-classical limit of the quantum harmonic oscillator. Before doing so, we start with a brief review of the well-known results of the classical harmonic oscillator (2.14) Here position x and momentum p are phase-space observables obeying the canonical Poisson bracket {p, x} = 1. (2.15) The solution to the classical harmonic oscillator is given in terms of trigonometrical functions p = √ 2mE cos(ωt + φ), x = 1 ω 2E m sin(ωt + φ), (2.16) where phase φ is determined by the initial condition, and E is the value of energy corresponding to the Hamiltonian.
The classical harmonic oscillator can be thought of as the simplest example of a classical integrable system. The associated action-angle pair of variables is given simply by S = E ω , ϕ = ω t, (2.17) satisfying canonical Poisson relation {S, ϕ} = 1.
Lax representation. Canonical notion of integrability rests on the concept of the Lax representation [36]. Here we review the construction for the classical harmonic oscillator. The Lax representation is not unique. One suitable choice for the Lax pair of matrices L and M reads Lax matrix L(v; t) depends on t through x(t) and p(t). As we clarify in turn, additional analytic dependence on the so-called spectral parameter v ∈ C is of pivotal importance for algebraic integrability. One proceeds by recasting the equation of motion of classical harmonic oscillator as the equation of motion of the Lax matrix L(u; t), At first glance it may appear that we have not achieved much. The significance of this reformulation however shows up in the fact the characteristic polynomial is time independent. It defines a spectral curve Σ(w, v) which in our case explicitly reads [36] Σ : Eigenvalues w ± (v) can be conveniently parametrised in terms of a single double-valued complex function called the classical quasi-momentum, satisfying 2 cos p cl (v) = tr (exp L(v)). The quasi-momentum p cl (v) exhibits a pair of square-root branch points at ±1. Function p cl (v) can be interpreted as a two-sheeted Riemann surface with a branch cut along the interval I ≡ [−1, 1] and a pole at infinity, v p = ∞. Due to analyticity it fulfils I dp cl (v) = 0.
Semi-classical limit. We next explain how the classical spectral curve emerges out of a highly-excited quantum eigenstates in which zeros of the Q-function (2.12) densely distribute along the real axis. This corresponds to the semi-classical limit of large excitation numbers n → ∞, with ∼ 1/n. At the classical level, the resulting condensate shows up as a squareroot type branch cut which characterizes the classical spectral curve (2.21). The opposite process -in which the cut disintegrates into a collection of poles -is none other than the familiar Bohr-Sommerfeld quantisation rule Taking the (quantum) quasi-momentum, making the identification u = √ 2nv, and finally taking the limit n → ∞, one finds where we have introduced the distribution of zeros of Hermite polynomial In the large-n limit, one retrieves the classical quasi-momentum Canonical angle. We now describe how to find the angle ϕ associated the spectral curve. This step formally necessitates to find the solution of the auxiliary linear problem in the form of the Lax equations The standard recipe is to identify the canonical angles with dynamical poles of the appropriately normalised eigenvectors ψ ± , called in the literature the Baker-Akhiezer functions [36]. Below we employ a slightly different but equivalent approach using the squared eigenfunction approach. The purpose of this is mostly to manifestly eliminate the normalisation ambiguity of the Baker-Akhiezer vectors. Introducing a 2 × 2 matrix of eigenvectors ψ ± (v, t), we introduce the 'squared eigenfunction' as Presently for classical harmonic oscillator, the solution to the above differential equation admits the following explicit form According to the general rule, the dynamical variable(s) γ j is/are defined as the zero(s) of off-diagonal element of squared eigenfunction Ψ. In the language of algebraic geometry, a set of dynamical variables {γ j } constitutes the so-called dynamical divisor on a Riemann surface. The equations of motion for point of the divisor go under the name of Dubrovin equations. In the case of hyperelliptic algebraic curves of genus g, the number of dynamical variables equals g + 1. In the case of harmonic oscillator, we deal with a surface of genus g = 0 (Riemann sphere), and there is thus a single dynamical variable γ 1 = γ 1 (t). Presently it satisfies a simple evolution law To see this, we perform the follow transformation of the spectral variable which resolves the square-root singularities at v ± and renders λ 2 a rational function of z, (2.35) The square-root branch cut with branch points v ± = ±1, cf. Eq. (2.21), has mapped into 2 regular points on the Riemann surface located at z ∈ {0, ∞}, whereas the original pole at infinity has two pre-images at two punctures z ± p = ±i, one on each Riemann sheet. The dynamical variable z * (t) = exp (ϕ(t)) satisfies has periodic motion with a linearly-evolving angle variable ϕ(t) = ω t, as shown in Fig. 1.

Classical limit of quantum correlation functions
We shall examine how classical quantities arise from the expectation values of quantum observables compute in semi-classical eigenstates, thereby providing an exact classical-quantum correspondence of the correlation functions in the case of harmonic oscillator. The statement is simple: the semi-classical limit of quantum averages corresponds to ergodic averages (Whitham averages [37]) at the classical level where T = 2π/ω is the fundamental period andÔ an arbitrary operator in Heisenberg picture. We shall demonstrate the correspondence with an exact computation on particular examples. Let us start with the operatorÔ =x 2 (0). From the normalized quantum eigenstates, we have which in the semi-classical limit becomes The correspondence continues to hold exactly even thenÔ comprises of non-commuting operators. This illustrate this fact, we consider two observablesÔ 1 =x(0)p(t 0 ) andÔ 2 = p(t 0 )x(0), whose exact quantum averages read Upon taking the semi-classical limit we arrive at The upshot of this elementary calculation is that expectation values of quantum observables in semi-classical limit can be effectively replaced by the corresponding classical observables, with no worry of ordering ambiguities with only enter as 'quantum corrections', i..e subleading contributions (∼ O( )) to the correlation functions.
This section has served mostly for pedagogical purpose, demonstrating the main notations of algebro-geometric methods. The remainder of the article is devoted to a bone-fide interacting many-body paradigm -the quantum Heisenberg spin-1/2 chain with an axially anisotropic interaction (the XXZ model). In distinction to harmonic oscillator, the analytic treatment gets more complicated and in practice the classical limit of correlation functions is rather difficult to access in an exact analytic fashion. We will come back to this in Section 8 after first introducing the formalism and the necessary ingredients.
One distinguishes three phases (regimes): ∆ > 1, ∆ = 1 and ∆ < 1, conventionally called gapped, isotropic and gapless regimes, in respective order, referring to properties of the antiferromagnetic ground state. We shall focus exclusively on the gapped regime and as customary parametrise anisotropy as where q stands for 'deformation parameter' of quantum symmetry algebra U q (su (2)). The model can be diagonalized by Bethe Ansatz techniques. The ground state is two-fold degenerate and we shall consider the ferromagnetically ordered state |↑ ⊗L as the reference pseudovacuum to build the entire spectrum of eigenstates. 1 Each finite-volume eigenstate in the sector with M down-turned spins, denoted by |{ϑ j } M j=1 , is uniquely characterised by a set of (generally complex-valued) rapidities ϑ j called Bethe roots, and they satisfy a coupled system of transcendental equations known as the Bethe equations (BE). They assume the following physical interpretation: the term in the square bracket on the left equals e ip , where p = p(θ) is the magnon momentum , (3.4) whereas the product of trigonometric functions represents a sequence of U (1) scattering phases acquired by a magnon upon a sequence of elastic collision with all the remaining magnons. The total momentum and energy are also given by an additive expression, reading To perform the asymptotic analysis, it is convenient to recast Eqs. (3.3) in terms of the Baxter Q-functions [41] Q(ϑ; which is presently a 'trigonometric polynomial' of degree M whose zeros (in the fundamental domain) corresponds to Bethe roots. Writing shortly Q 0 (ϑ) ≡ sin L (ϑ), and introducing standard compact notations for imaginary shifts, f [±k] (ϑ) ≡ f (ϑ ± iη/2), the Bethe equations can be put in the form or, equivalently, in the logarithmic form 2 Thermodynamic limit. In physical applications one is typically interested in the thermodynamic spectrum of eigenstates. This amounts to study the structure of solutions to the Bethe equations in the L → ∞ limit with M ∼ O(L) → ∞ magnons. These comprise of individual magnons and bound states thereof known as the Bethe strings, referring to compounds of complex Bethe roots whose constituent magnon rapidities are equidistantly separated (modulo finite-sized deviation exponentially small in L) by one imaginary unit, with the centre ϑ (k) ∈ R.
Low-energy scaling limit. There is however another and distinct notion of the thermodynamic scaling limit which allows us to access the low-energy spectrum at long wavelengths. In this limit one takes L → ∞ with M ∼ O(L) → ∞ such that all magnons have low momenta. This way we are left with m ∼ O(1) macroscopic bound states. It is important to stress that in this limit the Bethe's original argument for the string formation is no longer valid. What happens indeed, as we describe in detail below, is that bound state get less dense and generically appreciably deformed. This limit has been first investigated by Sutherland [42]. Further it has been examined in Ref. [43], where the excitations are referred to as 'quantum Bloch walls'. The classical interpretation of such state has been made precise in Ref. [10] for the case of the isotropic Heisenberg chain, establishing explicit connection to classical (nonlinear) spin waves of the Landau-Lifshitz continuous ferromagnet.
In the presence of the interaction anisotropy η, the notion of the "semi-classical limit" amounts in addition to take the weakly-anisotropic limit ∆ 1. More specifically, reinstating lattice spacing a and writing ≡ L a ∈ O(1), anisotropy parameter η ∼ O(1/L) → 0 has to be rescaled as with parameter kept fixed whilst taking the continuum thermodynamic limit, L → ∞ and a → 0. Later on, in section 3.1, we show that length parameter corresponds precisely to circumference of the emergent classical phase space. Expanding the logarithm of Q ± 0 we find where we have assumed that ϑ ∈ O(1). This suggests the following convenient change of variable which yields (3.14) The part involving Q(ϑ) can be treated in a similar fashion. The derivation, which is slightly more subtle, and other details are spelled out in Appendix B. By combining the two contributions, we finally arrive at the following compact representation in the µ-variable In the continuum limit the Bethe roots µ j condense along certain (in general disjoint) onedimensional segments (contours), that is C ≡ ∪ j C j , and one accordingly obtains a singular integral equation of the form with integral kernel of the form As it turns out, this correction to asymptotic Bethe equation (3.16) can only be safely omitted provided π (µ 2 + δ)ρ(µ) = iπn, n ∈ Z, (3.19) which holds when the density of roots near the real axis is relatively low. In contrast, for high enough densities the above perturbative analysis is no longer justified and one has to take into account the quantum Bethe equations in a non-perturbative fashion, at least in the region near the real axis. As discussed later on in Section 6.2, when this scenario takes place one find certain special types of solution called condensates [11].
Riemann-Hilbert problem. The asymptotic Bethe equations (3.16) can be universally formulated as a Riemann-Hilbert problem. To this end, we define the spectral resolvent (3.20) and write On each density contour C j function p(µ) experiences a jump discontinuity proportional to the density (of Bethe roots) ρ(µ) at that point, namely with ±i0 denoting infinitesimal shifts on either side of the contour. 3 In this view, it is natural to adopt the convention that C j is the jth square-root type branch cut of a two-sheeted Riemann surface where end points of C j are identified with its branch points. Function p(µ) is thus a double-valued complex function which is analytic in the finite part of µ-plane except on C j . Crossing a square-root type branch amounts to flip a sign of p(µ), i.e. jumping to another Riemann sheet. In the process p(µ) in general picks up an integral multiple of 2π, namely In the next section we shall explicitly demonstrate that p(µ) is indeed precisely classical quasimomentum which encodes the spectrum of monodromy matrices associated to a completely integrable classical anisotropic ferromagnet. We notice that occasionally it can be more convenient to write the Riemann-Hilbert problem in terms of spectral parameter ζ = 1/µ. This proves especially useful in the Sections 6 and 6.2, in order to compare the branch cuts to the distribution of Bethe roots with finite but large system sizes. A dictionary between the two sets of equation with different spectral parameters is given in Appendix A.

Anisotropic Landau-Lifshitz field theory
The task at hand is to derive the classical equation of motion which governs the low-energy spectrum of the XXZ Heisenberg chain. This can be systematically achieved by taking the semi-classical expansion of the quantum monodromy matrix and inferring the spatial component of the classical Lax connection. This will allow us to explicitly establish that p(µ) introduced earlier is indeed classical quasi-momentum that encodes the spectrum of solutions to the classical axially anisotropic Landau-Lifshitz model for an S 2 -valued classical spin vector field S(x, t).
Here and subsequently we have employed compact notation S t ≡ ∂ t S(x, t), S x ≡ ∂ x S(x, t), and likewise for higher partial derivatives.

Equation of motion (3.24) is generated by the Hamiltonian
(3.25) 3 We adopt the orientation of integration contours that go from the branch point with positive imaginary part to the branch point with negative imaginary part.
Zero-curvature representation. The integrability manifests itself in representing the equations of motion as the zero-curvature condition (3.26) which plays the role of the compatibility condition of the auxiliary linear problem, with [14] U(µ; x, t) = 1 2i representing the spatial and temporal components of the classical connection, and denotes the spin current density at δ = 0.
Semi-classical limit. The starting point of the derivation is the quantum Lax operator which serves as the building block for commuting transfer matrices and algebraic diagonalisation of the XXZ Hamiltonian (3.1). The fundamental Lax operator L(ϑ) acts on the physical Hilbert space V p ∼ = C 2 of a spin-1/2 degree of freedom and an auxiliary space V a associated to the defining representation of the quantum group U q (su(2)) (with deformation parameter q = e η ), reading depending on the complex (spectral) parameter ϑ. The auxiliary spin generators enclose the following q-deformed commutation relations The fundamental row transfer matrices of the quantum XXZ spin chain is defined as an ordered product of Lax matrices (called the quantum monodromy matrix) partially traced over the auxiliary space V a ∼ = C 2 , where we have adopted the right-to-left ordering convention and used L (k) (ϑ) to denote the embedding into the kth factor in an L-fold product Hilbert space H ∼ = V ⊗L p . Quantum Yang-Baxter relation ensures that matrices T (ϑ) are in mutual involution, [T (ϑ), T (ϑ )] = 0 for all ϑ, ϑ ∈ C; they thus serve as the generating operator for all local and nonlocal conserved charges [10,39]. One can moreover construct analogous transfer matrices T j (ϑ) with (j + 1)dimensional auxiliary unitary representations of U q (su(2)), commuting amongst each other, and providing the quasilocal conservation laws of the XXZ model. While these are rather important for the conventional thermodynamic considerations (see e.g. Refs. [44][45][46][47]), they will nonetheless not play a role in taking the semi-classical limit.
We now perform the asymptotic scaling limit L → ∞ as previously at the level of the Bethe equations, which amounts to take η = /L → 0 and q → 1. We are thus allowed to substitute the q-deformed spin generators with fundamental su(2) spins, S α →Ŝ α = 1 2σ α for α ∈ {x, y, z}. In this limit, the diagonal elements of the quantum Lax operator become Rewriting it in terms spectral parameter µ = tan ϑ , (3.36) and reinstating the lattice spacing a = /L, we have The asymptotic scaling limit of the quantum monodromy matrix M(µ) is thus given by the following path-ordered product Finally, by regarding the quantum spinsŜ α with classical spin field variables S α j ≡ S α (x = j a), and subsequently taking the continuum limit, the semi-classical approximation of the quantum monodromy matrix turns into a path-ordered integral (a Wilson loop which winds once around the system) where U(µ; x, t) is none other than the spatial connection associated to Eq. (3.28).

Finite-gap integration method
After having established that classical quasi-momentum p(µ), formally representing a singlevalued complex function on a two-sheeted Riemann surface In the previous section we have explained how it emerges out of the semi-classical spectrum of the anisotropic Heisenberg chain, describing the low-energy solutions of the asymptotic Bethe equations. In this section we shall proceed to explain how its analytic structure encodes physical properties of interacting nonlinear phases of classical spin-wave solution. We confine ourselves to solutions that involve a finite number of excited nonlinear modes (in correspondence with Riemann surfaces of finite genus), commonly known in the literature as the finite-gap solutions [29,48,49]. The complete programme for the algebro-geometric integration of nonlinear partial differential equations permits to solve the initial value problem consists of the following main steps: 1. prescribing an appropriately normalized meromorphic differential dp on a finite-genus (hyperelliptic) Riemann surface, 2. constructing the fundamental solution to the associated auxiliary linear problem, 3. identifying the dynamically separated variables and deriving their equations of motion, 4. constructing canonical action-angle variables that undergo linear evolution law on a Jacobian with aid of the Abel-Jacobi transformation, 5. expressing physical fields through the solution of the inverse problem.
The finite-gap integration scheme applies universality, and has already been previously carried out for the particular case of the Landau-Lifshitz field theory in Refs. [50,51]. We shall nonetheless describe a slightly different formulation which we find conceptually simpler, avoiding the conventional use of the Baker-Akhiezer functions. We discuss some general aspects and provide two explicit realisations for g = 0 and g = 1. We do not focus on the action-angle variables but rather give an explicit time dependence of the physical spin fields in terms of the Riemann Θ-functions.

Auxiliary linear problem
The task at hand is to solve the auxiliary linear problem (3.27). Formal integration along the spatial direction at a fixed time-slice yields By imposing periodic boundary conditions S(x + ) = S(x), we define the monodromy ma- According to the Bloch theorem, there are two linearly-independent solutions that linearise this translation, where Bloch multiplier Λ(µ) = diag(Λ + (µ), Λ − (µ)) is given by eigenvalues of M cl (µ), which (by virtue of the zero-curvature condition (3.26)) are conserved in time, ∂Λ ± (µ)/∂t = 0. Since Tr M cl (µ) = 0, monodromy matrix M cl (µ) is unimodular and hence its eigenvalues can be parametrised in terms of a single quasi-momentum p(µ) in the form Λ ± (µ) = exp (±ip(µ)). Subsequently we make use of the following general parametrisation M cl (µ) = cos (p(µ)) + i sin (p(µ))Ψ(µ; x 0 ), (4.4) where analogously to Section 2 we introduce the squared eigenfunction The latter is periodic, Ψ(µ; x + ) = Ψ(µ; x) and satisfies the adjoint linear system Uniformisation. Introducing an uniformised spectral parameter z, the solution to the above linear problem can be sought as a formal Laurent series the first few terms can be written compactly in the form and so forth.

Local conserved charges.
Quasi-momentum can be express through the squared eigenfunctions via [14] p where coefficients Q n provide (extensive) local conserved charges of a finite-gap solution. The first two take the form (modulo total derivatives) which yields

Finite-gap solutions
In order to fully reconstruct the spin field S(x, t) from the spectral data we need to know more than just the quasi-momentum p(µ) (storing information about conserved charges of the solution). We also need dynamical degrees of freedom associated with the algebraic curve.
We outline the algebro-geometric integration procedure which provides time evolution of the spin field S(x, t). The corresponding field configurations are called finite-gap solutions.
To achieve this, we parametrise the squared eigenfunctions as specifies a hyperelliptic algebraic curve Σ : The curve is fully characterized by 2(g + 1) branch points µ j (or, equivalently, symmetric polynomials r k thereof, with r 2g+2 = 1). Furthermore, by the unimodularity constraint det Ψ g = −1, functions a(µ) and b(µ) are not independent but satisfy an algebraic relation From the trace identity (4.11) one can readily infer the following compact expression for the quasi-momentum . (4.18) This form guarantees the correct asymptotic expansion about µ → ∞. It will produce only g + 1 functionally independent integrals of motion Q n . They can be expressed as certain functions of coefficients r j of polynomial R 2g+2 (µ). Moreover, the information about total filling fraction ν with respect to ferromagnetic vacuum S z vac = 1, namely can be obtained as (4.20)

Dynamical divisor
Off-diagonal elements of Ψ g (µ), i.e. dynamical zeros of the b-function, provide dynamical degrees of freedom of finite-gap solutions 4 , and allow to reconstruct the physical field S(x, t).
To satisfy quadratic constraint (4.17), we parametrise Here the leading coefficient has been fixed by inspecting the asymptotics at µ → ∞.
The set D = {γ j } g j=1 is known as the dynamical divisor of the Riemann surface Σ. We note that S − (x, t) is also an independent dynamical variable, so we have in total g + 1 dynamical degrees of freedom, matching the number of action variables (also equivalent to the number of forbidden zones in the finite-gap spectrum).

Dubrovin equations.
Space-time evolution of dynamical zeros γ j = γ j (x, t) takes place on a Riemann surface Σ. It proves useful in practice to extend the dynamical divisor D ext by adjoining to D two extra non-dynamical variables γ ± ≡ ±i , denoted as γ g+1 and γ g+2 respectively. Then, using Lagrange interpolation formula, we can restore functions a g+1 (µ) from Eq. (4.17), yielding .
The equations of motion (4.6) in terms of γ-variables from D ext take the form This system of differential equations, governing the motion of the dynamical divisor on a Riemann surface, goes commonly under the name of Dubrovin equations [29,48,53]. These equations are universal, i.e. do not depend on the model under consideration. What is model-specific is how γ-variables relate to physical fields, called the reconstruction formulae. Dynamics of the S 2 -values spin field, say S z and S − components, can be restored from Eqs. (4.6), yielding , (4.25) and , .

(4.27)
Abel-Jacobi transformation. In order to solve Dubrovin equations (4.23), we define the standard basis of 2g closed cycles, comprising of A-cycles and the associated canonically conjugate B-cycles, using the following prescription: A j -cycle encircles the the jth branch cut C j on the upper Riemann sheet of Σ, whereas B jk denotes a cycle that passes through cut C j on the upper sheet and connects back it itself through C k , as shown in Fig. 2.
Dubrovin equations for the dynamical divisor on Σ can be explicitly integrated with aid of the Abel-Jacobi transformation, where ω j form the basis of holomorphic differentials of the Riemann surface. The above mapping provides a variable transformation from γ-variables to the angle-type ϕ-variables, The basic holomorphic differentials ω j are formally the form for an appropriate choice of coefficients C jk , normalised with respect to A-cycles Taking into account Eq. (4.28), the equations of motion of the divisor, Eqs. (4.23) linearises. Indde, making use of the Lagrange interpolation formulae, one can find with wavenumbers k j are frequencies w j reading Dynamics of linear phases ϕ j (x, t) exercises a quasiperiodic motion on a Liouville torus T 2g of real dimension 2g. For any physically admissible initial condition, γ-variables evolve along closed trajectories which are homotopically equivalent to A-cycles A j , such that there is precisely one variable per cycle C j .
Periodicity constraints. Periodic (i.e. closed) solutions S(x) = S(x + ) are further distinguished by the periodicity of angles, ϕ j (x + , t) = ϕ j (x, t) + 2π n j where integers n j ∈ Z specify the winding numbers of assigned to each branch cut. These are mode numbers associated to finite-gap solutions. Similarly, invariance under translation for a temporal period T implies quantisation of frequencies w j . Under these extra conditions, coefficients C j1 and (r 1 /2)C j1 + C j2 become integer-valued, thus restricting the admissible algebraic curves. Figure 2: Assignment of cycles on the two-sheeted Riemann surface, illustrated on two cuts C j and C k . After tunnelling to the other Riemann sheet, the integration orientation is reversed [10].
Extra degree of freedom. There is an important subtlety in the construction above: the separated γ-variables which form the divisor D ≡ {γ j } g j=1 indeed do not comprise the complete set of dynamical variables of the finite-gap solutions. There is an additional degrees of freedom which is not included in the canonical set of angle variables (4.28) obtained from dynamical zeros γ j of the b-functions. This fact becomes clearly apparent already in the simplest class of g = 0 solution which, as we shall demonstrate in a moment, in fact corresponds to a nontrivially evolving spin field with a single linearly-evolving angle variable. Generally, that is in the case of algebraic curves of genus g, there are thus g+1 phases (that is precisely the number of cuts, i.e. half the number of branch points) and accordingly the linearised dynamics takes place in the Liouville torus of real dimension 2(g + 1).
By substituting the form of the b-functions, Eq. (4.21), into Eq. (4.18), and using the restoration formulae (4.25), (4.26), (4.27), one can establish the identify By introducing an auxiliary function we can therefore present the quasi-momentum as a period integral of the form . (4.36) For g = 0 expression, we put Eqs. (4.21) and (4.22) into Eq. (4.18), which yields Comparing Eq. (4.36) with Eq. (4.26) revealing that the following integral provides an extra winding number n g+1 . Taking into account that γ-cycles are topologically equivalent to A-cycles that loop around the cuts of Σ, we can write where n k is the winding number of the kth cut. While on a surface Σ quasi-momentum p(µ) is a single-valued function, when viewed as a function on the complex spectral plane, namely for µ ∈ C, it is a double-valued function which experiences jumps upon traversing one of the cuts; one ends up on another Riemann sheet, and thus R(µ) flip sign. Denoting points on the either side of the cut C k by µ ± i0, the crossing condition for each cut reads which are known as the Riemann-Hilbert problem. The additional winding number n g+1 is encoded in the asymptotic condition for p(µ) at µ = ∞ ± , Choosing a new basis of open B-cycles B j (for j = 1, 2, . . . , g + 1), starting on the upper sheet at ∞ + and connecting to ∞ − on the lower sheet by passing through the cut C k , as demonstrated in Fig. 2, the equations (4.40) alongside Eq. (4.41) can be compactly stated as An alternative approach. Eq. (4.39) is by construction a solution to the Riemann-Hilbert problem (4.40). However, the behaviour at µ → ∞ does not automatically have the form of Eq. (4.13). Instead, one has to demand certain conditions on the spectral curve, which are identical to the integrality condition of wavenumbers (4.33). Alternatively, we can construct quasi-momentum p(µ) that has the correct asymptotics at µ → ∞ by considering the meromorphic (second kind) differential with the prescribed asymptotics at µ → ∞.
As we have presented above, quasi-momentum p(µ) as a solution to the Riemann-Hilbert problem (4.40) can be obtained as integrals in the real space related to the dynamical divisor (4.36). In fact, we could construct p(µ) purely from the analytic properties and asymptotics. And this fixes two coefficients of P 2g+2 (µ), namely The remaining m = g − 1 coefficients are fixed by And the spectral curve is fixed by the condition (4.42). One can show that there is a unique differential with the properties above.

One-cut rational solution
The simplest class of solutions of the Riemann-Hilbert equations (3.23) is associated with algebraic curves of genus zero, corresponding to Riemann surfaces with a single branch cut. We fix the corresponding polynomial to be The branch points µ 1 ,μ 1 ∈ C are chosen based on reality condition. The solution is therefore involves two real degrees of freedom.
In the following derivation we set the classical period to = 1. With this choice, the admissible value of the wavenumbers is k = 2π n for integral n, and consequently the branch points are get "quantised" too.
Quasi-momentum p(µ) is given by Eq. (4.37). To satisfy "quantisation condition" and to obtain the prescribed filling fraction (4.20), we have allowing us to express the branch points as and the algebraic curve in terms of mode number n and filling fraction ν thus reads The associated quasi-momentum, cf. Eq. (4.37), can be written in the form Alternatively, the quasi-momentum can be found as the integral of the meromorphic differential on the rational curve (cf. Eq. (4.43)). In g = 0 case, such differential is described with coefficients c 0 and c 1 , dp 1 (µ) dµ = − 1 2 Both coefficients are fixed by the asymptotics, see Eq. (4.44), providing c 1 = 1 and c 0 = −(µ 1 +μ 1 )/2. By taking an integral of dp 1 (µ), we correctly recover the quasi-momentum (5.5).
Notice that presently there is no canonical A-cycles. The wave number can be retrieved by evaluating the 'open' B-cycle Here n ∈ Z is the mode number of the solution, physically representing the winding number. Phase-space averaged of local charges evaluated on a one-cut solution can be read off from asymptotic expansion of quasi-momentum p 1 (cf. Eq. (4.13)) µ → ∞ : The initial two non-trivial charges are momentum and energy, reading explicitly Having parametrised the one-cut quasi-momentum, we can now reconstruct the spin field S(x, t). We remind that, since the genus of the algebraic curve is zero, there is no dynamical γ-variable involved. The only dynamical degree of freedom is thus the transversal component S − (x, t). We first use Eq. (4.25) to deduce the longitudinal component where we have used the frozen (i.e. non-dynamical) γ-variables γ 1 = i and γ 2 = −i . Next, we solve Eqs. (4.26) and (4.27) to find the transversal component of the spin field, By imposing normalisation constraint | S(x, t)| = 1, we finally arrive at the following general form of the one-cut solution where the wavenumber and frequency are k = 2πn, w = (4π 2 n 2 + δ) cos θ 0 . (5.13) The momentum and energy can be alternatively computed by direct integration using Eq. (4.13)

One-cut solution from asymptotic Bethe ansatz
In the next section 6 we shall be concerned with semi-classical quantisation of finite-gap solutions. At this stage, we wish to discuss the the one-cut solutions from the vantage point of low-energy modes in the Heisenberg anisotropic spin chain. Our aim is to identify solutions to the Bethe equations (3.3) which show up classically as one-cut solution. We achieve this with the method outlined in Appendix E.1. Specifically, we are after low-energy solutions which comprise of macroscopically O(L) excited magnons which condense into a single mode. Fixing the mode number n = 1, we compute the total momentum P and total energy E of such state for different system sizes L but with a filling fraction ν = M/L fixed, see Eqs. (3.5).
The results are shown in Fig. 3.
For the classical 1-cut solution we compute the momentum and energy as in Eq. (5.9) with classical period = 1.
While quantum momentum in the classical limit yields directly its classical counterpart, to compare energies there is an extra rescaling factor of L to take care of, namely H · L ∼ H. At the classical level, one-cut solutions describe simple precessional dynamics which may be view as a finite-density analogue of Goldstone modes (that is elementary ferromagnetic excitations known as magnons). Remarkably, later on in Section 6, we demonstrate that quantisation of such solutions can gives rise to certain non-perturbative effects (first discussed in Refs. [11]) which can only be explained by property taking quantum corrections into account.

Bions from degenerate two-cut solutions
The next class of solutions in terms of complexity is periodic elliptic magnetisation waves which go under the name of 'cnoidal waves', encoded in elliptic algebraic curves (described by Riemann surfaces with two branch cuts, i.e. of genus g = 1). In the limit of unit elliptic modulus, the elliptic profiles reduce to trigonometric ones and one retrieves the well-known soliton solutions. In the easy-axis anistropic Landau-Lifshitz model, special types of two-cut solutions are known as bions, representing two-mode bound states of a kink and an anti-kink. A special degeneration of such a bion in the limit of large circumference produces a stable kink mode. The latter is for instance responsible for the freezing of a magnetic domain wall, as demonstrated in [16]. Subsequently in Section 6.3 we will study semi-classical quantisation of the classical bion solution. Before that, we derive it here using the outlined integration procedure.
The elliptic curve encoding the bion spectrum has the form parametrised by two pairs of complex-conjugate branch points located on the imaginary axis Re(µ) = 0 at µ j ∈ {±iξ 1 , ±iξ 2 }, satisfying conditions This means that the bion solutions can only be found in the regime δ > 0, i.e. easy-axis regime. In what follows we put, mostly for simplicity, δ = = 1. In fact, from the classical equation of motion for the spin field S(x, t) given by Eq. (3.24), the solution at δ = 1 can be mapped to another solution with δ > 0 by a simple rescaling We next outline how to reconstruct of the classical spin field from the algebraic curve. Before that, we introduce the standard full elliptic integrals, cf. Appendix C, In addition, when the argument of the elliptic function is omitted we adopt that k = ξ 2 /ξ 1 . Since the surface involves 2 branch cuts, there is a single dynamical γ-variable γ 1 (x, t). As described earlier, there are also 2 extra non-dynamical ones pinned at locations γ 2 = i and γ 3 = −i (recall that = 1). From the Dubrovin equations (4.23) we find where we have used sn(x + iK 2 )sn(x)k = 1, and use to denote the integration constant.
Remarkably, it turns out that γ 1 is static.
Using the reconstruction formulae (4.25) for S z component, we thus have The integration constant can be fixed by imposing the reality condition S z ∈ R, which yields a time-independent profile The solution is periodic with period where ±n are the mode numbers associated with the two cuts. Variable γ 1 can be therefore expressed as The other component of the spin field, say S − (x, t), can be reconstructed with aid of formulae (4.26) and (4.27), i.e.
. (5.26) Using the properties of elliptic functions, we have where an auxiliary function Θ(x) is defined as (cf. Appendix C) Finally, constant C is can be fixed by requiring S 2 = 1, yielding where φ 0 ∈ R is a phase that is determined by the initial condition. A representative instance of a bion solution is shown in Fig. 4. Static kink. We now consider a special degeneration of a bion solution by 'blowing up' the period . Such a limit is typically called the 'soliton limit'. This is achieved by sending both branch points of R 4 (µ) to i , that is ξ 1,2 → (presently = 1), which corresponds to collapsing the two branch cuts. In order to perform the soliton limit, we first shift the argument of the elliptic function by quarter period K 1 /ξ 1 ; for instance the S z field takes the form (cf. Eqs. (C.6)) S z (x) = −ξ 1 cn(ξ 1 (x + K 1 /ξ 1 )) dn(ξ 1 (x + K 1 /ξ 1 )) = ξ 1 sn(ξ 1 x). (5.30) Now taking the limits ξ 1,2 → 1 and accordingly k = ξ 2 /ξ 1 → 1, we find which is a static kink! The transversal components can be easily deduced from the equation of motion (3.24), yielding where φ 0 ∈ R is a phase determined by the initial condition.

Semi-classical quantisation
We have outlined how classical spectral curves encode complete information about the nondynamical part of the finite-gap solutions. Moduli of algebraic curves are completely fixed by locations of the square-root branch points, i.e. roots of polynomial R 2g+2 (µ), which (by the reality condition) must always appear in complex-conjugate pairs. Local conserved charges are expressible as functions of symmetric polynomials r j of branch points that are coefficients of R 2g+2 (µ). Before we continue, an important remark is in order. While branch points {µ j } are directly linked physical quantities, the square-root branch cuts pairwise connecting them are merely a matter of convention, as there is plenty of freedom in assigning branch cuts to given branch points. The prevalent choice in the finite-gap literature [36] is to place the cuts along straight (vertically aligned) lines connecting every complex conjugate pair of branch points, which are in turn encircled by the canonical basis A-cycles. It suffices to say that purely from the standpoint of classical finite-gap solution this convention is perfectly adequate. Semi-classical quantisation of finite-gap solutions nevertheless reveals a physically distinguished choice for branch cuts, where the surface is cut precisely along disjoint segments which represent the forbidden zones of the classical transfer function [10]. In the view point of the quantum spin chain, these are contours on this magnons (Bethe roots) condense. As demonstrated in turn, the physical cuts can significantly differ from the conventional straight cuts and even undergo certain condensate phenomenon. Computing the Bethe root densities along the physical contours is thus the essential step to perform a semi-classical quantisation.

Physical contours
We shall now describe a general procedure to determine the physical contours. We employ the algorithm proposed and implemented in Ref. [11]. Our aim is to to facilitate a direct comparison with exact low-momentum quantum eigenstates at finite L in the weakly-anisotropic regime. We find it most convenient to carry out this computation in ζ-plane 5 by applying the following anti-holomorphic transformation 6 µ → ζ(µ) : Firstly, we introduce the complex density function where n j ∈ Z is the mode number associated with cut C j . The tasks is to determine C j . In fact, by virtue of the second equality in Eq. (6.2), the density function ρ(ζ) can be defined on the entire Riemann surface. As a direct consequence of p(ζ = ζ ) = πn j , at the square root branch points ζ ∈ {ζ j ,ζ j } the densities vanishes, ρ(ζ ) = 0. In a small neighbourhood around it, one finds ρ(ζ = ζ + ε) = O( √ ε). Away from the branch points ρ(ζ) in general takes complex values.
The defining condition for the physical contour C j is prescribed as follows: starting from a given branch point, the integrated density differential ρ(ζ)dζ always remains real, Physical interpretations of this prescription is rather evident; physically ρ(ζ)dζ corresponds to the number of excitations (Bethe roots) within an infinitesimal interval in ζ-plane, which is a positive definite quantity by definition. Notice that this condition alone is not yet sufficient.
In particular, one finds that there are three distinct contours that emanate out of each branch point which comply with such a positivity requirement [11]. It turns out that one of the solutions has infinite filling fraction, and one can immediately discard it on this basis. From the remaining two contours, only one is physical. The defining condition is that the total filling fraction does not exceed the threshold value of maximal total filling ν max = 1/2, that is Semi-classically quantising a reference finite-gap solution amounts to dissolve each cut into a large (but finite) number of individual magnons. This is achieved by requesting that for a spin chain of large system length L, individual excited low-energy (nonlinear) spin-waves carry integer quanta of magnetisation M j ∼ O(L), such that (i) M j /L ≈ ν j and (ii) the Bethe roots are distributed approximately with density ρ j (ζ) along C j . Here we wish to contrast this with exact quantisation which instead takes quantum fluctuations fully into account to all orders. In other words, the obtained semi-classical solutions only serve as an approximation of finite-volume quantum-mechanical eigenstates. For a full non-perturbative (i.e. exact) quantisation, solving the Bethe ansatz equations (3.3) cannot be avoided.

Single contour at low density.
In order to give the simplest benchmark of the described procedure, we proceed by illustrate how a classical one-cut solution emerges as a semi-classical eigenstate in the anisotropic gapped Heisenberg spin-1 2 chain. For simplicity, we first consider the case when the filling fraction of a physical cut is low, mostly to ensure that the finite-size effects (cf. Eq. (3.18) and (B.17)) become negligible at large systems size. We then find that in the continuous limit (L → ∞) the Bethe roots distribute along the computed contour as shown in Fig. 5. More specifically, we set anisotropy to δ = 1 and choose filling fraction ν = 0.1 with mode number n = 1. Using the above prescription we next compute the density contour satisfying Eqs. (6.3) and (6.4), as shown in Fig. 5 (the numerical solutions to the Bethe equations (3.3) corresponding to the quantised 1cut solution are explained in Appendix E.1). Upon taking the L → ∞ limit, the semi-classical eigenstate will eventually be given by a dense arrangement of roots which distribute densely along the contour specified by the conditions (6.3) and (6.4).
Increasing the filling fraction ν, we observe that 'quantum fluctuations' (contained in higher order terms in the ABE) progressive amplify until eventually begin to produce condensates.

Formation of condensates
As announced previously, a certain type of critical phenomenon sets it when the maximal density along the contour exceeds the critical value, responsible for formation of condensates of uniform density. The critical condition is signalled by a divergence of the finite-size correction given by Eq. (3.18) (see also Eq. (B.17)).
We shall first examine the phenomenon of the one-cut solution working in ζ-plane. By starting from a relatively low filling fraction ν and gradually increasing it, the value of function ρ(ζ)(1 + δζ 2 ) on real axis approaches the value of i. Once this value is reached at critical value For fillings ν > ν crit , the density contour acquires a straight vertically segment of unit uniform density which is first produced on the real axis and then begins to expand outwards. To the best of our knowledge, such on object has been first identified in Ref. [10], where it is referred to as a condensate. From the viewpoint of the spin-chain spectrum, the spacing between constituent magnons is always equal to i; hence one can think of a condensate as a giant regular Bethe strings. In the complex spectral space of finite-gap solutions, condensates are associated to branch cuts of logarithmic type. 7 Despite of the appearance of additional condensate for ν > ν crite , the location of the physical contour can still be determined solely by knowledge of a classical finite-gap solution by taking into account conditions (6.3) and (6.4).
Emergence of uniform condensate is intimately related to the notion of so-called 'fluctuation points', which play a key role in classical modulation stability theory [55]. To be concrete, one can imagine a reference finite-gap solution with m cuts, and denote the occupied mode numbers {n 1 , n 2 · · · n m }. To excite a mode with an unoccupied n, say n * , the following condition for the quasi-momentum p(ζ) has to be satisfied by virtue of periodic boundary conditions, p(ζ m,n * ) = n * π, (6.5) 7 Logarithmic branch cuts of the same type likewise get produced in a degeneration process of coalescing two square-root type branch points (which conventionally corresponds to merging to nearby standard cutsa process which yields soliton modes. In this situation, finite-gap quasi-momentum is no longer meromorphic. The condensates we describe here are different as the quasi-momentum differential remains of second-type all the way through). where ζ m,n * label distinct fluctuation points. Note that these can either be real, or may occur in complex-conjugate pairs (owing to the square-root branch cut nature of the quasimomentum). Fluctuation points have a clear physical meaning: they are small fluctuations of a reference finite-gap solution which possess infinitesimal filling fractions. Upon increasing their filling fraction, they continuously grow up into nonlinear finite-gap mode.
Alternatively, one can say that fluctuation points are merely "almost degenerate" branch cuts, so small that they do not affect the form of the quasi-momentum p(ζ). An interesting question emerges whether classically a given finite-gap solution is stable due to the existence of these fluctuations, see for instance discussion in Ref. [55]. The upshot here is that the stability condition coincides with the condition for the formation of condensate in the semi-classical limit.
In the following we shall first take a look at the basic case with a single cut. There we find the physical contours of which the Bethe roots distribute consist of three pieces: two parts of to the contours which connect to the square-root branch points are joined by a uniform condensate attached to the 'the middle', with two additional curved contours emanating from the intersection points connecting to the nearby fluctuation point. We give an explicit demonstration of this scenario in Section 6.2.1. Presence of multiple excited cuts makes the situation even more involved as cuts exert among themselves an effectively attractive interaction. This situation is described in Section 6.2.2.
We mention in passing that a somewhat reminiscent phenomenon is known to appear in the context of matrix models [56] which are described by similar Riemann-Hilbert problems, or in other context such as in large-N Yang-Mills theory [57][58][59], and random tiling models [60,61]. They commonly go under the name of the Douglas-Kazakov phase transition [59], a variant of a third-order phase transition. More interestingly, the condensate appears when taking the semi-classical limit of the attractive Lieb-Liniger model [62,63]. In that case, the quantum phase transition can be detected via the calculation of ground state correlation functions [63]. In our case however, we do not deal with any real phase transition: the quasimomentum p(ζ) and branch points do not undergo any discontinuous change, unlike in the case of Douglas-Kazakov transition where the free energy experiences a change after forming a "condensate" [61].

One-cut case with condensate
In the case of a single cut, we have with the associated filling fraction ν 1 < 1/2. Locations of fluctuation points, which we denote below by ζ 1,k , can be inferred from the density , p(ζ 1,k ) = kπ, (6.7) Mode number n = 1. For the basic illustration, we consider the one-cut solution with mode number n = 1. So far the filling fraction does not exceed the critical value, we end up finding a single smooth umbrella-shaped contour, as exemplified in panel (a) in Fig. 6. The nearby fluctuation point sits a finite distance away, i.e. somewhere left from the cut, and the density at the real axis satisfies condition Increasing the filling fraction will also increase the density at the intersection with the real axis. In the process, the nearest fluctuation point approaches the physical contour and at the critical filling collides with it at ζ * , given by ρ(ζ * )(1 + δζ 2 * ) = i. (6.9) This event is shown in panel (b) in Fig. 6. 8 If the filling fraction keeps increasing further, the fluctuation points after collision 'tunnel through' the cut and leaves behind a vertical condensate. The nearest fluctuation point on the real axis will thus be located to the right of the physical cut, as pictured in panel (c) in Fig. 6. One can nonetheless recover the same quasi-momentum p(ζ) by considering an additional pair of contours which emanate out of the fluctuation point(s), satisfying ρ 2 (ζ)dζ ∈ R with density defined through Eq. (6.7) (depicted by brown dashed lines in panel (c) in Fig. 6). Interestingly, upon emergence of the condensate, the original contour cannot accommodate for all the Bethe roots, and some "excess" Bethe roots will lie along those additional contours. We emphasise again that the quasi-momentum p(ζ) remains intact.
The above picture offers the following suggestive explanation. One can see the present one-cut solution as a limiting (degenerate) case of a more general two-cut solution when one of the cuts becomes infinitesimally small, namely get 'switched off' to a fluctuation point. This is neatly captured by panel (a) in Fig. 7, where the blue solid lines represent parts of the original physical cut attached to branch points (ζ 1 ,ζ 1 ), while the red dashed line shows the condensate of Bethe roots with uniform density. The extra green solid line can be understood as part of one of the "unphysical contours" 9 associated with the infinitesimal branch cut (ζ F ,ζ F ). Combining the three ingredients, we are therefore able to determine the densities of Bethe roots along these contours. This amount to account for the leading-order quantum corrections to ABE (3.16) in non-perturbatively fashion.
To explicitly corroborate the above picture, we compared the computed contours with a numericals solution to the Bethe equations with finite but large system sizes, shown in Fig. 7. Moreover, we have also checked the proposed contours quantitatively through the calculation of the leading order of the Gaudin norm of Bethe states for the condensate formation, depicted in Fig. 13. We emphasise that physical contours are of vital importance for the functional integral approach to computing overlaps and the norms between semi-classical eigenstates, which yields an opportunity to verify whether the described contours are indeed correct. The numerical evidence is collected in Appendix D.1.
Mode number n ≥ 2. One can encounter an even more exotic situation for the mode number n ≥ 2 such as two fluctuation points forming a complex conjugate pair [11]. To understand why this happens, notice that p(ζ * ) = (n + 1)π only has one real solution for mode number n = 1. In contrast, for mode numbers n ≥ 2 it is possible to have a complex conjugate pair of solutions. In the latter scenario, the additional contour can be found via the same condition, i.e. ρ n+1 (ζ)dζ ∈ R, with uniform condensate located between the intersection points, similarly as in the case discussed above. This time however, the contours can be classically seen as arising from a three-cut solution with one large physical cut and two nearly degenerate infinitesimal cuts located at the complex-conjugated fluctuation points ζ F and ζ F . For instance, in Fig. 8 we demonstrate the situation for the isotropic Heisenberg spin case with mode number n = 2. It is unfortunate that the employed numerical method for generating analogous solutions in the anisotropic ferromagent does not work in the case n ≥ 2, cf. Appendix E.1. Despite that, since the Bethe roots do not change much when increasing η from 0 to L , we expect the phenomenon to survive in the presence of weak interaction anisotropy with η ∼ O 1 L .

Multiple cuts
To make the matter more complicated, in the case of multiple cuts the condensates can appear not only out of fluctuation points mechanism but also due to the interaction between two physical cuts. Here we focus for simplicity on the two-cut case, since scenarios with more cuts can be described in parallel with the phenomenology of the two-cut case. An exhaustive survey on the two-cut case at isotropic point (δ = 0) has been conducted previously in Ref. [11]. The anisotropic case analysed below with anisotropy η = L > 0 is pretty much similar. Nonetheless, we still would like to offer a few intuitive physical explanations to some scenarios. Firstly, when two physical cuts are well separated, each branch cut can produce a condensate when their filling fractions attain certain critical value due collisions with the fluctuation points, in analogy with the 1-cut case [11]. However, when the two physical cuts get even closer to each other, they can become glued to each other via a condensate.
A basic example of the above phenomenon occurs when there are two cuts have consecutive mode numbers, i.e. n 2 = n 1 +1. When they get sufficiently close to one another, their physical contours intersect at points ζ int andζ int , at which their combined density satisfies the following identity [ρ n 1 (ζ int ) + ρ n 1 +1 (ζ int )] 1 + δζ 2 int = ρ n 1 (ζ int ) + ρ n 1 +1 (ζ int ) 1 + δζ 2 int = i. (6.10) when a condensate gets born. Indeed, adding the condensate between the two such intersection points yield the same quasi-momentum (and hence leave the filling fraction intact). We have confirm this by numerically solving the Bethe equations at comparably large system sizes, as demonstrated in Fig. 9 (again for the isotropic case). Starting with low filling fractions for both cuts, we observe an apparent "attraction" between the cuts (cf. the second cut connecting (ζ 1 ,ζ 1 ) in panel (a) in Fig. 9). The four branch points in panel (a) in Fig. 9, reading ζ 1 = 0.10884679 + 0.047665716i and ζ 2 = 0.07330641 + 0.04152184i (with filling fractions ν 1 = 0.02, ν 2 = 0.06, = 1 and mode numbers n 1 = 1, n 2 = 2), have been found numerically using the recipe given in Appendix E.2. Increasing the filling fractions of the both cuts, they eventually merge and the end points of a condensate (i.e. a logarithmic branch cut), exemplified in panel (b) in Fig. 9. The four branch points in panel (b) in Fig. 9 are ζ 1 = 0.09587725 + 0.05961115i and ζ 2 = 0.07169871 + 0.04814544i with filling fractions ν 1 = 0.025, ν 2 = 0.075, = 1 and mode numbers n 1 = 1, n 2 = 2. The anisotropic case is plagued by the same numerical difficulties as encountered earlier for the one-cut solution with n ≥ 2, cf. Appendix E.1. Once again, we do not expect there to be any qualitative difference compared to the isotropic model.

Special case: bion
We have discussed earlier that the easy-axis regime of δ > 0 supports a special type of twocut solutions that are not allowed in the other two (that is isotropic and easy-plane) regimes: we are referring here to the bion solution which we have worked out in Section 5.2. Our motivation for investigating this case comes partly due to the curious domain-wall freezing phenomenon in a classical quench in the Landau-Lifshitz field theory [16,17].
To quantise the classical bion solution, we require the reality condition as in the one-cut case. Notice however that bions represent maximally saturated state with the total filling ν = 1 2 . At half filling one commonly finds condensates. For the bion solution there is no loss of generality in fixing the mode numbers of the two cuts to be +1 and −1, namely ∆n = 2. Recall that if the two cuts are well separated, each of them can indeed for a condensate on their own. In the case of bion however, we find the two branch cuts quite close to each other with a condensate shared amongst the two. Inspired by the emergence of a "double condensate" in the isotropic case [11,64] (precisely for the case of a two-cut solution with mode numbers +1 and −1), we conjecture that the very same phenomenon takes place in the present situation of a quantised bion. In particular, a "double condensate" arises when two branch cuts with modes numbers +1 and −1 intersect and thereby form a logarithmic branch cut with of 'doubled' uniform density 2i. Such objects, twice as dense as ordinary Bethe strings, have been observed in several cases of solutions to the isotropic Bethe equations with small system sizes, such as in Ref. [64].
The conjectured bion contours with a double condensate, see Fig. 10, permit to perform semi-classical quantisation also for the bion solution. The red dashed line denotes the "double condensate", i.e. ρ(ζ)(1 + δζ 2 ) = 2i. (6.11) We have checked the filling fraction ν = ν 1 + ν 2 = 1 2 , which coincides with the classical result, as well as the total momentum P = 0 (mod 2π) and total energy E = 3.96045 obtained by numerically integrating along the proposed contours, while for the classical bion P = 0 and E = 3.960358 (6). The results give a strong indication that the identified contours indeed correspond to a quantised bion.
Similar to the soliton limit for the classical bion solution that results in a static kink, cf. Section 5.2, we conjecture that similar limit also exists for the semi-classical limit. We expect that a genuine quantum static kink would form in the limit ζ 1 , ζ 2 → i/ , and the density of Bethe roots is distributed between −i/ and i/ at imaginary axis as double condensate. We need more evidence to confirm the conjecture, which we leave for the future.
It goes without saying that the ultimate check would be to find the corresponding solution to the Bethe equations for an increasingly large system size. This task appears quite challenging numerically, since for a large system size the solutions to the anisotropic Bethe equations (3.3) at half filling are difficult to obtain.
Despite of the aforementioned difficulties, it is possible to look for all of the exact eigenstates for a small system size (typically of order L ∼ 10), see for instance, Ref. [65,66] in the isotropic case. It would be worthwhile finding an exact quantum eigenstate that corresponds to the quantised bion even for a relatively small system size.

Semi-classical norms and overlaps
In this section we describe how to obtain the classical limits of overlaps between Bethe eigenstates. Specifically, our objective is to give closed formulae for (i) the Gaudin norm and (ii) the Slavnov inner product between an on-shell and off-shell Bethe state [67]. There exist two methods to achieve this. The first one, proposed by Gromov et al. in [31], is to start from the explicit determinant Bethe Ansatz formulae and deduce their semi-classical limits by evaluating them for the density resolvent of a particular finite-gap configuration. The other one, developed for the isotropic (XXX) Heisenberg model by Kostov and Bettelheim in Refs. [32][33][34], uses techniques of functional integral and complex analysis. These methods provide the leading (i.e. classical) contributions to the norms and overlaps. In Sections D.1 and D.2 we demonstrate how this leading contribution can be inferred numerically from the finite-size analysis.

Gaudin norm
The method proposed by Kostov in Refs. [32,33] has already been adapted for the specific case of the anisotropic Heisenberg model in [68,69]. Here we first focus on the computation of the Gaudin norm by following a more direct approach of [31] and generalise it for the presence of interaction anisotropy. The core idea of [31] is indeed elementary and makes use of the coarse-graining of the Gaudin determinant. By taking its logarithm, casting it as a Riemann sum and subsequently taking the L → ∞ limit, one is left with a complex integration along the density contours.
The Gaudin norm of a finite-volume Bethe eingenstate, where we have expressed the volume-law coefficient C 1 in terms of the resolvent density ρ(ζ) with support on a union of contours C = ∪ j C j .
We note that the dominant subleading correction to the above formula is quite subtle and has the form log N (L) = C 1 L + 1 2 log L + O(L 0 ) (for C 1 ∈ O(1)), as discussed in Ref. [31]. The numerical verifications (without or with condensate) are presented in Appendix D.1.

Slavnov overlap
An advantage of the functional approach developed in Refs. [32,33] is that, unlike the method of Ref. [31], it does not rely on the clustering properties of Gaudin determinant. Here we quote the final result of [68] (cf. formula (1.5) in there) for the anistropic Landau-Lifshitz field theory where p 1 (ζ) and p 2 (ζ) are classical quasimomenta associated to two semi-classical Bethe eigenstates |{ϑ 1 } and |{ϑ 2 } 10 , respectively. Φ b (z) denotes the Faddeev's quantum dilogarithm function [70], defined through the following integral representation dt t e zt sin(b 2 t) sinh(πt) , (7.4) This function can be though of as a 'quantum deformation' of the ordinary dilogarithm function Li 2 (z) to which it reduces in the isotropic limit δ → 0. 11 The contour prescription in Eq. (7.3) is such that C 1,2 tightly enclose the supports of the respective density resolvents, cf.
Ref. [33]. Further simplification can be made in the semi-classical limit η → /L which implies b → 0. In this limit the quantum dilogarithm simplifies to [71] Li 2 e i p 1 (ζ)+p 2 (ζ) (7.6) Similar to the Gaudin norm, the dominant subleading correction is of form log {ϑ} 1 |{ϑ} 2 = C 2 L + 1 2 log L + O(L 0 ). For the coinciding sets of rapidities, we retrieve the leading order expression for the Gaudin norm, reading This expression is indeed in agreement with Eq. (7.2), which can be established after explicit integration upon identifying the resolvent density with quasi-momentum, iπρ(ζ) = πn − p(ζ), Li 2 e 2ip(ζ) + π 2 ρ 2 (ζ)(1 + δζ 2 ) 2 − π 12 , (7.8) Eq. (7.6) comes with a practical limitation which has been acknowledged previously in Ref. [33], namely, the integration contours C must avoid crossing any branch cut. In the present case the issue is even more severe due to additional logarithmic branch cuts of the dilogarithm function. The limitation makes the numerical verification challenging. Still we manage to check it in the case of overlaps with a vacuum descendant state (see Appendix D.2).

Correlation functions: exact classical-quantum correspondence
As discussed in Section 7, the knowledge of physical branch cuts not only facilitates the semi-classical quantisation of finite-gaps solution but also proves useful in calculating the spectrum, overlaps and norms of semi-classical Bethe eigenstates. In this section, we would like to examine more general observables given by correlation functions of physical observables. The task of computing exact values of correlation functions of physical observables is indeed quite challenging task despite integrability. There have been numerous works on this subject, using form-factor expansion or bootstrap methods, e.g. [39,[72][73][74][75]. Here we are however interested particularly in the semi-classical regime where the known methods do not apply.
We recall that for semi-classical states of the quantum harmonic oscillator, the exact classical-quantum correspondence for the correlation functions can be established analytically, as outlined in the opening Section 2 of the paper. Remarkably, we noticed that correlation functions evaluated on the semi-classical states are in close relation with (phase-space) averaged correlation functions of classical spin fields governed by the Landau-Lifshitz field theory. Based on these observations, we can put forward a conjecture on an exact classical-quantum correspondence for the correlation functions for semi-classical states in generic quantum interacting integrable systems. Such a correspondence has been previously envisioned by Smirnov in his paper on the semi-classical limit of form factors [76].
In order to illustrate the correspondence, we consider correlation functions of quantum one-cut eigenstates with relatively small system sizes. Higher-point correlation functions can be computed directly from coordinate Bethe ansatz. We focus on the static correlation functions, that equal-time one-point ( σ x j and σ z j ) and 2-point ( σ x jσ x k and σ z jσ z k ) functions. Specifically, we consider systems with filling fraction ν = 1 3 , mode number n = 1 and system sizes L = 15 and 18 (with number of down-turned spins being 5 and 6 respectively). We construct the coordinate Bethe ansatz wavefunction in the local basis ofσ z by solving the Bethe equations (3.3). Notice that the quantum wavefunctions are translationally invariant, meaning that one-point functions do not depend on the spatial coordinate (lattice index) and that the two-point functions only depend on the difference of the two coordinates.
In addition, compared the quantum results with the phase-space averaged classical solutions. The classical spin field is clearly non-uniform in space, unlike translationally invariant quantum eigenstate. This alone already suggests that the classical correlation functions require additional averaging over the classical phase space to eliminate spatial dependence before making a comparison between the two, cf. Ref. [76]. For instance, for equal-time correlation functions we expect that where the classical spin field S x (x) has classical period , whereas |{ϑ} is the corresponding semi-classical quantum eigenstate in the limit L → ∞.
In our examples, the classical spin field S(x, t) with n = 1 and ν = 1 3 is with classical period . For illustration we given an explicit example. Let · · · 5 and · · · 6 denote correlation functions with the state with M = 5 down spins (L = 15 total spins) and M = 6 down spins (L = 18 total spins) respectively. For the one-point function, we find  Fig. 11, we observe some discrepancies (which decrease as the system size L increases.) We nonetheless expect that in the limit L → ∞ We believe that these numerical demonstrations are convincing enough to conjecture that in semi-classical states the correlation functions of quantum observables correspond to phasespace averaged classical correlation functions. It would be of great value to find a compact and general proof for this assertion.

Conclusion and outlook
In this work, we have studied the structure of the semi-classical spectrum of the anistoropic Heisenberg spin-1 2 chain with weak interaction anisotropy. Specialising to the easy-axis regime, we outlined how at the classical level the semi-classical eigenstates emerge as non-linear spin waves, governed by the Landau-Lifshitz field theory with uniaxial anisotropy. This identification has been achieved with aid of exact analytical methods of quantum and classical integrability.
Using the asymptotic Bethe ansatz technique we first derived the Riemann-Hilbert problem which prescribes the jump discontinuity conditions for a multi-valued complex function called quasi-momentum. Quasi-momentum encodes the information about the spectrum of nonlinear eigenmodes for a class of finite-gap solutions of the anisotropic Landau-Lifshitz ferromagnet.
We proceeded by outlining the main ingredients of the algebro-geometric integration procedure. We have recast the standard Lax representation in the form of an auxiliary linear transport problem with a flat connection as the adjoint linear equations of motion for squared eigenfunctions. The latter contain the dynamical information of the system. The dynamics can be linearised by the Abel-Jacobi transform on a finite-genus Riemann surface. The corresponding components of the physical spin field can be obtained using so-called reconstruction formulae.
To exemplify the integration procedure, we explicitly worked out the two simplest types of solutions: the one-cut solutions describing precessional motion about the anisotropy axis, and the two-cut finite-gap solutions which take the shape of elliptic waves. We have paid particular attention to the special case of an elliptic solution describing bions and a singular limit thereof which gives rise to the static kink solution.
The Riemann-Hilbert problem for the classical quasi-momentum offers a natural tool for performing semi-classical quantisation of finite-gap solutions, expressed as a singular integral equation for the spectral resolvent. The resolvent is supported on the branch cuts of the finitegap Riemann surface. By following the recipe devised in Ref. [11], we describe an algorithm for numerically obtaining the physical cut (i.e. density contours). Physically speaking, every cut represents a macroscopic condensate of magnons (Bethe roots) which combine in one occupied non-linear mode in the classical spectrum; semi-classical quantisation amounts to dissolve each cut into a large (but finite) number of unit quanta (carried by magnons), according to the resolvent density supported on the cut (which can be achieved with prescribed precision). When densities exceed a critical value, the physical contours undergo certain 'non-perturbative effect', which leads to the formation of special condensates of constant unit density. To properly resolve these situations, one has to take into account quantum corrections.
In this work we concentrate primarily on formal properties of semi-classical eigenstates and various mathematical underpinnings. Our hope is that this can provide a solid foundation for further development and ultimately pave the way to various physics applications, especially in the domain of out-of-equilibrium dynamics. In recent years, there has been a tremendous progress with quantum integrability techniques which, among other, enabled to study latetime relaxation from a highly-excited many-body initial state (widely known nowadays as the 'quantum quench' [3], see also [77][78][79]). In this respect, a functional integral approach called the Quench Action [80,81] has established itself as a central tool. The method relies on the ability to compute exact overlaps and norms which are in general difficult to deal with. Even though general explicit expressions for the overlap between a semi-classical eigenstate with an off-shell state are known due to Ref. [32,33,68], we have not succeeded yet to take full advantage of them. Despite that, we have benchmarked the semi-classical Gaudin norms and Slavnov overlaps on a few simplest finite-gap solutions and found good convergence.
One of the main difficulties that we encountered is that such formulae require a suitable contour prescription to make sure that branch cuts of the integrand are avoided which, in our experience, is not easy to overcome. An open question remains whether there exists an alternative formula which only invokes the resolvent densities, similarly to that for the semi-classical limit of the Gaudin norm [31]. Making these objects available and ready to use would be an important step closer toward the solution of the 'semi-classical' quench problem. For instance, one of the most prominent examples we mentioned is the 'melting of an initial magnetic domain', which exhibits intricate dynamical phases depending on the value of interaction anisotropy both in the classical and quantum cases. While a closed-form solution in the classical field theory has been recently obtained in Ref. [16] using the inverse scattering transform, obtaining an analytic solution to its quantum counterpart remains an open challenge. Given that the problem is inherently semi-classical in nature, solving the quench dynamics at the level of the semi-classical spectrum seems a natural route to establish the classical-quantum correspondence as discussed in Ref. [16] and further corroborated in Ref. [17].
Lastly, aiming at uncovering a general classical-quantum correspondence directly at the level of correlation functions [76], we briefly investigated the structure of correlation functions in semi-classically quantised finite-gap solutions and compare them to the classical counterparts. In this regard, the key question is whether quantum observables (including correlation functions) averaged on the semi-classical eigenstates simplify in any manner and how they are related to classical phase-space averages for the class of finite-gap solutions. Following an analogy from the harmonic oscillator, we formulated a conjecture for how static (i.e. equaltime) correlation function of local spin operators evaluated in semi-classical eigenstates show up as phase-space averages of classical spin-field configurations of the finite-gap spectrum, which we supplemented with non-trivial numerical checks. Our finding are aligned with the earlier results by Smirnov on the semi-classical limit of the form factors [76]. Providing a precise statement with a rigorous proof would, in our opinion, be highly valuable. Another potentially fruitful approach to explore is to consider the quantum correlation functions in terms of the quantum separated variables [52,82,83]. Comparing to the classical separated variables [36], it is quite opaque at the moment whether quantum dynamics can be described by the quantum analogue of the "angle variables". We believe that quantum integrable lattice models and spin chains provide paradigmatic examples to address these aspects.

Appendices A Riemann-Hilbert problem in ζ-plane
The Riemann-Hilbert problem can be formulated in terms of spectral parameter ζ = 1 µ , from which the phenomenon of condensate is more intuitive. For instance, where C j is the j-th branch cut in ζ plane, and quasi-momentum p(ζ) is defined as where the new integration kernel becomes Accordingly, the density (of the Bethe roots) becomes where the orientation of integration along the branch cut C j is inverted, i.e. from the branch point with negative imaginary part to the one with positive imaginary part.

B Finite size corrections to Riemann-Hilbert problem
As for the logarithm of Q , it is trickier to take the same semi-classical limit, and we need to split the term into the anomalous part and normal part [54,84], i.e. log Q where the cut-off parameter K satisfies the following property, and We denote the anomalous part as while the normal part is For the normal part, we can do this same expansion as in Eq. (3.12), i.e.
Combining two parts, we obtain Meanwhile, for the anomalous part, denoting m = k − j, we have .

(B.9)
We can expand Lµ j+m as where the "constants" can be expressed in terms of density function ρ(µ), i.e. and By combining the mth and −mth terms in the sum, we express the leading order of the sum as , (B.13) using the identity below (B.14) The first part can be combined with the sum for |m| > K, since Taking the limit K → ∞ (beware that K/L → 0), for the second part we have And putting back the values in Eq. (B.11), we will obtain the finite-size correction in Eq. (3.18). In addition, in terms of ζ variable, the finite-size correction takes the form of

C Useful formulae for elliptic functions
We define several useful functions and formulae that are extensively referred in Section 5.2.
To begin with, we define the elliptic integral of the first kind as , (C.1) and the Jacobi elliptic function sn(x, k 2 ) is defined as the inverse of the elliptic integral of the first kind, i.e.
and without ambiguity, we denote sn(x, k 2 ) =: sn(x). And similarly, we could define other Jacobi elliptic functions, , such that sn 2 x + cn 2 x = 1, k 2 sn 2 x + dn 2 x = 1. (C.5) When shifting the argument by quarter period K(k 2 ), we have In addition, we also make use of theta functions to express the spin field. The most important one here is More detailed explanations and properties of elliptic functions can be found in Ref. [85,86].

D Numerical tests D.1 Gaudin norm
We have performed numerical checks for two scenarios: firstly, the Gaudin norm of a 1cut solution without condensate, and secondly, the one with condensate, both in the regime δ = 0 and δ > 0. The isotropic case without condensate has been studied in Ref. [31], providing us a possible benchmark. The second check is more interesting, since it provides us a quantitative and non-trivial confirmation of our proposal for the location of the condensate and the existence of additional contours in Section 6.2.1.
We have calculated the Gaudin norm for 1-cut solutions with mode number n = 1 and filling fraction ν = 0.1, in either isotropic (δ = 0) and anisotropic (δ = 1) cases. There is no condensate in both cases with the parameter considered. The linearly fitted numerical results from the finite-size analysis are Comparing to the results of functional integral (7.2) (denoted as C 1 in the table below), we obtain the following table, As we can conclude from the results, the accuracy of the functional approach is very high in both isotropic and anisotropic cases, providing a demonstration on the powerfulness of the functional approach that only requires the knowledge of the density along the physical cuts.
In addition, we scrutinise the Gaudin norm of 1-cut solution with mode number n = 1 and filling fraction ν = 1 3 , in either isotropic (δ = 0) and anisotropic (δ = 1) cases, and in  both cases, there exist condensates. In the presence of condensate, the functional results can be obtained from computing the integral (7.2) along 3 contours i.e. C 1 , the parts that coincide with the original classical contour with density ρ 1 (ζ) in Eq. 6.7 (green dashed line in panel (b) of Fig. 7), C 2 , the additional contour with density ρ 2 (ζ) in Eq. 6.7 (yellow dashed line in panel (b) of Fig. 7) and C cond , the condensate part with density ρ cond = i 1+δζ 2 . The linearly fitted numerical results from the finite-size analysis are Comparing to the results of functional integral (7.2) (denoted as C 1 functional in the table below), we obtain the following table, 0.091273(6) 0.091121(9) δ = 1 0.082597(1) 0.081761 (2) Similarly to the previous cases, the functional results remain highly accurate. Above all, it provides an evidence of the location of condensate and the existence of addition contour(s). Indeed, if one would try a different contour, e.g. the classical contour without condensate, the functional result differs hugely from the results above. Before using Eq. (7.6), we need to calculate the "quasi-momentum" for the vacuum descendant state |φ . In terms of algebraic Bethe ansatz, the vacuum descendant state without inhomogeneities can be considered as [87] |φ = | ↓ · · · ↓ M ↑ · · · ↑ = lim ξ j →0 M j=1 B(ξ j )|0 , |0 = | ↑ · · · ↑ , (D. 13) where |0 is the pseudo-vacuum state | ↑ · · · ↑ , and B(λ) is the off-diagonal part of the transfer matrix. Therefore, we can express the density of the "off-shell Bethe roots" as ρ φ (ζ) = ν 1 δ(ζ), ξ 1 , · · · ξ M → 0. (D.14) Without losing generality, we set the classical period = 1 for both classical limit of the states. And for the quasi-momentum of |φ , we have Hence, we can compute the results of functional integral (7.6) with the contours chosen as panel (a) in Fig. 14 13 .

E Numerical recipe E.1 Numerical recipe to solve Bethe equations
In order to numerically determine the Bethe roots for the Bethe equations (3.3) with finite but large system size L that correspond to the quantised 1-cut states in the thermodynamic limit, we first use the algorithm in Section 7 of Ref. [11]. This will give us the location of the Bethe roots that are solution to the rational Bethe equations (i.e. at isotropic point).
In addition, we use the algorithm in Appendix C of Ref. [68], where we start to use the solution of the isotropic case as the starting point of Newton-Raphson method and turn on the anisotropy bit by bit to solve the Bethe equations. It works quite well for the case of mode number 1. But the current numerical method fails to get good convergence for mode number n > 2 and solutions with more than 2 cuts in the classical limit.

E.2 Determining branch points from filling fractions and mode numbers
We propose a numerical recipe on how to determine the multiple-cut branch points based on the mode numbers and filling fractions. For the sake of simplicity, we demonstrate the procedure with 2-cut solutions.
Firstly, only the ratio of two sets of branch points (µ 1 ,μ 1 , µ 2 ,μ 2 ) is related to the mode numbers and filling fractions, if we do not require a fixed classical period . This can be understood by a simple counting of degrees of freedom. To begin with, we can give any nonzero integer number as the mode number of the first cut. And there are 4 real parameters 13 The reason for this choice of integration contours is to avoid the branch cuts derived from the dilogarithm function. More details are explained in [33]. remaining, namely, Re µ 1 , Im µ 1 , Re µ 2 and Im µ 2 . However, there are only 3 parameters to be fixed, i.e. the mode number for the other cut, and 2 filling fractions. Of course, the modules of each cuts are related to the classical period in a nonlinear way. However, it is always possible to divide the branch points by the classical period to compare to the quantum results, once the mode numbers and filling fractions are fixed.
In addition, we would like to remark that unlike the 1-cut case, the determination is highly nonlinear, related to elliptic functions and integrals. Hence, we do not expect a simple analytic closed-form formula, and we are going to use the following numerical procedure: • We set the real part of branch points of the first cut to be 1, i.e. Re µ 1 = 1.
• For a range of Imµ 1 ∈ (b 1 , b 2 ) and Re µ 2 ∈ (c 1 , c 2 ), calculate Im µ 2 that satisfies the designated mode number n 1 and n 2 . The ranges are chosen through estimation. More precisely, with any Im µ 2 , we can obtain the form of dp requiring a vanishing A-cycle. And we have = 2πn 1 B 1 dp/ , (E.1) and we can numerical fix Im µ 2 satisfying B 2 dp = 2πn 2 . (E.2) • Hence we can find out the following quantities, and their dependence on Im µ 1 and Re µ 2 .
• By requiring thatν 1 = ν 1 , we could obtain a "curve" in the Imµ 1 , Re µ 2 -plane. And by further requiring thatν 1 = ν 2 , we are left with one point in the Im µ 1 , Re µ 2 -plane, and let us denote the point (b, c). With the information before, we also know the corresponding Im µ 2 , denoted as d here.
• From the information above, we get the following branch points µ 1 = 1 + bi, µ 2 = c + di and their complex conjugates, with classical period , which is the 2-cut solution with mode number n 1 , n 2 and filling fraction ν 1 , ν 2 .