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

Yuan Miao1, Enej Ilievski1,2 and Oleksandr Gamayun1

1 Institute for Theoretical Physics, University of Amsterdam,

Science Park 904, 1098XH Amsterdam, the Netherlands

2 Faculty of Mathematics and Physics, University of Ljubljana,

Jadranska 19, 1000 Ljubljana, Slovenia

<EMAIL_ADDRESS>
## Abstract

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 ground state. Special emphasis is devoted to the simplest types
of solutions, describing precessional motion and elliptic magnetisation waves.
The internal magnon structure of classical spin waves is resolved by
performing the semi-classical quantisation using the Riemann–Hilbert problem
approach. We present an expression for the overlap of two semi-classical
eigenstates and discuss how correlation functions at the semi-classical level
arise from classical phase-space averaging.

###### Contents

  1. 1 Introduction
  2. 2 Harmonic oscillator: introducing the spectral curve
    1. 2.1 Semi-classical limit
    2. 2.2 Classical limit of quantum correlation functions
  3. 3 Asymptotic Bethe Ansatz: the XXZ Heisenberg model
    1. 3.1 Anisotropic Landau–Lifshitz field theory
  4. 4 Finite-gap integration method
    1. 4.1 Auxiliary linear problem
    2. 4.2 Finite-gap solutions
      1. 4.2.1 Dynamical divisor
  5. 5 Examples of finite-gap solutions
    1. 5.1 One-cut rational solution
      1. 5.1.1 One-cut solution from asymptotic Bethe ansatz
    2. 5.2 Bions from degenerate two-cut solutions
  6. 6 Semi-classical quantisation
    1. 6.1 Physical contours
    2. 6.2 Formation of condensates
      1. 6.2.1 One-cut case with condensate
      2. 6.2.2 Multiple cuts
    3. 6.3 Special case: bion
  7. 7 Semi-classical norms and overlaps
    1. 7.1 Gaudin norm
    2. 7.2 Slavnov overlap
  8. 8 Correlation functions
  9. 9 Conclusion and outlook
  10. A Riemann-Hilbert problem in $\zeta$-plane
  11. B Finite size corrections to Riemann-Hilbert problem
  12. C Useful formulae for elliptic functions
  13. D Numerical tests
    1. D.1 Gaudin norm
    2. D.2 Slavnov overlap
  14. E Numerical recipes
    1. E.1 Numerical solution to Bethe equations
    2. E.2 Determining branch points from filling fractions and mode numbers

##  1 Introduction

Thermodynamic systems of interacting quantum particles present an outstanding
challenge theoretical physics. In spite of their inherent complexity,
tremendous progress has been made recently in understanding various facets of
quantum many-body systems, including thermalisation [1, 2], far-from-
equilibrium dynamics [3, 4], quantum transport [5, 6] and entanglement
dynamics [7], especially after the inception of generalised hydrodynamics [8,
9]. In this regard, dimensionally constrained models played a prominent role
as they not only permit for moderately efficient numerical simulations but
also provide a playground for developing and testing analytical approaches.
Integrable models are of special interest as they enable us to obtain non-
perturbative closed-form solutions that are otherwise rarely available.

This paper is devoted to study the structure and properties of semi-classical
eigenstates and the emergence of classical dynamics in the quantum Heisenberg
chain, an archetypal example of a quantum many-body system. Emergent
integrability in the $\mathcal{N}=4$ superconformal Yang–Mills theory [10, 11]
generated some amount of interest in the semi-classical limits in integrable
quantum spin chains. Particularly, the complete semi-classical spectrum of
isotropic Heisenberg model has been first described in [12] and subsequently
studied in great detail in [13]. The Heisenberg model however appears in a
variety of physics applications in the domain of statistical mechanics and is
of particular relevance for studying basic principles of out-of-equilibrium
many-body dynamics. Following the spirit of Ref. [13], we devote this work to
study the semi-classical eigenstates in the easy-axis regime of the
anisotropic Heisenberg spin-$1/2$ chain.

A separate motivation for studying the semi-classical part of the spectrum
originates from recent interest in magnetisation transport in integrable
quantum chains, where the anisotropic Heisenberg chain plays a prominent role.
Several conspicuous similarities with a purely classical magnetisation
dynamics governed by the Landau–Lifshitz ferromagnet [14] have been found,
both at qualitative and quantitative levels, which firmly point towards a
particular type of a classical–quantum correspondence [15, 16, 17]. On the
classical side, a key piece of evidence rests on an exact solution to the
initial value problem with a domain-wall initial profile [16] which discerns
three different dynamical regimes depending on the value of anisotropy: (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. This matches the phenomenology of the
quantum isotropic Heisenberg spin-$1/2$ chain inferred previously in [15]. In
this paper we specialise to the easy-axis regime where absence of transport
has been linked with the presence of stable kink in the spectrum [16]. Despite
that, in what precise manner do such kinks arise as coherent superposition of
magnon excitations of the underlying quantum chain remains unknown. The task
at hand is to perform semi-classical quantisation of classical nonlinear spin
waves that arise in the weakly-anisotropic ferromagnetic Heisenberg chain. The
future hope is to study non-equilibrium dynamics directly at the level of
semi-classical eigenstates and thus put the classical–quantum correspondence
on a firm footing.

It deserves to be emphasised that the aforementioned classical–quantum
correspondence of spin transport is not a particularity of the domain wall
physics but likewise manifests itself in thermal equilibrium states (i.e. at
finite density of magnon excitations). There however it comprises different
types of dynamics regimes. Most prominently, in the _isotropic_ Heisenberg
quantum chain, finite-temperature spin transport (in the zero magnetisation
section) is superdiffusive [18, 19, 20], belonging to the KPZ universality
class (characterised by dynamical exponent $z=3/2$ [21, 22]). Such an
anomalous behaviour has been attributed to interacting (and thermally dressed)
‘giant magnons’, referring to semi-classical eigenstates which at the
classical level show up as soliton modes. In contrast, genuine quantum
excitations associated to bound magnons states (Bethe strings) become
suppressed in this regime [23]. Curiously, this anomalous feature nevertheless
entirely disappears upon introducing any amount of interaction anisotropy: on
the easy-plane side, one finds ballistic transport characterised by a finite
spin Drude weight [24, 25], whereas easy-axis anisotropy restores normal
diffusion [26, 27, 21]. In our perspective, these findings offer an extra
motivation to carefully examine the structure of semi-classical eigenstates in
integrable quantum lattice models.

In this work we use the asymptotic Bethe ansatz approach [28] to identify and
describe the semi-classical part of the spectrum in the Heisenberg spin-$1/2$
chain with uniaxial anisotropy. We shall in large part follow the footsteps of
Refs. [12, 13] by employing an algebro-geometric integration technique [29,
30], the Riemann–Hilbert formalism. The main object of interest is the so-
called classical spectral curve which provides complete information about the
classical finite-gap spectrum of the anisotropic Landau–Lifshitz ferromagnet.
Our aim is to make the exposition self-contained and pedagogical. Our
attention will be largely devoted to the emergence of special types of semi-
classical solutions that represent bions and kinks, for which anisotropy
proves crucial 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 begin by briefly
introducing the notion of the spectral curve, giving the most basic example of
the harmonic oscillator. Next, in Sec. 3 we outline the asymptotic Bethe
ansatz technique by applying it to the anisotropic Heisenberg chain. We
proceed by solving the classical finite-gap solutions and explicitly working
out two specific examples in Sec. 4 and 5. The semi-classical quantisation is
then carried out in Sec. 6 and calculations of semi-classical overlaps are
given in Sec. 7. Lastly, in Sec. 8 we formulate a conjecture for a
classical–quantum correspondence of correlation functions. We finish in Sec. 9
with a conclusion and an outlook.

##  2 Harmonic oscillator: introducing the spectral curve

Prior to delving deep into the realm of many-body systems, we would like to
first familiarise the reader with several technical tools that constitute the
foundations of the algebro-geometric method to solve differential equations.
To this end, we shall describe the semi-classical spectrum of a single
quantum-mechanical degree of freedom, the good old quantum harmonic
oscillator,

$\hat{H}=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}\hat{x}^{2}.$ (2.1)

As usual, we use a canonical pair of position and momentum operators,
satisfying

$[\hat{x},\hat{p}]={\rm i}\hbar.$ (2.2)

In the language of second quantisation, the Hamiltonian of the quantum
harmonic oscillator takes the diagonal form

$\hat{H}=\hbar\omega\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right),$ (2.3)

in terms of the bosonic annihilation operator

$\hat{a}=\frac{1}{\sqrt{2\hbar}}\left(\sqrt{m\omega}\,\hat{x}+\frac{{\rm
i}}{\sqrt{m\omega}}\,\hat{p}\right).$ (2.4)

The ground state $\mathinner{|{0}\rangle}$ is the Fock vacuum and satisfies
$\hat{a}|0\rangle=0$. Excited eigenstates $\mathinner{|{n}\rangle}$, which
obey
$\hat{a}^{\dagger}\hat{a}\mathinner{|{n}\rangle}=n\mathinner{|{n}\rangle}$,
are produced by iterative application of the creation operator
$\hat{a}^{{\dagger}}$ on the ground state $\mathinner{|{0}\rangle}$,

$|n\rangle=\frac{\big{[}\hat{a}^{\dagger}\big{]}^{n}}{\sqrt{n!}}|0\rangle,$
(2.5)

such that

$\hat{H}|n\rangle=E_{n}|n\rangle=\hbar\omega\left(n+\frac{1}{2}\right)|n\rangle.$
(2.6)

Using a dimensionless coordinate $u=\sqrt{\frac{m\omega}{\hbar}}x$,
eigenfunctions $\psi(u)=\langle u|n\rangle$ can be expressed in terms of their
normalised logarithmic derivative called quasi-momentum (see e.g. [35]),

$\mathfrak{p}(u)=\frac{\sqrt{\hbar
m\omega}}{2}\,\frac{\partial_{u}\psi(u)}{\psi(u)}.$ (2.7)

suppressing the subscript $n$ mostly for the sake of clarity. The Schrödinger
equation for $\psi(u)$ accordingly transforms into a Ricatti equation

$\mathfrak{p}^{2}(u)-{\rm
i}\sqrt{\hbar}\partial_{u}\mathfrak{p}=2m(E-V(u)),\qquad
V(u)=\frac{\hbar\omega}{2}u^{2},$ (2.8)

whereas nodes (i.e. zeros) of the excited wavefunctions
$\mathinner{|{n}\rangle}$ have now turned into simple poles. After a short
exercise, one can deduce the following representation [35]

$\mathfrak{p}(u)={\rm
i}\sqrt{m\hbar\omega}\left(u-\sum_{j=1}^{n}\frac{1}{u-u_{j}}\right),$ (2.9)

with poles $u_{j}$ of the quasi-momentum satisfying a simple system of
equations

$u_{j}=\sum_{k\neq j}\frac{1}{u_{j}-u_{k}}.$ (2.10)

In this description, every eigenstate $\mathinner{|{n}\rangle}$ gets assigned
a unique set of poles $u_{j}$, with $j=1,\ldots,n$ and accordingly Eqs. (2.10)
bear a direct analogy to the celebrated Bethe ansatz equations arising in
integrable _interacting_ quantum systems. By maintaining this analogy, it is
furthermore convenient to introduce the $Q$-polynomial

$Q_{n}(u)=\prod_{j=1}^{n}(u-u_{j}),$ (2.11)

which in integrable quantum spin chains corresponds to the Baxter’s
$Q$-function. By integrating the quasi-momentum we can readily retrieve the
wavefunction for the $n$th excited state,

$\psi_{n}(u)=\frac{1}{\sqrt{\mathcal{N}}}\,e^{-u^{2}/2}Q_{n}(u),$ (2.12)

with energy $E_{n}=\left(n+\frac{1}{2}\right)\hbar\omega$ and normalisation
$\mathcal{N}$. The solutions to the Bethe-like equations (2.10) are none other
than Hermite polynomials $H_{n}(u)$ of order $n$, that is

$Q_{n}(u)=H_{n}(u).$ (2.13)

Poles $u_{j}$ are thus identified with zeros (roots) of Hermite polynomials.

###  2.1 Semi-classical limit

We next analyse the semi-classical eigenstates of the quantum harmonic
oscillator. Prior to that, we shall briefly remind on some of the well-known
results of the classical harmonic oscillator

$\mathcal{H}=\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}x^{2}.$ (2.14)

Here position $x$ and momentum $p$ are phase-space coordinates with the
canonical Poisson bracket

$\\{p,x\\}=1.$ (2.15)

Solutions to the classical harmonic oscillator are conventionally expressed in
terms of trigonometric functions

$p=\sqrt{2mE}\cos(\omega t+\phi),\qquad
x=\frac{1}{\omega}\sqrt{\frac{2E}{m}}\sin(\omega t+\phi),$ (2.16)

where $E$ is the value of energy and offset angle $\phi$ is determined by the
initial condition.

The classical harmonic oscillator is possibly the simplest example of a
dynamical system that is integrable in the Liouville–Arnol’d sense. The
associated action-angle pair of variables is simply

${\rm S}=\frac{E}{\omega},\qquad\varphi=\omega\,t,$ (2.17)

satisfying the canonical Poisson relation $\\{{\rm S},\varphi\\}=1$.

##### Lax representation.

We proceed by outlining an algebraic reformulation of the above construction.
The notion of algebraic integrability rests on the concept of the Lax
representation [36]. Here we shortly describe how this construction works by
working it out for the toy model – the classical harmonic oscillator. To begin
with, we emphasise that the Lax representation for an integrable dynamical
system is not unique. In the present case, a suitable choice for the _Lax
pair_ of $2\times 2$ matrices $L$ and $M$ is as follows,

$L(v;t)=\begin{pmatrix}p+{\rm i}m\omega xv&m\omega x-{\rm i}pv\\\ m\omega
x-{\rm i}pv&-p-{\rm i}m\omega xv\end{pmatrix},\qquad
M=\frac{\omega}{2}\begin{pmatrix}0&-1\\\ 1&0\end{pmatrix}.$ (2.18)

Lax matrix $L(v;t)$ depends on time $t$ through $x(t)$ and $p(t)$. In
addition, there is _analytic_ dependence on the so-called spectral parameter
$v\in\mathbb{C}$ which, as we clarify in a moment, is of pivotal importance
for algebraic integrability.

The equation of motion for the classical harmonic oscillator is equivalent to
the following equation of motion of the Lax matrix $L(v;t)$,

$\frac{\mathrm{d}}{\mathrm{d}t}L(v;t)=\left[M,L(v;t)\right].$ (2.19)

Although at first glance it may appear that we have not gained much at all, in
the Lax formulation one can immediately recognise the fact that the
characteristic polynomial of the Lax matrix,

${\rm det}(w-L(v))=0,$ (2.20)

is a time independent quantity. Formally speaking, it defines a complex curve
$\Sigma(w,v)\subset\mathbb{C}$ called the spectral curve. In the current case
it is simply an algebraic curve of the form [36]

$\Sigma:\qquad w^{2}(v)=2mE\,(1-v^{2}).$ (2.21)

Eigenvalues $w_{\pm}(v)$ can be conveniently parametrised in terms of a single
double-valued complex function

$\mathfrak{p}_{\rm cl}(v)=\sqrt{2mE\,(1-v^{2})},$ (2.22)

called the _classical_ quasi-momentum, satisfying $2\cos\mathfrak{p}_{\rm
cl}(v)=\mathrm{tr}\left(\exp L(v)\right)$. The quasi-momentum
$\mathfrak{p}_{\rm cl}(v)$ features a pair of square-root branch points at
$\pm 1$. One can thus interpret the spectral curve as a two-sheet Riemann
surface with a branch cut along the interval $\mathcal{I}\equiv[-1,1]$ and a
pole at infinity, $v_{\rm p}=\infty$. By analyticity, one furthermore has
$\oint_{\mathcal{I}}{\rm d}\mathfrak{p}_{\rm cl}(v)=0$.

##### Semi-classical limit.

We now consider the quantum harmonic oscillator and examine the highly-excited
eigenstates whose nodes densely distribute along the real axis. Our purpose
here is to demonstrate how the classical spectral curve emerges out of a
condensate of zeros of the $Q$-polynomial (2.12). This indeed corresponds to
the semi-classical limit of large mode numbers $n\to\infty$, with $\hbar\sim
1/n$. From the viewpoint of the classical system, the resulting condensate
shows up as a square-root branch cut of the classical spectral curve (2.21).
One can likewise think of the reverse process in which the branch cut
disintegrates into a collection of simple poles, which is none other than the
familiar Bohr–Sommerfeld quantisation rule

$n\gg 1:\qquad\frac{1}{2\pi}\oint_{\mathcal{C}}{\rm
d}v\,\mathfrak{p}(v)=\hbar\,n.$ (2.23)

The (quantum) quasi-momentum, upon substituting $u=\sqrt{2n}v$ and
subsequently taking the limit $n\to\infty$, becomes

$\begin{split}\lim_{n\to\infty}\mathfrak{p}(v)&=\lim_{n\to\infty}{\rm
i}\sqrt{m\hbar\omega}\left[\sqrt{2n}v-\sum_{j=1}^{n}\frac{1}{\sqrt{2n}v-u_{j}}\right]\\\
&={\rm
i}\sqrt{2mn\hbar\omega}\left[v-\int_{-1}^{1}dy\frac{\rho(y)}{v-y}\right],\end{split}$
(2.24)

where we have introduced the distribution of zeros of Hermite polynomial

$\rho(y)=\lim_{n\to\infty}\frac{1}{n}\sum_{j=1}^{n}\delta\left(y-\frac{u_{j}}{\sqrt{2n}}\right).$
(2.25)

We have thus retrieved the classical quasi-momentum

$\mathfrak{p}_{\rm cl}(v)={\rm
i}\lim_{n\to\infty}\sqrt{2mn\hbar\omega}\sqrt{v^{2}-1}=\sqrt{2mE(1-v^{2})},$
(2.26)

characterised by the classical energy

$E=\lim_{n\to\infty}\lim_{\hbar\to E/(\omega n)}n\hbar\omega.$ (2.27)

(a)

(b)

(c)

Figure 1: (a) Normalised wavefunction $\psi_{20}(v)$ with mode number $20$ and
rescaled coordinate $v=u/\sqrt{40}$. Its zeros are marked by red crosses. At
large mode numbers, zeros of wavefunctions approach each other and eventually
condense. The resulting condensate can be viewed as a square-root branch cut
of a spectral curve, cf. Eq. (2.21). Classical spectral curve of the harmonic
oscillator in the $v$-plane (b) and $z$-plane (c). The branch cut in (b) gets
mapped into two punctures of the Riemann surface. The dashed line indicates
the motion of dynamical variable $z^{\star}$.

##### Canonical angle.

The spectral curve can be associated a dynamical angle-type variable
$\varphi$. To find $\varphi$, one requires the solution to the auxiliary
linear problem. The latter takes the form of the Lax equations

$L\,\psi_{\pm}(v;t)=\lambda_{\pm}(v)\psi_{\pm}(v;t),\qquad\frac{\mathrm{d}}{\mathrm{d}t}\psi_{\pm}(v;t)=M\,\psi_{\pm}(v;t).$
(2.28)

The standard recipe is to identify the canonical angle variables with
dynamical poles of the appropriately normalised eigenvectors $\psi_{\pm}$,
normally known in the literature as the Baker–Akhiezer vectors [36]. Below we
employ a slightly different (but equivalent) approach using the squared
eigenfunction, mainly to circumvent the normalisation ambiguity inherent to
the Baker–Akhiezer vectors.

Introducing a $2\times 2$ matrix of eigenvectors $\psi_{\pm}(v,t)$,

$\boldsymbol{\psi}(v;t)=\big{(}\psi_{+}(v;t),\psi_{-}(v;t)\big{)},$ (2.29)

we defined the ‘squared eigenfunction’ as

$\boldsymbol{\Psi}=\boldsymbol{\psi}\,\sigma^{z}\,\boldsymbol{\psi}^{-1},$
(2.30)

satisfying $\det\boldsymbol{\Psi}=-1$. Notice that $\mathbf{\Psi}$ differs
from $L$ by the time-independent normalisation. Matrix $\boldsymbol{\Psi}$
itself thus also evolves according to the Lax equation of motion,

$\frac{{\rm d}}{{\rm
d}t}\boldsymbol{\Psi}(v;t)=\left[M,\boldsymbol{\Psi}(v;t)\right].$ (2.31)

In the present example of the classical harmonic oscillator, the solution to
the above differential equation admits the following explicit form

$\boldsymbol{\Psi}(v;t)=\frac{1}{\sqrt{2mE(1-v^{2})}}\begin{pmatrix}{\rm
i}m\omega x\left(v-\frac{{\rm i}p}{m\omega x}\right)&-{\rm
i}p\left(v+\frac{{\rm i}m\omega x}{p}\right)\\\ -{\rm
i}p\left(v+\frac{im\omega x}{p}\right)&-{\rm i}m\omega x\left(v-\frac{{\rm
i}p}{m\omega x}\right)\end{pmatrix},$ (2.32)

where $x=x(t)$ and $p=p(t)$.

Based on the general rule, the dynamical variables $\gamma_{j}$ correspond to
zeros of the off-diagonal element of the squared eigenfunction
$\boldsymbol{\Psi}$. In the language of algebraic geometry, the full set of
dynamical variables $\\{\gamma_{j}\\}$ constitutes the so-called dynamical
divisor of a Riemann surface. Their equations of motion on a surface are
governed by a system of differential equations that go under the name of
Dubrovin equations. In the case of hyperelliptic algebraic curves of genus
$\mathfrak{g}$, the total number of dynamical variables equals
$\mathfrak{g}+1$. Our toy example involves a single degree of freedom and
therefore we deal with a surface of genus $\mathfrak{g}=0$ (i.e. a Riemann
sphere). Accordingly, there is a single dynamical variable
$\gamma_{1}=\gamma_{1}(t)$, satisfying a simple evolution law

$\gamma_{1}(t)=-i\tan\omega t=-\frac{{\rm i}p}{m\omega x}.$ (2.33)

To explicitly see this, we perform the following variable transformation,

$v\mapsto z(v):\qquad v=-\frac{z-1/z}{z+1/z},$ (2.34)

which resolves the square-root type singularities at $v^{\pm}_{\star}$ and
renders $\lambda^{2}$ a rational function of $z$,

$\Sigma:\qquad\lambda^{2}(z)=2mE\frac{4}{\left(z+1/z\right)^{2}}.$ (2.35)

The original square-root branch cut with branch points $v^{\pm}_{\star}=\pm
1$, see Eq. (2.21), have turned into two regular points on the Riemann surface
located at $z_{\star}\in\\{0,\infty\\}$, whereas the original pole at infinity
has two pre-images at two punctures $z^{\pm}_{\rm p}=\pm{\rm i}$, one on each
Riemann sheet. The dynamical variable $z^{\ast}(t)=\exp{(\varphi(t))}$
satisfies has periodic motion with a linearly-evolving angle variable
$\varphi(t)=\omega\,t$, as shown in Fig. 1.

###  2.2 Classical limit of quantum correlation functions

We finally examine the expectation values of quantum observables computed in
semi-classical eigenstates. We explain how the latter manifest themselves as
classical quantities, thereby establishing an exact classical–quantum
correspondence at the level of correlation functions. In the semi-classical
limit, the average of a quantum observable $\hat{O}$ should be identified with
ergodic phase-space averages

$\lim_{n\to\infty,\hbar\to\frac{E}{\omega n}}\langle
n|\hat{O}(t_{1},\cdots,t_{m})|n\rangle=\frac{1}{T}\int_{0}^{T}{\rm
d}t\,O(t+t_{1},\cdots,t+t_{m}),$ (2.36)

where $T=2\pi/\omega$ denotes the fundamental period.

The correspondence can be readily exemplified on a few particular examples. We
first consider, as a simple example, operator $\hat{O}=\hat{x}^{2}(0)$.
Expectation values on a normalised quantum eigenstates read

$\langle
n|\hat{x}^{2}(0)|n\rangle=\frac{\hbar}{m\omega}\left(n+\frac{1}{2}\right),$
(2.37)

which in the semi-classical limit yields

$\lim_{n\to\infty,\hbar\to\frac{E}{\omega n}}\langle
n|\hat{x}^{2}(0)|n\rangle=\frac{E}{m\omega^{2}}=\frac{\omega}{2\pi}\int_{0}^{\frac{\omega}{2\pi}}{\rm
d}t\,x^{2}(t).$ (2.38)

The outlined correspondence continues to hold even when $\hat{O}$ consists of
several non-commuting operators. To illustrate this, we consider two
observables $\hat{O}_{1}=\hat{x}(0)\hat{p}(t_{0})$ and
$\hat{O}_{2}=\hat{p}(t_{0})\hat{x}(0)$, whose quantum averages yield

$\langle n|\hat{O}_{1}|n\rangle=\langle
n|\hat{x}(0)\hat{p}(t_{0})|n\rangle=\frac{i\hbar}{2}\left[e^{i\omega
t_{0}}(n+1)-e^{-i\omega t_{0}}n\right],$ (2.39) $\langle
n|\hat{O}_{2}|n\rangle=\langle
n|\hat{p}(t_{0})\hat{x}(0)|n\rangle=\frac{i\hbar}{2}\left[e^{i\omega
t_{0}}n-e^{-i\omega t_{0}}(n+1)\right],$ (2.40)

Taking again the semi-classical limit, we arrive at

$\begin{split}&\lim_{n\to\infty,\hbar\to\frac{E}{\omega n}}\langle
n|\hat{O}_{1}|n\rangle=\lim_{n\to\infty,\hbar\to\frac{E}{\omega n}}\langle
n|\hat{O}_{2}|n\rangle\\\ &=-\frac{E}{\omega}\sin(\omega
t_{0})=\frac{\omega}{2\pi}\int_{0}^{\frac{\omega}{2\pi}}{\rm
d}t\,x(t)p(t+t_{0}).\end{split}$ (2.41)

The upshot of this rather elementary calculation is that expectation values of
quantum observables in the semi-classical limit of large mode numbers can be
effectively replaced by the corresponding classical observables. There is no
ordering ambiguity at the leading order, which only enter in the form of
‘quantum corrections’, i.e. subleading contributions
($\sim\mathcal{O}(\hbar)$) to the correlation functions.

This concludes our pedagogical introduction to the basic notions of algebro-
geometric methods. In the remainder of this article we shall employ the same
methodology on a genuine many-body interacting model – the quantum Heisenberg
spin-$1/2$ chain with anisotropic interaction (the XXZ model). Due to multiple
degrees of freedom involved, analytical treatment becomes much more
challenging and how correlators simplify in the classical limit is no longer a
trivial task. Before we come back to this issue in Sec. 8 we first introduce
the formalism and other key ingredients.

##  3 Asymptotic Bethe Ansatz: the XXZ Heisenberg model

The anisotropic (XXZ) quantum spin-$1/2$ chain is often considered as an
archetypal model of many-body quantum dynamics. Besides that, it is one of the
best studied solvable models by Bethe ansatz [37, 38, 39]

$\hat{H}=-J\left[\sum_{j=1}^{L}\hat{S}^{\rm x}_{j}\hat{S}^{\rm
x}_{j+1}+\hat{S}^{\rm y}_{j}\hat{S}^{\rm y}_{j+1}+\Delta\left(\hat{S}^{\rm
z}_{j}\hat{S}^{\rm z}_{j+1}-\frac{1}{4}\right)\right],$ (3.1)

using the generators $\hat{S}^{\alpha}=\hat{\sigma}^{\alpha}/2$ (where
$\hat{\sigma}^{\alpha}$ are the Pauli matrices) of the $\mathfrak{su}(2)$ spin
algebra, satisfying commutation relations
$[\hat{\sigma}^{\alpha},\hat{\sigma}^{\beta}]=2{\rm
i}\sum_{\gamma}\epsilon_{\alpha\beta\gamma}\hat{\sigma}^{\gamma}$. We assume
the periodic boundary condition. We shall be interested in the ferromagnetic
regime ($J>0$, $\Delta\geq 0$) and for our convenience set (with no loss of
generality) $J=1$.

There are three phases (regimes) to be distinguished: $\Delta>1$, $\Delta=1$
and $0\leq\Delta<1$, conventionally called as the gapped, isotropic and
gapless regimes, in respective order. This nomenclature refers to properties
of the antiferromagnetic ground state. Our focus in this article will be
exclusively on the gapped regime. As customary, we parametrise anisotropy as

$\Delta=\tfrac{1}{2}(q+q^{-1})=\cosh\eta,\qquad
q=\exp{\eta},\qquad\eta\in\mathbb{R},$ (3.2)

where $q$ stands for ‘deformation parameter’ of quantum symmetry algebra
$\mathcal{U}_{q}(\mathfrak{su}(2))$.

The model can be diagonalised by the Bethe Ansatz, a celebrated method
invented by Hans Bethe [37]. In what follows, we adopt the ferromagnetic
eigenstate $\mathinner{|{\uparrow}\rangle}^{\otimes L}$ as the reference
particle pseudovacuum to construct the entire spectrum of eigenstates. In
fact, the ferromagnetic states are two-fold degenerate and both ferromagnetic
vacua are required to obtain the full spectrum of eigenstates.111This is no
longer the case if the $SU(2)$ invariance is broken by a boundary twist, in
which case a single vacuum suffices. Since the $z$-component of total spin is
conserved, the Hamiltonian block-decomposes into magnetisation sectors
labelled by quantum number $M$, pertaining to the number of down-turned spins
with respect to pseudovacuum. For fixed $M$, every finite-volume eigenstate,
$|\\{\vartheta_{j}\\}^{M}_{j=1}\rangle$, is uniquely characterised by a set of
(in general complex-valued) rapidities $\vartheta_{j}$ called Bethe roots,
corresponding to solutions to a coupled system of equations

$\left[\frac{\sin(\vartheta_{j}+i\frac{\eta}{2})}{\sin(\vartheta_{j}-i\frac{\eta}{2})}\right]^{L}\prod_{k\neq
j}^{M}\frac{\sin(\vartheta_{j}-\vartheta_{k}-i\eta)}{\sin(\vartheta_{j}-\vartheta_{k}+i\eta)}=1,$
(3.3)

known as the _Bethe equations_ (BE). A suggestive physical interpretation
underneath the Bethe equations is as follows: the term in the square bracket
on the left-hand side equals $e^{{\rm i}p}$, with $p=p(\vartheta)$ being the
(bare) momentum of a single magnon excitation

$p(\vartheta)=-{\rm i}\log\frac{\sin(\vartheta+{\rm
i}\frac{\eta}{2})}{\sin(\vartheta-{\rm i}\frac{\eta}{2})},$ (3.4)

whereas the product of quotients of trigonometric functions appearing on the
right-hand side is interpreted as a net $U(1)$-valued scattering amplitude
acquired by an individual magnon upon undergoing elastic collisions with all
the remaining magnons. The total momentum and energy of an eigenstate are
obtained by summing over all the constituent magnons, yielding manifestly
additive expressions of the form

$P\big{(}\\{\vartheta_{j}\\}_{M}\big{)}=\sum_{j=1}^{M}p(\vartheta_{j}),\quad
E\big{(}\\{\vartheta_{j}\\}_{M}\big{)}=\sum_{j=1}^{M}\frac{\sin^{2}{\rm
i}\eta}{\cos 2\vartheta_{j}-\cos{\rm i}\eta}.$ (3.5)

To facilitate the asymptotic analysis of Eqs. (3.3), it is convenient to
express them in terms of the Baxter $Q$-functions [40]

$Q(\vartheta;\\{\vartheta_{j}\\})=\prod_{j=1}^{M}\sin\left(\vartheta-\vartheta_{j}\right),$
(3.6)

representing ‘trigonometric polynomials’ of degree $M$ whose zeros (in the
fundamental domain) correspond precisely to the Bethe roots of a given
eigenstate. Denoting $Q_{0}(\vartheta)\equiv\sin^{L}(\vartheta)$, and making
use of compact notations for imaginary shifts, $f^{[\pm k]}(\vartheta)\equiv
f(\vartheta\pm{\rm i}\eta/2)$, the Bethe equations can be presented in the
form

$\frac{Q_{0}^{+}(\vartheta_{j})}{Q_{0}^{-}(\vartheta_{j})}=\frac{Q_{j}^{[+2]}(\vartheta_{j})}{Q_{j}^{[-2]}(\vartheta_{j})},$
(3.7)

where $Q_{j}(\vartheta)\equiv\prod_{k\neq
j}^{M}\sin(\vartheta-\vartheta_{k})$, or in an equivalent logarithmic form

$\log Q_{0}^{+}(\vartheta_{j})-\log Q_{0}^{-}(\vartheta_{j})=2n_{j}\pi{\rm
i}+\log Q_{j}^{[+2]}(\vartheta_{j})-\log Q_{j}^{[-2]}(\vartheta_{j}).$ (3.8)

The above form is universal and provides a useful starting point to obtain the
semi-classical limits in many other quantum integrable models.

##### Thermodynamic limit.

In physics applications one is mostly interested in thermodynamic properties.
To this end, one has to infer the structure of solutions to the Bethe
equations for large system size, i.e. in the $L\to\infty$ limit with
$M\sim\mathcal{O}(L)\to\infty$ magnons. A typical thermodynamic eigenstate
comprises of elementary magnon excitations and bound states thereof. The
latter show up as compounds of complex Bethe roots (known as the Bethe
strings) whose constituent magnon rapidities with unit (imaginary) equidistant
separation (neglecting finite-size deviations that are typically exponentially
small in $L$)

$\vartheta^{(k)}_{j}=\vartheta^{(k)}+\frac{{\rm i}\eta}{2}(k+1-2j),\qquad
j=1,\ldots,k,$ (3.9)

with the centre $\vartheta^{(k)}\in\mathbb{R}$.

##### Low-energy scaling limit.

There exists another thermodynamic limit that is distinct to the one described
above. To extract the low-energy spectrum of spin fluctuations at long
wavelengths, one has to consider a thermodynamic _scaling_ limit by taking
both $L$ and $M\sim\mathcal{O}(L)$ large, while additionally demanding that
all magnons have low momenta scaling as $\mathcal{O}(1/L)$. This way we are
left with a finite number $m\sim\mathcal{O}(1)$ macroscopic bound states. We
stress that, in this regime, the Bethe’s original argument for the string
formation is no longer valid. What happens instead is that bound states become
less dense and in general appreciably deformed.

This low-energy eigenstates have been first investigated by Sutherland [41]
and afterwards also in Ref. [42], where they are dubbed as ‘quantum Bloch
walls’. Their classical nature has however been elucidated only later in Ref.
[12], establishing an explicit connection with (nonlinear) spin waves governed
by the continuous Landau–Lifshitz ferromagnet with isotropic interaction.

The method outlined in Ref. [12] is rather general and can be extended to
accommodate also for the interaction anisotropy. As shown below, in the
presence of the interaction anisotropy $\eta$, the thermodynamic scaling limit
that governs the semi-classical spectrum of low-energy eigenstates requires to
additionally assume that anisotropy is weak, $\Delta\gtrapprox 1$. To be more
specific, after reinstating lattice spacing $a$ and writing $\ell\equiv
L\,a\in\mathcal{O}(1)$, anisotropy parameter $\eta\sim\mathcal{O}(1/L)\to 0$
has to be rescaled as

$\eta=\frac{\epsilon\,\ell}{L}+\mathcal{O}\left(\frac{1}{L^{2}}\right),$
(3.10)

with parameter

$\epsilon\equiv\sqrt{\delta},\qquad\delta=\frac{2(\Delta-1)}{a^{2}}\in\mathcal{O}(1),$
(3.11)

kept fixed whilst taking the continuum thermodynamic limit, $L\to\infty$ and
$a\to 0$. Here parameter $\ell$ plays the role of a length and later on, in
Sec. 3.1, we show that it corresponds precisely to the circumference of the
emergent classical phase space.

By first expanding the logarithm of $Q_{0}^{\pm}$ we find

$\log Q_{0}^{\pm}(\vartheta_{j})=\log
Q_{0}(\vartheta_{j})\pm\frac{i\eta}{2}\frac{{\rm d}}{{\rm d}\vartheta}\log
Q_{0}(\vartheta)\rvert_{\vartheta=\vartheta_{j}}-\frac{\eta^{2}}{8}\frac{{\rm
d}^{2}}{{\rm d}\vartheta^{2}}\log
Q_{0}(\vartheta)\rvert_{\vartheta=\vartheta_{j}}+\mathcal{O}\left(\frac{1}{L^{2}}\right),$
(3.12)

where we have assumed that $\vartheta\in\mathcal{O}(1)$. At this stage it is
convenient to perform a change of variable by introducing

$\mu=\frac{\eta}{\ell}\frac{{\rm d}}{{\rm d}\vartheta}\log
Q_{0}(\vartheta)=\frac{\epsilon}{\tan\vartheta},\qquad\mu_{j}=\lim_{\vartheta\to\vartheta_{j}}\mu,$
(3.13)

which then yields

$\log Q_{0}^{+}(\vartheta_{j})-\log Q_{0}^{-}(\vartheta_{j})={\rm
i}\ell\mu_{j}+\mathcal{O}\left(\frac{1}{L^{2}}\right).$ (3.14)

The part involving $Q(\vartheta)$ can be treated analogously. The details are
given in Appendix B. By combining the two contributions, we finally arrive at
the following compact representation in the $\mu$-variable 222Exact solutions
to this algebraic equation for a single mode number can be found in Ref. [43].
A similar construction for the isotropic case appeared in Ref. [44].

$\mu_{j}=\frac{2\pi}{\ell}n_{j}-\frac{2}{L}\sum_{k\neq
j}^{M}\frac{\mu_{j}\mu_{k}+\delta}{\mu_{j}-\mu_{k}}+\mathcal{O}\left(\frac{1}{L}\right).$
(3.15)

Upon taking the thermodynamic scaling limit, the Bethe roots $\mu_{j}$
condense along certain one-dimensional segments (contours). In general, there
are several disjoint contours $\mathcal{C}\equiv\cup_{j}\mathcal{C}_{j}$.
Accordingly, the Bethe equations turn into singular integral equations of the
form

$\frac{\ell\mu}{2}=\pi
n_{j}-\ell\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{\mathcal{C}}{\rm
d}\lambda\,\mathcal{K}_{\delta}(\mu,\lambda)\rho(\lambda),\qquad\mu\in\mathcal{C}_{j},$
(3.16)

with integral kernel

$\mathcal{K}_{\delta}(\mu,\lambda)\equiv\frac{\mu\lambda+\delta}{\mu-\lambda}.$
(3.17)

Eqs. (3.16) are known as the _asymptotic Bethe equations_ (ABE). To satisfy
the reality constraint, the Bethe roots must appear in complex-conjugate
pairs, implying that contours $\mathcal{C}_{j}$ are symmetric under reflection
about the real axis. The leading correction term to Eq. (3.15) is of the order
$\mathcal{O}(1/L)$ (cf. Appendix B),

$\frac{1}{L}\pi\rho^{\prime}(\mu)\ell^{2}(\mu^{2}+\delta)^{2}\coth\left[\pi\ell(\mu^{2}+\delta)\rho(\mu)\right].$
(3.18)

It turns out that this correction to Eqs. (3.16) can only be safely discarded
when

$\pi\ell(\mu^{2}+\delta)\rho(\mu)\neq{\rm i}\pi n,\qquad n\in\mathbb{Z},$
(3.19)

which is satisfied when the density of roots near the real axis is
sufficiently low. In contrast, when the density of Bethe roots is high enough,
the assumptions underlying the above perturbative expansion are no longer
justified and one has to resort to the full quantum Bethe equations non-
perturbatively (at least in the region near the real axis where the effect is
most pronounced). We shall return to this subtlety in Sec. 6.2 and examine it
closely on a specific class of solutions. As it turns out, the effect is
responsible for emergence of certain special features in the solutions called
_condensates_ [13].

##### Riemann–Hilbert problem.

The asymptotic Bethe equations (3.16) can be formulated as a Riemann-Hilbert
problem. To this end, we define the _spectral resolvent_

$G(\mu)=\ell\int_{\mathcal{C}}{\rm
d}\lambda\mathcal{K}_{\delta}(\mu,\lambda)\rho(\lambda),$ (3.20)

and define

$\mathfrak{p}(\mu)=G(\mu)+\frac{\ell\mu}{2}.$ (3.21)

At every point $\mu$ along the density contour $\mathcal{C}_{j}$ function
$\mathfrak{p}(\mu)$ experiences a jump discontinuity that is proportional to
the density (of Bethe roots) $\rho(\mu)$,

$\mathfrak{p}(\mu+{\rm i}0)-\mathfrak{p}(\mu-{\rm i}0)=2{\rm
i}\pi\ell(\mu^{2}+\delta)\rho(\mu),\qquad\mu\in\mathcal{C}_{j},$ (3.22)

with $\pm{\rm i}0$ denoting infinitesimal displacements to either side of the
contour. 333Orientation of integration contours is in the direction towards
the branch point with a negative imaginary part. Individual contours
$\mathcal{C}_{j}$ can thus be pictured as the $j$th branch cut (of square-root
type) of a two-sheeted Riemann surface. The end points of $\mathcal{C}_{j}$
thus correspond to branch points. In this view, function $\mathfrak{p}(\mu)$
is a double-valued complex function which, apart from contours
$\mathcal{C}_{j}$, is analytic everywhere on the complex $\mu$-plane. A branch
cut of square-root type implies that upon crossing it we end up on another
Riemann sheet and $\mathfrak{p}(\mu)$ flips its sign. In addition,
$\mathfrak{p}(\mu)$ in general picks up an integer multiple of $2\pi$, namely

$\mathfrak{p}(\mu+{\rm i}0)+\mathfrak{p}(\mu-{\rm i}0)=2\pi
n_{j},\qquad\mu\in\mathcal{C}_{j}.$ (3.23)

Remarkably, $\mathfrak{p}(\mu)$ is precisely the classical quasi-momentum
pertaining to the completely integrable classical anisotropic ferromagnet. As
explained and demonstrated in the next section, quasi-momentum encodes the
eigenvalues of the classical monodromy matrix obtained from a path-ordered
exponential of the classical Lax connection.

In the Sec. 6 we shall make a direct comparison between the branch cuts and
discrete distributions of Bethe roots for finite system sizes. There we find
it useful to rewrite the Riemann-Hilbert problem in terms of variable
$\zeta=1/\mu$. In Appendix A we provide a dictionary between the two
parametrisations used.

###  3.1 Anisotropic Landau–Lifshitz field theory

The task at hand is to derive the classical equation of motion for the low-
energy spectrum of the Heisenberg XXZ chain. We present how to achieve this in
a systematic fashion. The first step is to infer the spatial component of the
classical Lax connection from the semi-classical expansion of the quantum
monodromy matrix. This will also enable us to explicitly establish that
$\mathfrak{p}(\mu)$ from the previous section is truly the classical quasi-
momentum that belongs to the classical _axially anisotropic Landau–Lifshitz
model_. Written in terms of a $S^{2}$-valued classical spin field
$\vec{\mathcal{S}}(x,t)$, the latter is described by the following partial
differential equation

$\vec{\mathcal{S}}_{t}=\vec{\mathcal{S}}\times\vec{\mathcal{S}}_{xx}+\vec{\mathcal{S}}\times\mathrm{J}\vec{\mathcal{S}},\qquad\mathrm{J}\,\equiv\mathrm{diag}(0,0,\delta).$
(3.24)

Here and subsequently we make use of compact notations
$\vec{\mathcal{S}}_{t}\equiv\partial_{t}\vec{\mathcal{S}}(x,t)$,
$\vec{\mathcal{S}}_{x}\equiv\partial_{x}\vec{\mathcal{S}}(x,t)$, and similarly
for higher partial derivatives. The Hamiltonian that generates Eq. (3.24) is
of the form

$\mathcal{H}=\frac{1}{2}\int_{0}^{\ell}{\rm
d}x\left[\vec{\mathcal{S}}_{x}(x)\cdot\vec{\mathcal{S}}_{x}(x)+\delta-\delta(\mathcal{S}^{\rm
z}(x))^{2}\right].$ (3.25)

##### Zero-curvature representation.

Complete integrability of Eq. (3.24) can be made manifest by recasting it in
the form of a _zero-curvature condition_

${\bf U}_{t}(\mu;x,t)-{\bf V}_{x}(\mu;x,t)+[{\bf U}(\mu;x,t),{\bf
V}(\mu;x,t)]=0.$ (3.26)

The latter ensures compatibility of the following _auxiliary linear problem_ ,

$\Psi_{x}(\mu;x,t)={\bf U}(\mu;x,t)\Psi(\mu;x,t),\qquad\Psi_{t}(\mu;x,t)={\bf
V}(\mu;x,t)\Psi(\mu;x,t),$ (3.27)

where [14]

${\bf U}(\mu;x,t)=\frac{1}{2{\rm i}}\begin{pmatrix}\mu\,\mathcal{S}^{\rm
z}&\sqrt{\mu^{2}+\delta}\,\mathcal{S}^{-}\\\
\sqrt{\mu^{2}+\delta}\,\mathcal{S}^{+}&-\mu\,\mathcal{S}^{\rm
z}\end{pmatrix}.$ (3.28) ${\bf V}(\mu;x,t)=\frac{{\rm
i}}{2}\begin{pmatrix}(\mu^{2}+\delta)\mathcal{S}^{\rm
z}&\mu\sqrt{\mu^{2}+\delta}\mathcal{S}^{-}\\\
\mu\sqrt{\mu^{2}+\delta}\mathcal{S}^{+}&-(\mu^{2}+\delta)\mathcal{S}^{\rm
z}\end{pmatrix}\\\ -\frac{1}{2{\rm i}}\begin{pmatrix}\mu\,\mathcal{J}_{0}^{\rm
z}&\sqrt{\mu^{2}+\delta}\,\mathcal{J}_{0}^{-}\\\
\sqrt{\mu^{2}+\delta}\,\mathcal{J}_{0}^{+}&-\mu\,\mathcal{J}_{0}^{\rm
z}\end{pmatrix},$ (3.29)

are the spatial and the temporal component of the classical Lax connection,
respectively, with

$\mathcal{J}_{0}\equiv\vec{\mathcal{S}}_{x}\times\vec{\mathcal{S}},$ (3.30)

denoting the spin current density at $\delta=0$.

##### Semi-classical limit.

We next show how the above (classical) Lax connection can be retrieved from
the Lax operator of the quantum chain. The quantum Lax operator is the
elementary building block of commuting transfer matrices which facilitate
algebraic diagonalisation of the XXZ Hamiltonian (3.1). The fundamental Lax
operator $\mathbf{L}(\vartheta)$ acts on a one-site (physical) Hilbert space
$\mathcal{V}_{p}\cong\mathbb{C}^{2}$ of a spin-$1/2$ degree of freedom and an
auxiliary space $\mathcal{V}_{a}$ associated to the fundamental representation
of the quantum group $\mathcal{U}_{q}(\mathfrak{su}(2))$ (with deformation
parameter $q=e^{\eta})$, and depends analytically on the complex (spectral)
parameter $\vartheta$. It reads explicitly

${\bf L}(\vartheta)=\frac{1}{\sinh\eta}\begin{pmatrix}\sin(\vartheta+{\rm
i}\eta\mathbf{S}^{\rm z})&i\sinh\eta\mathbf{S}^{-}\\\
i\sinh\eta\mathbf{S}^{+}&\sin(\vartheta-{\rm i}\eta\mathbf{S}^{\rm
z})\end{pmatrix},$ (3.31)

in terms of auxiliary spin generators which enclose the $q$-deformed
commutation relations

$[\mathbf{S}^{+},\mathbf{S}^{-}]=[2\mathbf{S}^{z}]_{q},\qquad
q^{2\mathbf{S}^{\rm z}}\mathbf{S}^{\pm}=q^{\pm
2}\mathbf{S}^{\pm}q^{2\mathbf{S}^{\rm z}},$ (3.32)

using the standard notation $[x]_{q}=\sinh(\eta\,x)/\sinh\eta$.

We proceed by constructing the fundamental row transfer matrix of the quantum
XXZ spin chain, obtained as a partial trace (over the fundamental auxiliary
space $\mathcal{V}_{a}\cong\mathbb{C}^{2}$) of the monodromy matrix, i.e. an
ordered product of Lax matrices

$T(\vartheta)=\mathrm{Tr}_{\mathcal{V}_{a}}\,\mathbf{M}(\vartheta)=\mathrm{Tr}_{\mathcal{V}_{a}}{\bf
L}^{(1)}(\vartheta){\bf L}^{(2)}(\vartheta)\cdots{\bf L}^{(N)}(\vartheta).$
(3.33)

Here we have adopted the right-to-left ordering convention and used ${\bf
L}^{(k)}(\vartheta)$ to denote the embedding into the $k$th factor in an
$L$-fold product Hilbert space $\mathcal{H}\cong\mathcal{V}^{\otimes L}_{p}$.
By virtue of the quantum Yang–Baxter relation, matrices $T(\vartheta)$
mutually commute, $[T(\vartheta),T(\vartheta^{\prime})]=0$ for all
$\vartheta,\vartheta^{\prime}\in\mathbb{C}$. Therefore, commuting transfer
matrices serve as the generating operator for the local (and nonlocal)
conserved charges [38, 12]. An infinite tower of commuting fused transfer
matrices $T_{j}(\vartheta)$ with $(j+1)$-dimensional auxiliary unitary
representations of $\mathcal{U}_{q}(\mathfrak{su}(2))$ can be constructed in a
similar manner, providing additional quasilocal conservation laws of the XXZ
model. While these are of utmost importance for thermodynamic properties at
finite energy density (see e.g. Refs. [45, 46, 47, 48]), they do not play any
role upon taking the semi-classical limit.

We have thus far analysed the asymptotic scaling limit $L\to\infty$ at the
level of the Bethe equations by parametrising the interaction parameter as
$\eta=\epsilon\,\ell/L\to 0$ and subsequently taking $\eta\to 0$ ($q\to 1$).
Now we do the same at the level of the transfer matrix, where we are allowed
to substitute the $q$-deformed spin generators with the fundamental
$\mathfrak{su}(2)$ spins,
$\mathbf{S}^{\alpha}\to\hat{S}^{\alpha}=\tfrac{1}{2}\hat{\sigma}^{\alpha}$ for
$\alpha\in\\{{\rm x,y,z}\\}$. In this limit, the diagonal elements of the
quantum Lax operator become

$\lim_{\eta\to 0}\frac{\sin(\vartheta\mathds{1}+{\rm i}\eta\mathbf{S}^{\rm
z})}{\sinh\eta}=\frac{\sin\vartheta}{\eta}\mathds{1}+{\rm
i}\cos{(\vartheta)}\hat{S}^{\rm z}+\mathcal{O}(\eta),$ (3.34)

whence

${\bf L}(\vartheta)\simeq\frac{\sin\vartheta}{\eta}\left[\mathds{1}+{\rm
i}\eta\begin{pmatrix}\cot{(\vartheta)}\hat{S}^{\rm
z}&\csc{(\vartheta)}\hat{S}^{-}\\\
\csc{(\vartheta)}\hat{S}^{+}&-\cot{(\vartheta)}\hat{S}^{\rm
z}\end{pmatrix}\right].$ (3.35)

By reinstating the lattice spacing $a=\ell/L$ and using the spectral parameter

$\mu=\frac{\epsilon}{\tan{\vartheta}},$ (3.36)

the Lax matrix writes

${\bf L}(\mu)\simeq\frac{1}{a\sqrt{\mu^{2}+\delta}}\left[\mathds{1}+{\rm
i}a\begin{pmatrix}\mu\hat{S}^{\rm z}&\sqrt{\mu^{2}+\delta}\hat{S}^{-}\\\
\sqrt{\mu^{2}+\delta}\hat{S}^{+}&-\mu\hat{S}^{\rm z}\end{pmatrix}\right],$
(3.37)

whereas the asymptotic scaling limit of the associated monodromy matrix ${\bf
M}(\mu)$ is given by the following path-ordered product

$\begin{split}\frac{{\bf
M}(\mu)}{(a\sqrt{\mu^{2}+\delta})^{L}}&\sim\prod_{j=L}^{1}\left[\mathds{1}+{\rm
i}a\begin{pmatrix}\mu\hat{S}^{\rm
z}_{j}&\sqrt{\mu^{2}+\delta}\,\hat{S}^{-}_{j}\\\
\sqrt{\mu^{2}+\delta}\,\hat{S}^{+}_{j}&-\mu\,\hat{S}^{\rm
z}_{j}\end{pmatrix}\right].\end{split}$ (3.38)

In the final step we replace the quantum spins $\hat{S}^{\alpha}$ with
classical spin variables via
$\mathcal{S}^{\alpha}_{j}\equiv\mathcal{S}^{\alpha}(x=j\,a)$, and subsequently
take the continuum limit. We thus arrive as the following semi-classical
approximation of the quantum monodromy matrix

${\bf M}_{\rm cl}(\mu)\equiv\mathscr{P}\exp\left[\frac{{\rm
i}}{2}\oint_{0}^{\ell}\frac{{\rm d}x}{2\pi}{\bf U}(\mu;x,t)\right],$ (3.39)

with ${\bf U}(\mu;x,t)$ being the spatial component of the Lax connection
introduced earlier in Eq. (3.28).

##  4 Finite-gap integration method

Having retrieved the classical quasi-momentum $\mathfrak{p}(\mu)$ from the
semi-classical expansion of the Heisenberg quantum chain, we proceed to
explain how its analytic structure encodes the spectrum of interacting
nonlinear phases that characterise classical spin-field configurations. We
shall confined our considerations, as usual, to a class of solutions that
involve only a finite number of excited nonlinear modes (corresponding to
Riemann surfaces of finite genus), commonly known in the literature as the
finite-gap solutions [49, 50, 29].

We next describe the full programme for performing an algebro-geometric
integration of completely integrable nonlinear partial differential equations
which permits to solve the initial value problem. The main steps comprise of:

  1. 1.

prescribing an appropriately normalised meromorphic differential ${\rm
d}\mathfrak{p}$ on a hyperelliptic Riemann surface of finite genus,

  2. 2.

constructing the fundamental matrix solution to the associated auxiliary
linear problem,

  3. 3.

identifying the dynamical separated variables and deriving their equations of
motion,

  4. 4.

employing the Abel–Jacobi transformation to obtain canonical action-angle
variables satisfying a linear evolution law on a Jacobian hypertorus,

  5. 5.

expressing physical fields through the solution to the inverse problem.

We stress that the outlined finite-gap integration scheme is not particular to
the model at hand but rather applies universality. Moreover, the method has
been previously developed also for the Landau–Lifshitz model in Refs. [51,
52]. Nevertheless, we wish to offer a different formulation here that we find
conceptually simpler. Specifically, we shall avoid the conventional use of the
Baker–Akhiezer vectors. Our plan is to first discuss some general aspects and
then to provide two explicit realisations for $\mathfrak{g}=0$ and
$\mathfrak{g}=1$. We do not pay much attention to the construction of the
action-angle variables but rather give an explicit time dependence of the
physical spin fields in terms of the Riemann theta functions.

###  4.1 Auxiliary linear problem

In this section we describe how to solve the auxiliary linear problem (3.27).
By formally integrating along the spatial direction at a fixed time-slice, we
have

$\boldsymbol{\psi}(\mu;x)=\mathscr{P}\,\exp{\left(\int_{x_{0}}^{x}{\rm
d}x^{\prime}\,{\bf U}(\mu;x^{\prime})\right)}\equiv{\bf T}_{\rm cl}(x,x_{0}).$
(4.1)

By imposing periodic boundary conditions
$\vec{\mathcal{S}}(x+\ell)=\vec{\mathcal{S}}(x)$, we define the monodromy
matrix (3.39)

${\bf M}_{\rm cl}(\mu)={\bf T}_{\rm cl}(x_{0}+\ell,x_{0}).$ (4.2)

According to the Bloch theorem, there are two linearly-independent solutions
that linearise this translation,

$\boldsymbol{\psi}(\mu;x+\ell)=\boldsymbol{\Lambda}(\mu)\,\boldsymbol{\psi}(\mu;x),$
(4.3)

where Bloch multiplier $\boldsymbol{\Lambda}(\mu)={\rm
diag}(\Lambda_{+}(\mu),\Lambda_{-}(\mu))$ is given by eigenvalues of ${\bf
M}_{\rm cl}(\mu)$, which (by virtue of the zero-curvature condition (3.26))
are conserved under time evolution, $\partial\Lambda_{\pm}(\mu)/\partial t=0$.
Since ${\rm Tr}\,{\bf U}_{\rm cl}(\mu)=0$, monodromy matrix ${\bf M}_{\rm
cl}(\mu)$ is unimodular and hence its eigenvalues can be parametrised as
$\Lambda_{\pm}(\mu)=\exp{(\pm{\rm i}\mathfrak{p}(\mu))}$ in terms of a single
complex-valued function $\mathfrak{p}(\mu)$ called quasi-momentum.
Subsequently we make use of the following general parametrisation

${\bf M}_{\rm cl}(\mu)=\cos{(\mathfrak{p}(\mu))}+{\rm
i}\sin{(\mathfrak{p}(\mu))}\boldsymbol{\Psi}(\mu;x_{0}).$ (4.4)

In analogy with the harmonic oscillator presented in the introductory section
2, we introduce the _squared eigenfunctions_ through

$\boldsymbol{\Psi}(\mu;x_{0})=\boldsymbol{\psi}(\mu;x_{0})\,\boldsymbol{\sigma}^{\rm
z}\,\boldsymbol{\psi}(\mu;x_{0})^{-1},$ (4.5)

representing a periodic matrix solution (that is
$\boldsymbol{\Psi}(\mu;x+\ell)=\boldsymbol{\Psi}(\mu;x)$) to the _adjoint_
linear system

$\boldsymbol{\Psi}_{x}(\mu;x,t)=\big{[}{\bf
U}(\mu;x,t),\boldsymbol{\Psi}(\mu;x,t)\big{]},\qquad\boldsymbol{\Psi}_{t}(\mu;x,t)=\big{[}{\bf
V}(\mu;x,t),\boldsymbol{\Psi}(\mu;x,t)\big{]}.$ (4.6)

##### Uniformisation.

Introducing an uniformised spectral parameter $z$,

$\mu=\frac{1}{2}\left(z-\frac{\delta}{z}\right),\qquad\sqrt{\mu^{2}+\delta}=\frac{1}{2}\left(z+\frac{\delta}{z}\right),$
(4.7)

the solution to the above linear problem can be sought in the form of a formal
Laurent series

$\boldsymbol{\Psi}(\mu)=\sum_{n=0}^{\infty}\frac{\boldsymbol{\Psi}_{n}}{z^{n}}.$
(4.8)

Defining $S^{2}$-valued matrices

${\bf
S}\equiv\vec{\boldsymbol{\sigma}}\otimes\vec{\mathcal{S}}=\begin{pmatrix}\mathcal{S}^{\rm
z}&\mathcal{S}^{-}\\\ \mathcal{S}^{+}&-\mathcal{S}^{\rm z}\\\
\end{pmatrix},\qquad\widetilde{\bf S}\equiv\boldsymbol{\sigma}^{\rm z}\,{\bf
S}\,\boldsymbol{\sigma}^{\rm z},$ (4.9)

the first few terms can be written compactly in the form

$\boldsymbol{\Psi}_{0}={\bf S},\quad\boldsymbol{\Psi}_{1}={\rm i}[{\bf S},{\bf
S}_{x}],\quad\boldsymbol{\Psi}_{2}=-{\rm Tr}({\bf S}_{x})^{2}{\bf S}-[{\bf
S},[{\bf S},{\bf S}_{xx}]]+\frac{\epsilon^{2}}{4}[{\bf S},[\widetilde{\bf
S},{\bf S}]],$ (4.10)

and so forth.

##### Local conserved charges.

The quasi-momentum can be related to the squared eigenfunctions via [14]

$\mathfrak{p}(z)=-{\rm i}\int^{\ell}_{0}{\rm d}x\,{\rm Tr}\big{[}{\bf
U}(\boldsymbol{\Psi}+\boldsymbol{\sigma}^{\rm
z})^{-1}\big{]}=-\frac{z\,\ell}{4}+\sum_{n=0}^{\infty}\frac{Q_{n}}{z^{n}}.$
(4.11)

where coefficients $Q_{n}$ provide (extensive) local conserved charges for a
given solution. The first two take the form (modulo total derivatives)

$\begin{split}Q_{0}&=\frac{{\rm i}}{4}\int^{\ell}_{0}{\rm
d}x\,\frac{\mathcal{S}^{-}\mathcal{S}^{+}_{x}-\mathcal{S}^{+}\mathcal{S}^{-}_{x}}{1+\mathcal{S}^{\rm
z}}=-\frac{\mathcal{P}}{2},\\\ Q_{1}&=-\frac{1}{2}\int^{\ell}_{0}{\rm
d}x\,\Big{[}{\vec{\mathcal{S}}}^{2}_{x}+\epsilon^{2}\big{(}1-(\mathcal{S}^{\rm
z})^{2}\big{)}-\frac{\epsilon^{2}}{2}\Big{]}=-\mathcal{H}+\frac{\epsilon^{2}\ell}{4},\end{split}$
(4.12)

which yields

$\mathfrak{p}(\mu)=-\frac{\mu\ell}{2}-\frac{\mathcal{P}}{2}-\frac{\mathcal{H}}{\mu}+\mathcal{O}(\mu^{-2}).$
(4.13)

###  4.2 Finite-gap solutions

The quasi-momentum $\mathfrak{p}(\mu)$ contains only information about the
conserved quantities encoded in a particular spectral curve. In order to
reconstruct the spin field $\vec{\mathcal{S}}(x,t)$ from the spectral data, we
also need knowledge of the dynamical degrees of freedom. Below we outline the
algebro-geometric procedure to infer the time evolution of the spin field
$\vec{\mathcal{S}}(x,t)$ for the class of finite-gap solutions.

To achieve this, we parametrise the squared eigenfunction as

$\boldsymbol{\Psi}_{\mathfrak{g}}(\mu)=\frac{1}{\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\mu)}}\begin{pmatrix}a_{\mathfrak{g}+1}(\mu)&\sqrt{\mu^{2}+\epsilon^{2}}\,b_{\mathfrak{g}}(\mu)\\\
\sqrt{\mu^{2}+\epsilon^{2}}\,\bar{b}_{\mathfrak{g}}(\mu)&-a_{\mathfrak{g}+1}(\mu)\end{pmatrix},$
(4.14)

where functions $a_{d}(\mu)$, $b_{d}(\mu)$ and $\mathcal{R}_{d}(\mu)$ are
_polynomials_ in variable $\mu$ of degree $d$ In particular,
$\mathcal{R}_{2\mathfrak{g}+2}(\mu)$ is a polynomial of degree
$2\mathfrak{g}+2$,

$\mathcal{R}_{2\mathfrak{g}+2}(\mu)=\prod_{j=1}^{\mathfrak{g}+1}(\mu-\mu_{j})(\mu-\bar{\mu}_{j})=\sum_{k=0}^{2\mathfrak{g}+2}(-1)^{k}r_{2\mathfrak{g}+2-k}\mu^{k},$
(4.15)

which specifies a hyperelliptic algebraic curve in $\mathbb{C}^{2}$,

$\Sigma:\qquad y^{2}(\mu)=\mathcal{R}_{2\mathfrak{g}+2}(\mu).$ (4.16)

The curve is fully characterised by $2(\mathfrak{g}+1)$ branch points
$\mu_{j}$ (or, equivalently, symmetric polynomials $r_{k}$ thereof, with
$r_{2\mathfrak{g}+2}=1$). By the unimodularity constraint,
$\det\boldsymbol{\Psi}_{\mathfrak{g}}=-1$, functions $a(\mu)$ and $b(\mu)$ are
not independent but are instead subjected to an algebraic relation

$a^{2}_{\mathfrak{g}+1}(\mu)+\left(\mu^{2}+\delta\right)b_{\mathfrak{g}}(\mu)\bar{b}_{\mathfrak{g}}(\mu)=\mathcal{R}_{2\mathfrak{g}+2}(\mu).$
(4.17)

From the trace identity (4.11) one can readily infer the following compact
expression for the quasi-momentum

$\mathfrak{p}(\mu)=-\frac{\mu}{2}\int^{\ell}_{0}{\rm d}x\,\mathcal{S}^{\rm
z}(x)-\frac{\mu^{2}+\delta}{4}\int^{\ell}_{0}{\rm
d}x\,\frac{\mathcal{S}^{-}\bar{b}_{\mathfrak{g}}(\mu)+\mathcal{S}^{+}b_{\mathfrak{g}}(\mu)}{\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\mu)}+a_{\mathfrak{g}+1}(\mu)}.$
(4.18)

This form is compatible with the correct asymptotic expansion about
$\mu\to\infty$. A series expansion of $\mathfrak{p}(\mu)$ will thus involve
only $(\mathfrak{g}+1)$ functionally independent integrals of motion $Q_{n}$.
They can be expressed as certain functions of coefficients $r_{j}$ of
$\mathcal{R}_{2\mathfrak{g}+2}(\mu)$. Moreover, the total filling fraction
$\nu$ with respect to ferromagnetic vacuum $\mathcal{S}^{\rm z}_{\rm vac}=1$,

$\nu\equiv\frac{1}{2\ell}\int^{\ell}_{0}{\rm d}x\,(1-\mathcal{S}^{\rm z}(x)),$
(4.19)

can be obtained as

$\mathfrak{p}(\mu=\pm{\rm i}\epsilon)=\mp\frac{{\rm
i}\epsilon}{2}\int^{\ell}_{0}{\rm d}x\,\mathcal{S}^{\rm z}(x)=\mp\frac{{\rm
i}\epsilon\ell}{2}(1-2\nu).$ (4.20)

####  4.2.1 Dynamical divisor

Off-diagonal elements of the fundamental matrix solution
$\boldsymbol{\Psi}_{\mathfrak{g}}(\mu)$, i.e. dynamical zeros of the
$b$-function, provide _dynamical_ degrees of freedom of finite-gap solutions
444 Indeed, function $b_{\mathfrak{g}}(\mu)$ can be interpreted as the
classical analogue of the ${\bf B}$-operator, namely the off-diagonal element
of the quantum monodromy matrix ${\bf B}\equiv{\bf M}_{12}$, whose operator-
valued zeros are the quantum separated variables [53]., enabling the
reconstruction of the time-evolved spin field $\vec{\mathcal{S}}(x,t)$.

To satisfy quadratic constraint (4.17), we parametrise

$b_{\mathfrak{g}}(\mu;x,t)=\mathcal{S}^{-}(x,t)\prod_{j=1}^{\mathfrak{g}}\big{(}\mu-\gamma_{j}(x,t)\big{)}.$
(4.21)

where the leading coefficient has been fixed to match the asymptotics at
$\mu\to\infty$.

The set $\mathcal{D}=\\{\gamma_{j}\\}_{j=1}^{g}$ is known as the _dynamical
divisor_ of the Riemann surface $\Sigma$. Since $\mathcal{S}^{-}(x,t)$ should
also be regarded as an independent dynamical variable, we have in total
$(\mathfrak{g}+1)$ dynamical degrees of freedom. This number exactly matches
the number of action variables and corresponds to the number of forbidden
zones in the finite-gap spectrum.

##### Dubrovin equations.

Dynamics of variables $\gamma_{j}=\gamma_{j}(x,t)$ takes place on the Riemann
surface $\Sigma$. Below we employ an extended dynamical divisor
$\mathcal{D}_{\rm ext}$ by adjoining it two extra non-dynamical variables
$\gamma_{\pm}\equiv\pm{\rm i}\epsilon$ which we label by
$\gamma_{\mathfrak{g}+1}$ and $\gamma_{\mathfrak{g}+2}$, respectively. Using
Lagrange interpolation formula, we can then restore functions
$a_{\mathfrak{g}+1}(\mu)$ from Eq. (4.17), yielding

$a_{\mathfrak{g}+1}(\mu;x,t)=\sum_{j=1}^{\mathfrak{g}+2}\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\gamma_{j}(x,t))}\prod_{k\neq
j}^{\mathfrak{g}+2}\frac{\mu-\gamma_{k}(x,t)}{\gamma_{j}(x,t)-\gamma_{k}(x,t)}.$
(4.22)

The equations of motion (4.6) in terms of $\gamma$-variables from
$\mathcal{D}_{\rm ext}$ take the form

$\begin{split}\partial_{x}\gamma_{j}(x,t)&={\rm
i}\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\gamma_{j}(x,t))}\prod_{k\neq
j}^{\mathfrak{g}}(\gamma_{j}(x,t)-\gamma_{k}(x,t))^{-1},\\\
\partial_{t}\gamma_{j}(x,t)&={\rm
i}\,\Gamma(x,t)\,\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\gamma_{j}(x,t))}\prod_{k\neq
j}^{\mathfrak{g}}(\gamma_{j}(x,t)-\gamma_{k}(x,t))^{-1},\end{split}$ (4.23)

with

$\Gamma(x,t)\equiv\frac{r_{1}}{2}-\sum_{k\neq
j}^{\mathfrak{g}}\gamma_{k}(x,t),\qquad
r_{1}=\sum_{j=1}^{\mathfrak{g}+1}(\mu_{j}+\bar{\mu}_{j}).$ (4.24)

We have obtained a system of differential equations that governs the motion of
the dynamical divisor of a Riemann surface, commonly known in the literature
under the name of _Dubrovin equations_ [49, 54, 29]. The form of these
equations is universal, namely they do not depend on the model under
consideration. What is model-specific instead are the reconstruction formulae,
i.e. how $\gamma$-variables relate to physical fields. A spin field whose
target space is a $2$-sphere can be described by two degrees of freedom, e.g.
the $\mathcal{S}^{\rm z}$ and $\mathcal{S}^{-}$ components. These can be
restored from Eqs. (4.6), which yields

$\mathcal{S}^{\rm
z}(x,t)=\sum_{j=1}^{\mathfrak{g}+2}\frac{\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\gamma_{j}(x,t))}}{\prod_{k\neq
j}^{\mathfrak{g}+2}(\gamma_{j}(x,t)-\gamma_{k}(x,t))},$ (4.25)

and

${\rm
i}\frac{\mathcal{S}^{-}_{x}(x,t)}{\mathcal{S}^{-}(x,t)}=\sum_{j=1}^{\mathfrak{g}+2}\frac{\gamma_{j}\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\gamma_{j}(x,t))}}{\prod_{k\neq
j}^{\mathfrak{g}+2}(\gamma_{j}(x,t)-\gamma_{k}(x,t))},$ (4.26)

${\rm i}\frac{\mathcal{S}^{-}_{t}(x,t)}{\mathcal{S}^{-}(x,t)}={\rm
i}\frac{r_{1}}{2}\frac{\mathcal{S}^{-}_{x}(x,t)}{\mathcal{S}^{-}(x,t)}-\sum_{j=1}^{\mathfrak{g}+2}\frac{\gamma_{j}\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\gamma_{j})}\left(\sum_{k\neq
j}^{\mathfrak{g}+2}\gamma_{k}(x,t)\right)}{\prod_{k\neq
j}^{\mathfrak{g}+2}(\gamma_{j}(x,t)-\gamma_{k}(x,t))}.$ (4.27)

##### Abel–Jacobi transformation.

Dubrovin equations (4.23) allow for exact integration. To this end, we define
the standard basis of $2\mathfrak{g}$ closed cycles on $\Sigma$. They further
split into $\mathcal{A}$-cycles and their conjugate $\mathcal{B}$-cycles
according to the following prescription: $\mathcal{A}_{j}$-cycle encircles the
the $j$th branch cut $\mathcal{C}_{j}$ on the upper Riemann sheet of $\Sigma$,
whereas $\mathcal{B}_{jk}$ denotes a cycle that passes through cut
$\mathcal{C}_{j}$ on the upper sheet and closes back to itself through
$\mathcal{C}_{k}$, as shown in Fig. 2.

Dubrovin equations for the dynamical divisor on $\Sigma$ can be integrated
with aid of the Abel–Jacobi transformation,

$\varphi_{j}=2\pi\sum_{k=1}^{\mathfrak{g}}\int^{\gamma_{k}(x,t)}_{\gamma_{k}(0,0)}\omega_{j},$
(4.28)

where $\omega_{j}$ form the basis of holomorphic differentials of the Riemann
surface. The above mapping provides a variable transformation from
$\gamma$-variables to $\varphi$-variables of the angle-type,
$\\{\gamma_{j}(x,t)\\}\mapsto\\{\varphi_{j}(x,t)\\}$.

The holomorphic differentials $\omega_{j}$ are formally the form

$\omega_{j}=\sum_{k=1}^{\mathfrak{g}}\frac{C_{jk}\mu^{\mathfrak{g}-k}{\rm
d}\mu}{\sqrt{\mathcal{R}_{2\mathfrak{g}+2}(\mu)}},\qquad
j=1,2,\ldots\mathfrak{g}.$ (4.29)

Coefficients $C_{jk}$ are determined by requiring canonical normalisation with
respect to $\mathcal{A}$-cycles

$\oint_{\mathcal{A}_{j}}\omega_{k}=\delta_{jk}.$ (4.30)

Taking into account Eq. (4.28), equations of motion of the dynamical divisor
(4.23) linearise. Indeed, making use of the Lagrange interpolation formulae,
one can find

$\partial_{x}\varphi(x,t)=2\pi{\rm
i}\,C_{j1},\qquad\partial_{t}\varphi(x,t)=2\pi{\rm
i}\left[\frac{r_{1}}{2}C_{j1}+C_{j2}\right],$ (4.31)

implying

$\varphi_{j}(x,t)=k_{j}\,x+w_{j}\,t+\varphi_{j}(0,0),$ (4.32)

with wave numbers $k_{j}$ are frequencies $w_{j}$ reading

$k_{j}=2\pi{\rm i}C_{j1},\qquad w_{j}=2\pi{\rm
i}\Big{[}\frac{r_{1}}{2}C_{j1}+C_{j2}\Big{]}.$ (4.33)

Phases $\varphi_{j}(x,t)$ satisfy linear evolution in both space and time,
exercising a quasiperiodic motion on a Liouville torus
$\mathbb{T}^{2\mathfrak{g}}$ of real dimension $2\mathfrak{g}$. For any
physically admissible initial condition, $\gamma$-variables evolve along
closed trajectories which are homotopically equivalent to $\mathcal{A}$-cycles
$\mathcal{A}_{j}$, such that there is precisely one variable per cycle
$\mathcal{C}_{j}$.

Figure 2: Cycles on a two-sheeted Riemann surface, illustrated on an example
of two cuts $\mathcal{C}_{j}$ and $\mathcal{C}_{k}$. Upon tunnelling to the
other Riemann sheet, the integration orientation is reversed [12].

##### Periodicity constraints.

A family of periodic (i.e. closed) solutions
$\vec{\mathcal{S}}(x)=\vec{\mathcal{S}}(x+\ell)$ is further distinguished by
the periodicity of angles,
$\varphi_{j}(x+\ell,t)=\varphi_{j}(x,t)+2\pi\,n_{j}$ where integers
$n_{j}\in\mathbb{Z}$ specify the mode numbers assigned to each branch cut.
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,
imposing a non-trivial restriction on the admissible algebraic curves.

##### Extra degree of freedom.

Now we come to an important subtlety in the above construction: the dynamical
separated variables forming the divisor
$\mathcal{D}\equiv\\{\gamma_{j}\\}_{j=1}^{\mathfrak{g}}$ _do not_ provide the
complete set of dynamical degrees of freedom for the class of finite-gap
solutions. It turns out there is an additional degree of freedom that is not
amongst the canonical angle variables (4.28) obtained from dynamical zeros
$\gamma_{j}$ of the $b$-functions. This fact indeed becomes apparent already
in the simplest case of genus-zero solutions which, as we in turn demonstrate,
governs an evolution of a spin field described by a single angle variable.
More generally, for algebraic curves of genus $\mathfrak{g}$, there are thus
$(\mathfrak{g}+1)$ phases in total (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(\mathfrak{g}+1)$.

Substituting the form of the $b$-functions (4.21), into Eq. (4.18), and using
the restoration formulae (4.25), (4.26), (4.27), we deduce the following
identities

$\int^{\ell}_{0}{\rm
d}x\frac{\bar{b}(\mu;x,t)\mathcal{S}^{-}(x,t)}{\sqrt{\mathcal{R}(\mu)}+a(\mu)}=\int^{\ell}_{0}{\rm
d}x\frac{b(\mu;x,t)\mathcal{S}^{+}(x,t)}{\sqrt{\mathcal{R}(\mu)}+a(\mu)}.$
(4.34)

Introducing an auxiliary function

$r(\mu,\gamma)=\frac{\sqrt{\mathcal{R}(\mu)}-\sqrt{\mathcal{R}(\gamma)}}{\mu-\gamma}-\frac{1}{2}\sum_{\sigma\in\\{+,-\\}}\frac{\sqrt{\mathcal{R}(\sigma{\rm
i}\epsilon)}-\sqrt{\mathcal{R}(\gamma)}}{\sigma{\rm i}\epsilon-\gamma},$
(4.35)

we can present the quasi-momentum as a period integral of the form

$\mathfrak{p}(\mu)=-\frac{1}{2}\int^{\ell}_{0}{\rm
d}x\,\sum_{j=1}^{\mathfrak{g}}\frac{r(\mu,\gamma_{j})}{\prod_{k\neq
j}(\gamma_{j}-\gamma_{k})}.$ (4.36)

When $\mathfrak{g}=0$, we put Eqs. (4.21) and (4.22) into Eq. (4.18), which
yields

$\mathfrak{p}(\mu)=-\frac{\ell}{2}\left(\sqrt{\mathcal{R}_{2}(\mu)}-\frac{1}{2}\Big{(}\sqrt{\mathcal{R}_{2}({\rm
i}\epsilon)}+\sqrt{\mathcal{R}_{2}(-{\rm i}\epsilon)}\Big{)}\right).$ (4.37)

Comparing Eq. (4.36) with Eq. (4.26), we find that

$n_{\mathfrak{g}+1}\equiv\frac{1}{2\pi{\rm i}}\int^{\ell}_{0}{\rm
d}x\,\partial_{x}\log\mathcal{S}^{-},$ (4.38)

provides an extra mode number $n_{\mathfrak{g}+1}$.

Using the fact that $\gamma$-cycles are topologically equivalent to
$\mathcal{A}$-cycles that encircle the cuts of $\Sigma$, we can write

$\mathfrak{p}(\mu)=\frac{{\rm
i}}{2}\sum_{k=1}^{\mathfrak{g}}n_{k}\oint_{\mathcal{A}_{k}}\frac{\sqrt{\mathcal{R}(\mu)}-\sqrt{\mathcal{R}(\gamma)}}{\mu-\gamma}\frac{{\rm
d}\gamma}{\sqrt{\mathcal{R}(\gamma)}}-\pi n_{\mathfrak{g}+1},$ (4.39)

where $n_{k}$ is the mode number of the $k$th cut.

On Riemann surface $\Sigma$, quasi-momentum $\mathfrak{p}(\mu)$ is a single-
valued function. Alternatively, it can be understood as a double-valued
function on the complex spectral plane with $\mu\in\mathbb{C}$, experiencing
jumps on different Riemann sheets upon traversing branch cuts with
$\sqrt{\mathcal{R}(\mu)}$ flipping sign. The discontinuity condition for each
cut reads

$\mathfrak{p}(\mu+{\rm i}0)+\mathfrak{p}(\mu-{\rm
i}0)=2\pi(n_{k}-n_{\mathfrak{g}+1}),\qquad\mu\in\mathcal{C}_{k},$ (4.40)

where infinitesimal shifts $\pm{\rm i}0$ pertain to points on the either side
of the cut $\mathcal{C}_{k}$. The above equations are known as the
_Riemann–Hilbert problem_. An additional mode number $n_{\mathfrak{g}+1}$ can
be inferred from the asymptotic condition at $\mu=\infty_{\pm}$, reading

$\mathfrak{p}(\infty_{+})+\mathfrak{p}(\infty_{-})=-2\pi n_{\mathfrak{g}+1}.$
(4.41)

It proves convenient to introduce a new basis of open $\mathcal{B}$-cycles
$\mathcal{B}_{k}$ (for $k=1,2,\ldots,\mathfrak{g}+1$), connecting infinity
$\infty_{+}$ on the upper sheet with $\infty_{-}$ on the lower sheet by
passing through the cut $\mathcal{C}_{k}$, as demonstrated in Fig. 2, Eqs.
(4.40) alongside (4.41) can be compactly stated as

$\oint_{\mathcal{B}_{k}}{\rm d}\mathfrak{p}(\mu)=2\pi n_{k},\qquad
k=1,2,\ldots,\mathfrak{g}+1.$ (4.42)

##### An alternative approach.

Despite Eq. (4.39) provides a solution to the Riemann–Hilbert problem (4.40),
the behaviour at $\mu\to\infty$ does not comply with the form of Eq. (4.13).
An extra requirement on the spectral curve is needed, namely demanding
integrality of the wave numbers (4.33). One route to achieve this is to
construct quasi-momentum $\mathfrak{p}(\mu)$ that satisfies the required
asymptotics at $\mu\to\infty$ by considering a meromorphic (second-kind)
differential

${\rm
d}\mathfrak{p}=-\frac{\ell}{2}\frac{\mathcal{P}_{2\mathfrak{g}+2}(\mu){\rm
d}\mu}{\sqrt{\mathcal{R}(\mu)}},\qquad\mathcal{P}_{2\mathfrak{g}+2}(\mu)=\sum_{j=0}^{\mathfrak{g}+1}c_{j}\mu^{j}.$
(4.43)

Then $\mathfrak{p}(\mu)$ is unambiguously determined by demanding analyticity
and specifying its asymptotics; the latter readily fixes the highest two
coefficients of $\mathcal{P}_{2\mathfrak{g}+2}(\mu)$, namely

$c_{\mathfrak{g}+1}=1,\qquad
c_{\mathfrak{g}}=-\frac{r_{1}}{2}=-\frac{1}{2}\sum_{j=1}^{2\mathfrak{g}+2}\mu_{j},$
(4.44)

whereas the remaining $m=\mathfrak{g}-1$ coefficients are fixed by

$\oint_{\mathcal{A}_{j}}{\rm d}\mathfrak{p}(\mu)=0,\quad
j=1,2,\cdots\mathfrak{g}.$ (4.45)

Finally, the spectral curve is uniquely fixed by the condition (4.42).

##  5 Examples of finite-gap solutions

###  5.1 One-cut rational solution

The simplest solutions of the Riemann-Hilbert problem (3.23) belong to
algebraic curves of genus zero (Riemann surfaces with a single branch cut).
These correspond to quadratic polynomials of the form

$\mathcal{R}_{2}(\mu)=(\mu-\mu_{1})(\mu-\bar{\mu}_{1}),$ (5.1)

where branch points $\mu_{1},\bar{\mu}_{1}\in\mathbb{C}$ are conjugate to one
another in order to obey the reality condition. This leads to solutions that
involve two real degrees of freedom.

In what follows, we set the classical period to $\ell=1$. With this choice,
the admissible values of the wave numbers are $k=2\pi\,n$ with $n$ integer. As
a consequence, the branch points cannot be chosen arbitrarily but get
“quantised” as well.

Quasi-momentum $\mathfrak{p}(\mu)$ is given by Eq. (4.37). To satisfy the
“quantisation condition” and to obtain the prescribed filling fraction (4.20),
we have

$\frac{1}{2}\left(\sqrt{\mathcal{R}_{2}({\rm
i}\epsilon)}+\sqrt{\mathcal{R}_{2}(-{\rm i}\epsilon)}\right)=2\pi
n,\quad-\frac{1}{4}\left(\sqrt{\mathcal{R}_{2}({\rm
i}\epsilon)}-\sqrt{\mathcal{R}_{2}(-{\rm i}\epsilon)}\right)=-\frac{{\rm
i}\epsilon}{2}(1-2\nu),$ (5.2)

allowing us to parametrise the branch points as

$\mu_{1}+\bar{\mu}_{1}=4\pi
n(1-2\nu),\quad|\mu_{1}|^{2}=4\pi^{2}n^{2}+4\delta\nu(1-\nu).$ (5.3)

The algebraic curve can be expressed in terms of mode number $n$ and filling
fraction $\nu$, reading

$y^{2}(\mu)=\mathcal{R}_{2}(\mu)=\mu^{2}-4\pi
n(1-2\nu)\mu+4\pi^{2}n^{2}+4\delta(1-\nu),$ (5.4)

whereas the associated quasi-momentum, cf. Eq. (4.37), can be written as

$\begin{split}\mathfrak{p}_{1}(\mu)&=-\pi
n-\frac{1}{2}\sqrt{(\mu-\mu_{1})(\mu-\bar{\mu}_{1})}\\\ &=-\pi
n-\frac{1}{2}\sqrt{\mu^{2}-4\pi
n(1-2\nu)\mu+4\pi^{2}n^{2}+4\delta(1-\nu)}.\end{split}$ (5.5)

On the other hand, the quasi-momentum can also be constructed as an integral
of a suitable meromorphic differential on the rational curve (cf. Eq. (4.43)).
In $\mathfrak{g}=0$ case, such differential is described with coefficients
$c_{0}$ and $c_{1}$,

$\frac{{\rm d}\mathfrak{p}_{1}(\mu)}{{\rm
d}\mu}=-\frac{1}{2}\frac{c_{1}\mu+c_{0}}{\sqrt{\mathcal{R}_{2}(\mu)}}.$ (5.6)

In this case, both coefficients are readily fixed by the asymptotics, see Eq.
(4.44), yielding $c_{1}=1$ and $c_{0}=-(\mu_{1}+\bar{\mu}_{1})/2$. By taking
an integral of ${\rm d}\mathfrak{p}_{1}(\mu)$, we correctly recover the quasi-
momentum (5.5).

Notice that presently (i.e. in the zero-genus case) there is no canonical
$\mathcal{A}$-cycle. There is a single wave number which can be retrieved by
evaluating the ‘open’ $\mathcal{B}$-cycle

$k=\oint_{\mathcal{B}_{1}}{\rm
d}\mathfrak{p}_{1}=-\left(\mathfrak{p}_{1}(\infty_{+})+\mathfrak{p}_{1}(\infty_{-})\right)=2\pi
n.$ (5.7)

Here $n\in\mathbb{Z}$ is the mode number of the solution.

Coefficients of the asymptotic expansion of $\mathfrak{p}_{1}$ (cf. Eq.
(4.13))

$\mu\to\infty:\qquad\mathfrak{p}_{1}(\mu)\sim-\frac{\mu}{2}-\frac{\mathcal{P}}{2}-\frac{\mathcal{H}}{\mu}+\mathcal{O}\left(\mu^{-2}\right),$
(5.8)

provide phase-space averages of local charges evaluated on a particular one-
cut solution. The initial two coefficients correspond to total momentum and
energy, reading explicitly

$\mathcal{P}=2\pi
n\nu,\qquad\mathcal{H}=\frac{1}{2}(4\pi^{2}n^{2}+\delta)\nu(1-\nu).$ (5.9)

The knowledge of the one-cut quasi-momentum allows to express the dynamics of
the spin field $\vec{\mathcal{S}}(x,t)$ in terms of canonical separated
$\gamma$-variables. However, since we deal with an algebraic curve of genus
zero, we remind that there is no genuine $\gamma$-variable present. Instead,
the only dynamical degree of freedom is the transversal spin component
$\mathcal{S}^{-}(x,t)$. We first use Eq. (4.25) to deduce the longitudinal
component

$\mathcal{S}^{\rm
z}(x,t)=\frac{\sqrt{\mathcal{R}_{2}(\gamma_{1})}}{(\gamma_{1}-\gamma_{2})}+\frac{\sqrt{\mathcal{R}_{2}(\gamma_{2})}}{(\gamma_{2}-\gamma_{1})}=1-2\nu,$
(5.10)

where we have made use of the frozen (i.e. non-dynamical) $\gamma$-variables
$\gamma_{1}={\rm i}\epsilon$ and $\gamma_{2}=-{\rm i}\epsilon$.

In the next step, we can solve Eqs. (4.26) and (4.27) to find the transversal
component of the spin field,

${\rm i}\partial_{x}\log\mathcal{S}^{-}=2\pi n,\qquad{\rm
i}\left(\partial_{t}\log\mathcal{S}^{-}-\pi
n(1-2\nu)\partial_{x}\log\mathcal{S}^{-}\right)=\delta(1-2\nu).$ (5.11)

By imposing normalisation constraint $|\vec{\mathcal{S}}(x,t)|=1$, we finally
arrive at the following general form of the one-cut solution

$\mathcal{S}^{\rm
z}(x,t)=1-2\nu=\cos\theta_{0},\qquad\mathcal{S}^{\pm}(x,t)=\sin\theta_{0}\exp\left[\pm{\rm
i}(kx-wt)\right],$ (5.12)

where the wave number and frequency are

$k=2\pi n,\quad w=(4\pi^{2}n^{2}+\delta)\cos\theta_{0}.$ (5.13)

The momentum and energy can be alternatively computed by direct integration
using Eq. (4.13)

$\begin{split}&\mathcal{P}=\frac{1}{2{\rm i}}\int^{1}_{0}{\rm
d}x\,\frac{\mathcal{S}^{-}\mathcal{S}^{+}_{x}-\mathcal{S}^{+}\mathcal{S}^{-}_{x}}{1+\mathcal{S}^{\rm
z}}=2\pi n\nu,\\\ &\mathcal{H}=\frac{1}{2}\int^{1}_{0}{\rm
d}x\,\left[\vec{\mathcal{S}}_{x}\cdot\vec{\mathcal{S}}_{x}+\delta(1-(\mathcal{S}^{\rm
z})^{2})\right]=\frac{1}{2}(4\pi^{2}n^{2}+\delta)\nu(1-\nu),\end{split}$
(5.14)

in agreement with Eq. (5.9).

####  5.1.1 One-cut solution from asymptotic Bethe ansatz

(a)

(b)

Figure 3: (a) Comparison between momentum of a one-cut solution with mode
number $n=1$ (red dashed line, for both the isotropic and anisotropic
interaction $\delta$) and momenta of the corresponding quantum eigenstates
(shown for different system sizes $L=M/\nu$, with $M=30$ and filling fraction
$\nu$), with green (black) circles representing the isotropic (anisotropic,
with $\delta=1$) cases. (b) Comparison between energies of a one-cut solution
with mode number $n=1$ (red dashed line for isotropic case and blue dashed
line for anisotropic case with $\delta=1$) and rescaled energies ($E\cdot L$)
of the corresponding quantum eigenstates (the same system sizes as in panel
(a)).

Before describing how to perform semi-classical quantisation of finite-gap
solutions (cf. Sec. 6), we wish to briefly discuss the one-cut solutions from
the vantage point of low-energy eigenstates in the Heisenberg anisotropic spin
chain. Our aim is mainly to identify which solutions to the Bethe equations
(3.3) show up classically as one-cut solutions. A numerical method for
achieving this is outlined in Appendix E.1. Specifically, we are seeking for
low-energy solutions comprising of macroscopically $\mathcal{O}(L)$ excited
magnons that condense into a single coherent mode. By fixing the mode number
to $n=1$, we can compute the total momentum $P$ and total energy $E$ of such
state with different system sizes $L$ and filling fraction $\nu=M/L$ fixed,
see Eqs. (3.5). The results are shown in Fig. 3. On the other hand, we can
likewise compute the momenta and energies for classical one-cut solutions, cf.
Eq. (5.9), by setting the classical period $\ell=1$.

Despite that the quantum quasi-momentum in the classical limit reduces to its
classical counterpart, we remind that a direct comparison of the respective
energies requires to account for an extra rescaling factor of $L$, namely
$H\cdot L\sim\mathcal{H}$, where $H$ denotes the eigenvalue of the quantum
Hamiltonian defined in Eq. (3.1), while the classical Hamiltonian
$\mathcal{H}$ is given by Eq. (3.25).

At the classical level, one-cut solutions describe a simple precessional
motion. In some sense, one may view them as a finite-density analogue of
Goldstone modes (representing elementary ferromagnetic excitations called
magnons). Remarkably, it turns out (see Sec. 6) that quantisation of such
solutions can give rise to certain non-perturbative effects (first discussed
in Ref. [13]) which necessitate to incorporate quantum corrections into the
picture.

###  5.2 Bions from degenerate two-cut solutions

As a non-trivial example, we next focus on the class of two-cut solutions.
These are, physically speaking, periodic elliptic magnetisation waves which
are often referred to as ‘cnoidal waves’. The space of two-cut solution is
associated to elliptic algebraic curves corresponding to Riemann surfaces of
genus $\mathfrak{g}=1$, characterised by two branch cuts. In the special limit
of unit elliptic modulus, the profiles become trigonometric and one retrieves
the celebrated soliton solutions. For the particular case of the easy-axis
anisotropic Landau–Lifshitz model, there exist special types of two-cut
solutions known as _bions_ ; these represent two-mode bound states formed of a
kink and an anti-kink. Moreover, a special degeneration of such a bion
solution (upon decompactifying the circumference) produces a static kink. As
mentioned in the introduction, kinks are ultimately responsible for the
observed freezing of a magnetic domain wall, as demonstrated in [16]. In Sec.
6.3 below, we shall perform the semi-classical quantisation on a bion and
study its subtle features. Before that, we derive it here using the outlined
integration procedure.

The elliptic curve encoding the bion spectrum has the form

$\mathcal{R}_{\rm
bion}^{2}(\mu)=\mathcal{R}_{4}^{2}(\mu)=(\mu^{2}+\xi^{2}_{1})(\mu^{2}+\xi^{2}_{2}),$
(5.15)

parametrised by two pairs of complex-conjugate branch points located on the
imaginary axis ${\rm Re}(\mu)=0$ at $\mu_{j}\in\\{\pm{\rm i}\xi_{1},\pm{\rm
i}\xi_{2}\\}$, satisfying conditions

$\xi_{1}>\epsilon>0,\qquad 0<\xi_{2}<\epsilon.$ (5.16)

The upshot here is that the bion solutions can only be found in the regime
$\delta>0$ (the easy-axis regime). In what follows we put, mostly for
simplicity, $\delta=\epsilon=1$. In fact, from the classical equation of
motion for the spin field $\vec{\mathcal{S}}(x,t)$ given by Eq. (3.24), the
solution at $\delta=1$ can be mapped to another solution with
$\delta^{\prime}>0$ by a simple rescaling

$\vec{\mathcal{S}}(x,t)\quad\mapsto\quad\vec{\mathcal{S}}(x^{\prime}=\epsilon\,x,t^{\prime}=\delta\,t).$
(5.17)

We next outline how to reconstruct of the spin field from the algebraic curve.
To set the stage, we introduce the standard full elliptic integrals, cf.
Appendix C,

$K_{1}=K\left(\frac{\xi_{2}^{2}}{\xi_{1}^{2}}\right),\qquad
K_{2}=K\left(1-\frac{\xi_{2}^{2}}{\xi_{1}^{2}}\right).$ (5.18)

When the argument of the elliptic function is omitted we adopt that ${\rm
k}=\xi_{2}/\xi_{1}$.

Since the Riemann surface involves two branch cuts, this time we do have a
genuine dynamical $\gamma$-variable $\gamma_{1}(x,t)$. As said earlier, there
exist two extra non-dynamical variables pinned to locations $\gamma_{2}={\rm
i}$ and $\gamma_{3}=-{\rm i}$ (recall that $\epsilon=1$). From the Dubrovin
equations (4.23) we find

$\gamma_{1}(x,t)=\frac{{\rm i}\xi_{1}}{{\rm sn}(u)}=-{\rm i}\xi_{2}\,{\rm
sn}(u+{\rm i}K_{2}),\qquad u=\xi_{1}x+\varrho,$ (5.19)

where we have used ${\rm sn}(x+{\rm i}K_{2}){\rm sn}(x){\rm k}=1$, and used
$\varrho$ to denote the integration constant. Remarkably, it turns out that in
this particular subclass of two-cut solutions even $\gamma_{1}$ is static.

Using the reconstruction formulae (4.25) for $\mathcal{S}^{\rm z}$ component,
we thus have

$\mathcal{S}^{\rm z}=\frac{\sqrt{\mathcal{R}_{4}(\gamma_{1})}+{\rm
i}\sqrt{\mathcal{R}_{4}({\rm i})}\gamma_{1}}{\gamma_{1}^{2}+1}.$ (5.20)

To fix the integration constant $\varrho$, we impose the reality condition
$\mathcal{S}^{\rm z}\in\mathbb{R}$,

$\varrho=\mathrm{arcsn}\left({\rm
i}\frac{\xi_{1}}{\xi_{2}}\sqrt{\frac{1-\xi_{2}^{2}}{\xi_{1}^{2}-1}}\right),$
(5.21)

which yields a time-independent profile

$\mathcal{S}^{\rm
z}(x)=-\xi_{2}\frac{\mathrm{cn}(\xi_{1}x)}{\mathrm{dn}(\xi_{1}x)}.$ (5.22)

The solution has a spatial period

$\ell=\frac{4nK_{1}}{\xi_{1}},$ (5.23)

where $\pm n$ are the mode numbers associated with the two cuts. Variable
$\gamma_{1}$ can be therefore expressed as

$\gamma_{1}(x)=-\xi_{1}\xi_{2}\frac{-{\rm i}\xi_{1}\sqrt{\mathcal{R}_{4}({\rm
i})}\,\mathrm{cn}(\xi_{1}x)\mathrm{dn}(\xi_{1}x)+{\rm
i}(\xi_{1}^{2}-\xi_{2}^{2}){\rm
sn}(\xi_{1}x)}{\xi_{1}^{2}(1-\xi_{2}^{2})+(\xi_{1}^{2}-1)\xi_{2}^{2}\,{\rm
sn}(\xi_{1}x)^{2}}.$ (5.24)

The other independent component of the spin field, say $\mathcal{S}^{-}(x,t)$,
can be reconstructed with aid of formulae (4.26) and (4.27), i.e.

$\mathcal{S}^{-}(x,t)=\frac{1}{\sqrt{1+\gamma_{1}^{2}(x,t)}}\exp\left(-{\rm
i}\int{\rm d}x\frac{-{\rm i}\sqrt{\mathcal{R}_{4}({\rm
i})}}{1+\gamma_{1}^{2}(x,t)}\right).$ (5.25)

Using the properties of elliptic functions, we have

$\mathcal{S}^{-}(x)=C\frac{\Theta(\xi_{1}x+\varrho+{\rm
i}K_{2})}{\Theta(\xi_{1}x+K_{1})}\exp\left(\xi_{1}x\frac{\pi{\rm
i}}{2K_{1}}+\xi x\frac{\Theta^{\prime}(\beta)}{\Theta(\beta)}\right),$ (5.26)

where an auxiliary function $\Theta(x)$ is defined as (cf. Appendix C)

$\Theta(x)=\vartheta_{3}\left(\frac{\pi
x}{2K_{1}}+\frac{\pi}{2},-i\frac{K_{2}}{K_{1}}\right).$ (5.27)

Finally, constant $C$ is can be fixed by requiring normalisation
$\vec{\mathcal{S}}^{2}=1$, yielding

$\mathcal{S}^{-}(x)=\sqrt{1-\xi_{2}^{2}}\frac{\Theta(\xi_{1}x+\rho+{\rm
i}K_{2})}{\Theta(\xi_{1}x+K_{1})}\exp\left(\xi_{1}x\frac{\pi{\rm
i}}{2K_{1}}+\xi x\frac{\Theta^{\prime}(\beta)}{\Theta(\beta)}+{\rm
i}\phi_{0}\right),$ (5.28)

where $\phi_{0}\in\mathbb{R}$ is a phase that is determined by the initial
condition. We just found that for the bion solution $\mathcal{S}^{-}$ is time-
independent as well. A representative example of a bion solution is shown in
Fig. 4.

Figure 4: Spin-field components of a typical bion solution, depicted for
anisotropy parameter $\delta=1$ and branch points $\xi_{1}\simeq 1.0583559$
and $\xi_{2}=0.9$. Components $\mathcal{S}^{\rm z}(x)$, $\mathcal{S}^{\rm
x}(x)=\mathrm{Re}\,\mathcal{S}^{-}(x)$, and $\mathcal{S}^{\rm
y}(x)=-\mathrm{Im}\,\mathcal{S}^{-}(x)$ are shown by blue, yellow and red
curves, respectively. The branch points are determined using the method in
Appendix E. The bion solution is periodic with period
$\ell_{2}=4K_{1}/\xi_{1}\simeq 15.956517$.

##### Static kink.

As our final example, we consider a particular degeneration of a bion solution
which produces a kink. One can think of it as ‘soliton limit’ which generally
corresponds to ‘blowing up’ the period $\ell$. This requires to bring both
branch points of $\sqrt{\mathcal{R}_{4}(\mu)}$ in the upper-half plane
together to meet at ${\rm i}\epsilon$, that is sending $\xi_{1,2}\to\epsilon$
(presently $\epsilon=1$) and thus collapsing both branch cuts to a point. In
order to perform this soliton limit, we first shift the argument of the
elliptic function by quarter period $K_{1}/\xi_{1}$; for instance the
$\mathcal{S}^{\rm z}$ field takes the form (cf. Eqs. (C.6))

$\mathcal{S}^{\rm
z}(x)=-\xi_{1}\frac{\mathrm{cn}(\xi_{1}(x+K_{1}/\xi_{1}))}{\mathrm{dn}(\xi_{1}(x+K_{1}/\xi_{1}))}=\xi_{1}\mathrm{sn}(\xi_{1}x).$
(5.29)

Now taking the limits $\xi_{1,2}\to 1$ and accordingly ${\rm
k}=\xi_{2}/\xi_{1}\to 1$, we find

$\mathcal{S}^{\rm z}_{\rm kink}(x)=\tanh(x),$ (5.30)

which is none other than a _static kink_! The transversal components can be
easily deduced from the equation of motion (3.24), yielding

$\mathcal{S}^{\rm-}_{\rm kink}(x)=\operatorname{sech}(x)\,e^{{\rm
i}\phi_{0}},$ (5.31)

where $\phi_{0}\in\mathbb{R}$ is a phase determined by the initial condition.

##  6 Semi-classical quantisation

We have now fully prepared the stage to carry out the semi-classical
quantisation on finite-gap solutions. In this respect, the associated spectral
curves (encoding complete information about the conserved quantities) will be
of vital importance.

An important remark is in order at this point. First, recall that moduli of
algebraic curves are completely fixed by locations of the square-root branch
points, i.e. roots of polynomial $\mathcal{R}_{2\mathfrak{g}+2}(\mu)$, which
(by the reality condition) must always appear in complex-conjugate pairs.
Local conserved charges are expressible as functions of symmetric polynomials
of branch points, namely coefficients $r_{j}$ of
$\mathcal{R}_{2\mathfrak{g}+2}(\mu)$. While branch points $\\{\mu_{j}\\}$
themselves are directly linked physical quantities, the branch cuts (of
square-root type) obtained by pairwise connecting the branch points are not
physical but merely a matter of convention. There is indeed plenty of freedom
in assigning branch cuts to a given set of branch points. For instance, the
prevalent choice in the finite-gap literature [36] is to place the cuts along
straight vertical lines connecting every complex conjugate pair of branch
points, which are in turn encircled by the canonical basis
$\mathcal{A}$-cycles. Such a choice is, purely from the standpoint of
classical finite-gap solutions, perfectly adequate. One the other hand, if
classical solutions are instead thought of as emergent macroscopic bound
states of magnons of the underlying quantum spin chain, it is natural to cut
the surface along one-dimensional disjoint segments associated to forbidden
zones of the classical transfer function [12], corresponding to contours on
which magnons (Bethe roots) condense. This prescription for the branch cuts is
physically distinguished. As demonstrated in turn, such physical cuts not only
appreciably differ from the conventional straight cuts in general, but also
undergo the phenomenon of condensate formation. Computing the Bethe root
densities along the physical contours is thus the essential step to perform
the semi-classical quantisation.

###  6.1 Physical contours

In this section we describe a general procedure to determine the physical
contours. The algorithm we employ has been proposed and implemented in Ref.
[13]. We shall also facilitate a direct comparison with exact low-momentum
quantum eigenstates at finite $L$ in the weakly-anisotropic regime. For our
convenience, we carry out this computation in $\zeta$-plane 555Our variable
$\zeta$ is analogous to variable $x$ used in Refs. [12, 55, 13] in the case of
the isotropic Heisenberg model. by applying the following anti-holomorphic
transformation 666Upon this transformation, orientation of integration
contours gets reversed.

$\mu\mapsto\zeta(\mu):\qquad\zeta=\frac{1}{\mu}=\frac{\tan\vartheta}{\epsilon}.$
(6.1)

Firstly, we introduce the _complex_ density function

$\rho(\zeta)=\frac{\mathfrak{p}(\zeta+{\rm i}0)-\mathfrak{p}(\zeta-{\rm
i}0)}{2{\rm i}\pi\ell(1+\delta\zeta^{2})}=\pm\frac{\mathfrak{p}(\zeta\pm{\rm
i}0)-\pi n_{j}}{{\rm
i}\pi\ell(1+\delta\zeta^{2})},\qquad\zeta\in\mathcal{C}_{j},$ (6.2)

where $n_{j}\in\mathbb{Z}$ is the mode number associated to cut
$\mathcal{C}_{j}$. By virtue of the second equality in Eqs. (6.2), the density
function $\rho(\zeta)$ can be defined on the entire Riemann surface. As a
direct consequence of $\mathfrak{p}(\zeta=\zeta_{\star})=\pi n_{j}$ at the
square root branch points $\zeta_{\star}\in\\{\zeta_{j},\bar{\zeta}_{j}\\}$,
the density always vanishes, $\rho(\zeta_{\star})=0$. In a small neighbourhood
around it, one finds
$\rho(\zeta=\zeta_{\star}+\varepsilon)=\mathcal{O}(\sqrt{\varepsilon})$. Away
from the branch points $\rho(\zeta)$ in general takes complex values.

The task is to determine the physical contours $\mathcal{C}_{j}$ The latter
can be singled out by the following condition: starting from the branch point
$\zeta_{j}$, the integrated density differential $\rho(\zeta){\rm d}\zeta$
must always remain real,

$\int^{\zeta}_{\zeta_{j}}{\rm
d}\zeta^{\prime}\rho(\zeta^{\prime})\in\mathbb{R}\qquad{\rm
for}\qquad\zeta\in\mathcal{C}_{j}.$ (6.3)

This prescription has a transparent physical interpretation; physically
$\rho(\zeta){\rm d}\zeta$ corresponds to the number of excitations (Bethe
roots) within an infinitesimal interval in $\zeta$-plane, which is a positive
definite quantity by definition. This condition alone is however not
sufficient yet. In particular, it turns out that there are three distinct
contours emanating out of each branch point compatible with the above
positivity requirement [13]. Amongst those, one of the contour carries an
infinite filling fraction and can be thus immediately ruled out as unphysical.
Out of the remaining two contours, only one can be physical. The defining
condition is that the total filling fraction does not exceed the threshold
value of maximal total filling $\nu_{\rm max}=1/2$, that is

$\int_{\mathcal{C}}{\rm d}\zeta^{\prime}\rho(\zeta^{\prime})\leq\nu_{\rm
max},\qquad\mathcal{C}=\bigcup_{j}\mathcal{C}_{j}.$ (6.4)

Consider now a certain reference finite-gap solution. To quantise it at the
semi-classical level, every density contour (physical cut) has to be dissolved
into a large (but finite) number of individual magnons. This invariably
requires to reintroduce the length $L$ of the underlying spin chain, thus
rendering the total magnetisation carried by individual coherent states to
come in integer quanta of $M_{j}\sim\mathcal{O}(L)$. The precise requirements
are that (i) $M_{j}/L\approx\nu_{j}$ and (ii) that the Bethe roots are
distributed approximately with density $\rho_{j}(\zeta)$ along
$\mathcal{C}_{j}$. Here it is important to make a distinction with the exact
quantisation which instead takes quantum fluctuations fully into account to
all orders in the effective Planck constant. This means, in other words, that
the semi-classical solutions produced with the outlined procedure can be at
best an approximation of finite-volume exact quantum-mechanical eigenstates at
large wavelengths, while a full non-perturbative (i.e. exact) quantisation
would require solving the Bethe ansatz equations (3.3).

##### Single contour at low density.

To benchmark the above procedure, we proceed by illustrating first how one-cut
solutions emerge as semi-classical eigenstates in the anisotropic gapped
Heisenberg spin-$1/2$ chain.

We shall first suppose that the filling fraction of a physical cut is
sufficiently _low_ , ensuring that the finite-size effects (cf. Eqs. (3.18)
and (B.16)) can be safely neglected at large system sizes. We then find that
the Bethe roots patterns which solve the asymptotic Bethe equations distribute
along certain arcs in the complex rapidity plane, as exemplified in Fig. 5. To
be concrete, we put anisotropy to $\delta=1$ and set the filling fraction to
$\nu=0.1$ and the mode number to $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 procedure to numerically solve the Bethe ansatz equations (3.3) for the
case of a single quantised one-cut solution is described in Appendix E.1).
Upon taking the $L\to\infty$ limit and rescaling the rapidity variable, the
semi-classical eigenstate will eventually be described by a dense arrangement
of Bethe roots distributing along the contour specified by the conditions
(6.3) and (6.4).

By ramping up the filling fraction $\nu$, we observe that ‘quantum
fluctuations’ (contained in higher order terms in the ABE) progressive
amplify. As announced earlier, this eventually leads to a critical phenomenon
of condensate formation. This feature will be closely examined in the next
section.

Figure 5: Direct comparison between the physical density contour of a
classical one-cut solution (red dashed line, corresponding to anisotropy
$\delta=1$, filling fraction $\nu=0.1$ and mode number $n=1$) determined by
imposing the positivity condition (6.3), and the corresponding solution to
Bethe ansatz equations (3.3) with $M=30$ Bethe roots, system size $L=300$ and
anisotropy $\eta=\sqrt{\delta}/L=1/300$ (blue dots). The Bethe roots are given
by $\zeta_{j}=\tan(\vartheta_{j})/\sqrt{\delta}$.

###  6.2 Formation of condensates

_Condensates_ refer to segments of a uniform density as a part of a physical
contour. We borrowed this name from Refs. [56, 12] where (to the best of our
knowledge) such objects have been first identified. Condensates enter the
picture whenever the maximal density along one of the physical contours
exceeds a particular critical value which is signalled by a divergence of the
finite-size correction given by Eq. (3.18) (see also Eq. (B.16)).

We shall first examine the phenomenon on the simplest case of the one-cut
solution, using the $\zeta$-plane parametrisation. By starting at some
sufficiently low filling fraction $\nu$ we can observed that upon gradually
increasing the filling fraction, the value of $\rho(\zeta)(1+\delta\zeta^{2})$
on real axis approaches the value of ${\rm i}$. This value is reached at the
critical filling fraction $\nu_{\rm crit}$, precisely when quantum fluctuation
of order $\mathcal{O}(1/L)$ become divergent, as indicated by Eq. (B.16). For
larger fillings $\nu>\nu_{\rm crit}$, the density contour develops a
vertically straight segment of unit uniform density. Such a condensate appears
first on the real axis (right after crossing $\nu_{\rm cric}$) and
progressively expands outwards upon further increasing $\nu$. From the
viewpoint of the underlying quantum chain, the spacing between constituent
Bethe roots is always equal to ${\rm i}\eta$. One can therefore think of
condensates as giant regular Bethe strings. In the complex spectral plane
associated to finite-gap solutions, the ends points of a condensate correspond
to branch cuts of _logarithmic_ type. 777Logarithmic branch cuts get likewise
produced in a well-known soliton degeneration process, corresponding to
merging two nearby standard square-root branch cuts by coalescing their type
branch points in a pairwise manner. In effect, finite-gap quasi-momentum is no
longer meromorphic. Condensates are different in this respect, as their quasi-
momentum differential remains meromorphic all the way through.

(a)

(b)

(c)

Figure 6: Branch points ($\zeta_{1}$ and $\bar{\zeta}_{1}$) and fluctuation
points ($\zeta_{F}$) for the one-cut solution with anisotropy parameter
$\delta=1$ and various filling fractions $\nu=\\{0.1,0.206354963,1/3\\}$,
increasing from left to right. Panel (b) corresponds to the critical filling
fraction $\nu_{c}\simeq 0.206354963$, when the fluctuation point $\zeta_{F}$
collide with the physical cut.

In spite of appearance of an additional condensate above $\nu>\nu_{\rm crit}$,
it is still possible to extract the physical contour solely from the knowledge
of a finite-gap solution by taking into account conditions (6.3) and (6.4).

Emergence of condensates is intimately tied to the notion of ‘fluctuation
points’, playing a pivotal role in classical modulation stability theory [57].
Fluctuation points can perceived as small fluctuations of a reference finite-
gap solution, corresponding to tiny cuts that possess infinitesimal filling
fractions. Upon increasing their filling fraction, they grow up into nonlinear
finite-gap mode. Let us slightly expand on this point. Imagine a reference
finite-gap solution with $m$ cuts, and let $\\{n_{1},n_{2}\cdots n_{m}\\}$
denote the occupied mode numbers. To excite a mode with an unoccupied $n$,
call it $n_{*}$, the following condition for the quasi-momentum
$\mathfrak{p}(\zeta)$ has to satisfy the periodic boundary condition,

$\mathfrak{p}(\zeta_{m,n_{*}})=n_{*}\pi,$ (6.5)

using $\zeta_{m,n_{*}}$ to 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 quasi-momentum).

Fluctuation points are treated as “almost degenerate” branch cuts, so small
that they do not affect the form of the quasi-momentum $\mathfrak{p}(\zeta)$.
This raises an interesting question whether classically any given finite-gap
solution is modulationally stable under such fluctuations, see for instance
discussions in Ref. [57]. In this respect, we note that the stability
condition coincides with the condition for the formation of condensates.

In the following we shall first take a look at the basic case with a single
cut. We find the physical contours made out of Bethe roots 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 middle’, with two
additional bent contours emanating from the intersection points that connect
to the nearby fluctuation point(s). We give an explicit demonstration of this
scenario in Sec. 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 Sec. 6.2.2.

Let us also mention that a somewhat reminiscent phenomenon is known to appear
in the context of matrix models [58] (which are described by a similar type of
Riemann-Hilbert problems) and also elsewhere, e.g. in the large-$N$ Yang–Mills
theory [59, 60, 61] and random tiling models [62, 63]. They commonly go under
the name of the Douglas–Kazakov phase transition [61], a variant of a third-
order phase transition. Analogous condensates also appear in the semi-
classical regime of the Lieb-Liniger model with attractive interaction [64,
65] where a quantum phase transition can be detected through the calculation
of correlation functions in the ground state [65]. We emphasise however that
in our case there is no real phase transition going on in the sense that
branch points and the quasi-momentum $\mathfrak{p}(\zeta)$ itself do not
undergo any discontinuous change, in distinction with the case of
Douglas–Kazakov transition where the free energy becomes non-smooth after
formation of a “condensate” [63].

####  6.2.1 One-cut case with condensate

In the case of a single cut, we have

$\int_{\mathcal{C}_{1}}{\rm d}\zeta\rho(\zeta)=\nu_{1}\in\mathcal{O}(1),$
(6.6)

with an upper-bounded filling fraction $\nu_{1}<1/2$. Locations of fluctuation
points, denote below by $\zeta_{1,k}$, can be inferred from the density

$\rho_{k}(\zeta)=\pm\frac{\mathfrak{p}(\zeta)-\pi k}{2\pi
i(1+\delta\zeta^{2})},\qquad\mathfrak{p}(\zeta_{1,k})=k\pi,$ (6.7)

##### Mode number $n=1$.

The condensate phenomenon can be best illustrated on the basic example of a
one-cut solution with mode number $n=1$. Below the density threshold we find a
single smooth arc-shaped contour, as exemplified in panel (a) in Fig. 6. The
closest fluctuation point sits at a finite distance away from the cut
(somewhere to the left of it), whereas the density at the real axis satisfies

$\rho(\zeta_{\ast})(1+\delta\zeta_{\ast}^{2})<{\rm i}.$ (6.8)

Increasing the filling fraction will cause an increase in the density on the
real axis. During the process, the nearby fluctuation point approaches the
physical contour until at the critical filling it eventually collides with it
at $\zeta_{\ast}$,

$\rho(\zeta_{\ast})(1+\delta\zeta_{\ast}^{2})={\rm i}.$ (6.9)

This event is shown in panel (b) in Fig. 6. 888Comparing to the isotropic case
in [13], $\nu_{c}\simeq 0.2092896452$, the condensate appears with a slightly
smaller filling fraction.

(a)

(b)

Figure 7: (a) Fluctuation point $\zeta_{F}$, lying on the real axis, can be
viewed as an infinitesimal (collapsed) branch cut – a deactivated mode. (b)
Comparison between the classical contour (dashed line) with anisotropic
parameter $\delta=1$, filling fraction $\nu=0.3$, and mode number $n=1$
(computed based on reality condition (6.3)) with a condensate (red dashed
line) and an additional contour originating from $\zeta_{F}$ (brown dashed
line) and the corresponding solution to the Bethe equations (cf. Eq. (3.3))
with $M=48$, $L=144$ and $\eta=\frac{\sqrt{\delta}}{L}=\frac{1}{144}$ (blue
dots). Notice that the Bethe roots plotted are
$\frac{\tan\vartheta}{\sqrt{\delta}}$.

Upon increasing the filling fraction even further, the fluctuation points
after collision ‘tunnel through’ the cut. This leaves behind a straight
condensate positioned in a vertical direction. The nearest fluctuation point
on the real axis will then appear to the right of the physical cut, as
pictured in panel (c) in Fig. 6. One can nonetheless recover the same quasi-
momentum $\mathfrak{p}(\zeta)$ by considering an additional pair of contours
which emanate out of the fluctuation point(s), satisfying $\rho_{2}(\zeta){\rm
d}\zeta\in\mathbb{R}$ with density defined through Eq. (6.7) (depicted by
brown dashed lines in panel (c) in Fig. 6). Due to an extra condensate, the
original contour cannot accommodate for all the Bethe roots, and some “excess”
Bethe roots will lie along those additional contours. We wish to emphasise
again that the quasi-momentum $\mathfrak{p}(\zeta)$ remains intact.

There is a suggestive explanation behind the above picture if one pictures a
one-cut solution as a limiting (degenerate) case of a more general two-cut
solution with one of its cuts ‘switched off’ to a fluctuation point. This is
neatly captured in Fig. 7 in panel (a), where the blue solid lines represent
parts of the original physical connecting to the pair branch points
$(\zeta_{1},\bar{\zeta}_{1})$, while the red dashed line depicts the Bethe-
root condensate of uniform density. The extra green solid line belongs to one
of the “unphysical contours”999We call it “unphysical” because the green
contour alone does not yield the correct value for the filling fraction for
the infinitesimal cut. Yet the combination of all three parts here is clearly
physically meaningful. associated with the infinitesimal branch cut
$(\zeta_{F},\bar{\zeta}_{F})$. Combining all the ingredients, we are therefore
able to determine the densities of Bethe roots along these three contours.
This amounts to account for the leading-order quantum corrections to ABE
(3.16) in non-perturbatively fashion.

To better corroborate the above picture, we made a direct comparison with the
contours obtained numerically by solving the Bethe equations for large system
sizes, cf. Fig. 7. Moreover, we have performed another quantitative test for
the proposed contours through the calculation of the leading-order Gaudin norm
for the Bethe states, as shown in Fig. 13. We emphasise that physical contours
are a key ingredient for the functional integral approach to compute overlaps
(and norms) between semi-classical Bethe eigenstates, thus providing an
opportunity to verify whether the described contours are indeed suitable. The
numerical evidence is collected in Appendix D.1.

##### Mode number $n\geq 2$.

One can encounter even more exotic situations. While
$\mathfrak{p}(\zeta_{\ast})=(n+1)\pi$ has only one real solution for $n=1$,
higher mode numbers $n\geq 2$ permit for complex-conjugate pairs of
fluctuation points [13]. In this scenario, the same condition
$\rho_{n+1}(\zeta){\rm d}\zeta\in\mathbb{R}$ yields an additional contour with
a condensate appearing between the intersection points, along the lines of the
proceeding discussion. This time instead, such contours can be understood at
the classical level as arising from a three-cut solution with one large
physical cut and two almost degenerate tiny cuts located at the complex-
conjugated fluctuation points $\zeta_{F}$ and $\bar{\zeta}_{F}$. For instance,
in Fig. 8 we give an illustration of that for the isotropic Heisenberg spin
chain with mode number $n=2$. Unfortunately, for the anisotropic ferromagnet
the employed numerical method for producing analogous solutions does not work
for $n\geq 2$, cf. Appendix E.1. Given that the distributions of Bethe roots
do not appreciably change upon introducing a tiny deformation parameter
$\eta\sim\mathcal{O}\left(\frac{1}{L}\right)$, we expect the phenomenon to
survive the presence of weak interaction anisotropy.

(a)

(b)

Figure 8: (Left) Complex fluctuation points $\zeta_{F}$ and $\bar{\zeta}_{F}$.
They can be seen as collapsed cuts of a three-cut solution with
$\zeta_{F,1}\to\zeta_{F,2}\to\zeta_{F}$ and
$\bar{\zeta}_{F,1}\to\bar{\zeta}_{F,2}\to\bar{\zeta}_{F}$. (Right) Comparison
between the classical density contour (dashed line), including with the
condensate (red dashed line) and the additional contours emanating from
fluctuation points $\zeta_{F}$ and $\bar{\zeta}_{F}$ (brown and purple dashed
lines), obtained from reality condition (6.3) for the case of isotropic
interaction ($\delta=0$), with filling fraction $\nu=0.1$ and mode number
$n=2$, to the corresponding solution to Bethe equations (3.3) with $M=60$,
$L=600$ and $\eta=0$. The Bethe roots $\zeta_{j}$ (blue dots) are rescaled by
the system size $L$ and plotted in the inverse spectral plane, i.e.
$\zeta_{j}=1/(L\lambda_{j})$, where $\lambda_{j}$ solve the isotropic Bethe
equations, $\left(\frac{\lambda_{j}+{\rm i}/2}{\lambda_{j}-{\rm
i}/2}\right)^{L}=\prod_{k\neq j}^{M}\frac{\lambda_{j}-\lambda_{k}+{\rm
i}}{\lambda_{j}-\lambda_{k}-{\rm i}}$.

####  6.2.2 Multiple cuts

When multiple cuts get involved, the situation is far more complicated. In
that case, the condensates can appear not only out of fluctuation points but
also via an attractive interaction amongst the physical cuts. Here we focus
for simplicity on the two-cut case, since a general scenario with more cuts
can be largely described based on the phenomenology of the two-cut case. In
Ref. [13], the authors made an exhaustive survey on the two-cut case at
isotropic point ($\delta=0$). The anisotropic case with
$\eta=\frac{\epsilon}{L}>0$ can be analysed in a similar fashion. There are
several discernible features we wish to highlight.

Firstly, when two physical cuts are far apart from one another, each branch
cut can produce a condensate upon colliding with their nearby fluctuation
points, in analogy with the one-cut case discussed in Ref. [13]. However, when
the physical cuts approach closer their mutual attraction amplifies until they
eventually merge with one another. The result of this are two joined contours
glued via a condensate at the cusps.

(a)

(b)

Figure 9: (Left) Comparison between the density contour (dashed line) of a
classical two-cut solution (shown for the isotropic case ($\delta=0$), with
partial filling fractions $\nu_{1}=0.02$, $\nu_{2}=0.06$ and mode numbers
$n_{1}=1$, $n_{2}=2$), obtained from reality condition (6.3), and the
corresponding numerical solution to Bethe equations (3.3) with $M_{1}=10$,
$M_{2}=30$, $L=500$ and $\eta=0$ (blue dots). (Right) Comparison between the
classical contour (dashed line) (shown for the isotropic case ($\delta=0$),
with partial filling fractions $\nu_{1}=0.025$, $\nu_{2}=0.075$ and mode
numbers $n_{1}=1$, $n_{2}=2$), obtained from reality condition (6.3), to the
corresponding numerical solution to Bethe equations (3.3) with $M1=10$,
$M_{2}=30$, $L=400$ and $\eta=0$ (blue dots). The Bethe roots plotted are
$\zeta_{j}=1/(L\lambda_{j})$, same as in Fig. 8.

A basic instance of the above phenomenon involves two cuts with consecutive
mode numbers, namely $n_{2}=n_{1}+1$. The moment the two physical contours
intersect, say at points $\zeta_{\rm int}$ and $\bar{\zeta}_{\rm int}$, the
combined density satisfies

$\left[\rho_{n_{1}}(\zeta_{\rm int})+\rho_{n_{1}+1}(\zeta_{\rm
int})\right]\left(1+\delta\zeta_{\rm
int}^{2}\right)=\left[\rho_{n_{1}}(\bar{\zeta}_{\rm
int})+\rho_{n_{1}+1}(\bar{\zeta}_{\rm
int})\right]\left(1+\delta\bar{\zeta}_{\rm int}^{2}\right)={\rm i},$ (6.10)

giving birth to a condensate. Indeed, installing a condensate between the two
such intersection points does not alter the the quasi-momentum and hence the
filling fraction stays intact. We have confirmed this to be the case by
numerically solving the Bethe equations for moderately large system sizes, as
demonstrated in Fig. 9 (again for the isotropic case). In particular, at low
filling fractions for both cuts their mutual “attraction” becomes apparent
(cf. the second cut connecting $(\zeta_{1},\bar{\zeta}_{1})$ in panel (a) in
Fig. 9). The four branch points in panel (a) in Fig. 9, reading
$\zeta_{1}=0.10884679+0.047665716{\rm i}$ and
$\zeta_{2}=0.07330641+0.04152184{\rm i}$ (with filling fractions
$\nu_{1}=0.02$, $\nu_{2}=0.06$, $\ell=1$ and mode numbers $n_{1}=1$,
$n_{2}=2$), have been determined numerically using the recipe given in
Appendix E.2. At a certain critical value of the filling fractions the two
cuts merge together. The intersection point becomes a logarithmic branch point
of a condensate, as exemplified in panel (b) in Fig. 9. The four branch points
in panel (b) in Fig. 9 are $\zeta_{1}=0.09587725+0.05961115{\rm i}$ and
$\zeta_{2}=0.07169871+0.04814544{\rm i}$ with filling fractions
$\nu_{1}=0.025$, $\nu_{2}=0.075$, $\ell=1$ and mode numbers $n_{1}=1$,
$n_{2}=2$. For the anisotropic interaction we encountered the same numerical
difficulties as previously for the one-cut solution with $n\geq 2$, cf.
Appendix E.1. We nevertheless do not expect any qualitative difference
compared to the isotropic model.

###  6.3 Special case: bion

As discussed earlier, the easy-axis regime (i.e. for $\delta>0$) permits for a
distinguished subclass of two-cut solutions that do not take place in the
other two (that is isotropic and easy-plane) regimes. Here we have in mind the
bion solution which we have described and parametrised in Sec. 5.2. One part
of the motivation for investigating this case in detail is to elucidate the
microscopic origin and stability of kinks in Landau–Lifshitz field theory,
which we expect to have a pivotal importance for understanding the freezing
property of a domain-wall profile, investigated recently in Refs. [16, 17].

Figure 10: Quantised bion configuration with the proposed density contours
(physical cuts). The two nonlinear modes with $\delta=1$, parametrised by
pairs of branch points $(\zeta_{1},\bar{\zeta}_{1})$ and
$(\zeta_{2},\bar{\zeta}_{2})$ on imaginary $\zeta$-axis, have associated mode
numbers $n_{1}=1$ and $n_{2}=-1$. The branch points have been found
numerically and are located at $\zeta_{1}=(1/0.9){\rm i}$ and
$\zeta_{2}=(1/1.058355921){\rm i}$, The red dashed line represents the “double
condensate”. The corresponding classical solution is plotted in Fig. 4.

To quantise the classical bion solution we demand the same reality condition
as previously in the one-cut case. Notice however that bion solutions belong
to maximally saturated states with the total filling $\nu=\frac{1}{2}$.
Condensates appear to be a common feature at half filling. In describing a
bion solution, there is no loss of generality in fixing the mode numbers of
the two cuts to $\pm 1$, such that $\Delta n=2$. Recall that in a general
situation with two cuts being far apart, each cut can grow a condensate on its
own. The bion case is different in that the two branch cuts reside close to
each other and share a condensate in common. In fact, in the isotropic
ferromagnet studied in Refs. [13, 56] the two-cut solution with mode numbers
set to $\pm 1$ is known to result in a “double condensate”. Led by this
observation, we conjecture that the same phenomenon takes place presently in
the case of bions; a “double condensate” emerges when a pair of branch cuts
with mode numbers $\pm 1$ intersect with one another, thereby producing a
logarithmic cut with of ‘doubled’ uniform density $2{\rm i}$. Analogous
objects which are twice as dense as ordinary Bethe strings have been
previously found in Ref. [56] in the study of solutions to the isotropic Bethe
equations.

We proceed by semi-classically quantising the bion solution using the
conjectured form of its contours with a double condensate, depicted in Fig. 10
by the red dashed line satisfying

$\rho(\zeta)(1+\delta\zeta^{2})=2{\rm i}.$ (6.11)

We have been able to explicitly match the classical values of the filling
fraction, momentum and energy: the two partial filling fractions add up
exactly to one half, i.e. $\nu=\nu_{1}+\nu_{2}=\frac{1}{2}$, whereas total
momentum $P=0\,(\mathrm{mod}\,2\pi)$ and total energy $E=3.96045$ (obtained by
numerically integrating along the proposed contours) match those of a
classical bion configuration, with $P=0$ and $E=3.960358(6)$. These results
strongly indicate that we have indeed correctly identified the physical
contours associated to a quantised bion.

As discussed earlier in Section 5.2, kinks arise as a particular (soliton)
limit of a bion solution in which the two branch points $\zeta_{1},\zeta_{2}$
on the imaginary axis coalesce at ${\rm i}/\epsilon$. By inspecting this
degeneration process at the level of semi-classical eigenstates, we find a
uniform condensate with a double density of Bethe roots running along the
imaginary axis between $-{\rm i}/\epsilon$ and ${\rm i}/\epsilon$. We note
that (anti)kinks are not compatible with periodic boundary conditions. In an
infinitely extended quantum chain however, the kink and antikink eigenstates
represent extra degenerate ground states (with broken translational symmetry)
of the XXZ ferromagnet in the gapped phase. Kink eigenstates have been derived
in Refs. [66, 67] using a curious correspondence between the XXZ quantum chain
and the problem of a melting crystal corner. This derivation however does not
require any use of the Bethe quantisation condition and consequently cannot
reveal the internal magnon structure of the kink. It would be valuable to
devise a method to extract the corresponding numerical solution to the Bethe
equations for large system sizes. The task of solving the anisotropic Bethe
equations (3.3) in the vicinity of half filling remains quite challenging at
this moment. Perhaps one could get some hints by first scanning through the
complete list of exact eigenstates for relatively small system sizes
(typically of order $L\sim 10$, using e.g. the techniques proposed in Refs.
[68, 69]) and look for traces of finite-size bions. At this junction, our
statements regarding the kink solution thus remain to an extent conjectural.

##  7 Semi-classical norms and overlaps

In this section we outline how to compute an overlap between two semi-
classical Bethe eigenstates. We provide closed formulae for (i) the Gaudin
norm and (ii) the Slavnov inner product between an on-shell and off-shell
Bethe states [70]. There are two possible routes to achieve this. The first
one, proposed by Gromov et al. in [31], is to perform coarse-graining directly
at the level of the general determinant formula for a specific finite-gap
density resolvent. The other approach, developed for the isotropic (XXX)
Heisenberg model by Kostov and Bettelheim in Refs. [32, 33, 34], makes use of
functional integration techniques with a bit of complex analysis. Both methods
provide the leading (i.e. classical) contributions to the overlaps and norms.
We shall not repeat the derivations here but instead only succinctly summarise
the main formulae for the model of our interest. Moreover, in Sections D.1 and
D.2 we provide a direct numerical confirmation based on the finite-size
analysis.

###  7.1 Gaudin norm

The method proposed by Kostov in Refs. [32, 33] has already been generalised
for the specific case of the anisotropic Heisenberg model in [71, 72]. To
compute the Gaudin norm we instead employ a more direct approach of [31],
which we generalise here by including the interaction anisotropy. The idea is
to convert the logarithm of the Gaudin determinant into a Riemann summation
which, after taking the $L\to\infty$ limit, corresponds to complex integration
along the physical contours which support the Bethe roots.

The Gaudin norm of a finite-volume Bethe eigenstate,

$\mathcal{N}=\langle\\{\vartheta\\}|\\{\vartheta\\}\rangle,$ (7.1)

grows exponentially in system size, i.e. $\log\mathcal{N}\sim\mathcal{O}(L)$
to the leading order of $L$. In the $\zeta$-plane parametrisation, we find the
following explicit form

$\begin{split}\log\mathcal{N}&=C_{1}L+o(L),\\\ C_{1}&=\int_{\mathcal{C}}{\rm
d}\zeta\left[{\rm
i}\pi\ell(1+\delta\,\zeta^{2})\rho(\zeta)+2\int_{0}^{\rho(\zeta)(1+\delta\,\zeta^{2})}\frac{{\rm
d}\xi}{1+\delta\,\xi^{2}}\log\big{(}(2\sinh(\pi\xi)\big{)}\right],\end{split}$
(7.2)

where we have expressed the volume-law coefficient $C_{1}$ in terms of the
resolvent density $\rho(\zeta)$ with support on a union of contours
$\mathcal{C}=\cup_{j}\mathcal{C}_{j}$.

We note that the dominant subleading correction to the above formula is quite
subtle and has the form $\log\mathcal{N}(L)=C_{1}\,L+\frac{1}{2}\log
L+\mathcal{O}(L^{0})$ (for $\quad C_{1}\,\in\mathcal{O}(1)$), as discussed in
Ref. [31]. Several numerical verifications (both with or without a condensate)
are presented in Appendix D.1.

###  7.2 Slavnov overlap

To express the semi-classical overlaps we follow instead the functional
integral approach devised in Refs. [32, 33]. This method does not rely on the
clustering properties of Gaudin determinant as in Ref. [31]. Here we merely
quote the final result of [71] (cf. formula (1.5) in there) for the
anisotropic Landau–Lifshitz field theory

$\log\langle\\{\vartheta_{1}\\}|\\{\vartheta_{2}\\}\rangle\simeq\oint_{\mathcal{C}_{1\cup
2}}\frac{{\rm d}\zeta}{2\pi{\rm
i}}\log\Phi_{\sqrt{\eta}}\big{(}\mathfrak{p}_{1}(\zeta)+\mathfrak{p}_{2}(\zeta)+\pi\big{)},$
(7.3)

involving two classical quasimomenta $\mathfrak{p}_{1}(\zeta)$ and
$\mathfrak{p}_{2}(\zeta)$ that correspond to the semi-classical Bethe
eigenstates $\mathinner{|{\\{\vartheta_{1}\\}}\rangle}$ and
$\mathinner{|{\\{\vartheta_{2}\\}}\rangle}$ 101010Only one of the states has
to be on-shell, i.e. solution to Bethe ansatz equations (3.3)., respectively.
Function $\Phi_{b}(z)$ stands for quantum dilogarithm [73], defined through
the following integral representation

$\Phi_{b}(z)=\exp\left[\frac{{\rm i}}{2}\int_{\mathbb{R}+{\rm i}0}\frac{{\rm
d}t}{t}\frac{e^{zt}}{\sin(b^{2}t)\sinh(\pi t)}\right],$ (7.4)

This function can be understood of as a ‘quantum deformation’ 111111The word
‘quantum’ here refers to the $q$-deformation parameter of ‘quantum calculus’
which should not be confused with the $q$-parameter of the quantum algebra
$\mathcal{U}_{q}(\mathfrak{su}(2))$ of the underlying anisotropic Heisenberg
chain. of the ordinary dilogarithm function $\mathrm{Li}_{2}(z)$ to which it
reduces in the isotropic limit $\delta\to 0$. The contour prescription in Eq.
(7.3) is such that $\mathcal{C}_{1,2}$ wrap around tightly around the supports
of the respective density resolvents, cf. Ref. [33].

Further simplification of the above formula can be made in the semi-classical
limit $\eta\to\epsilon/L$ which implies $b=\sqrt{\eta}\to 0$. In this limit
the quantum dilogarithm simplifies to [74] 121212Beware that the definition of
$\Phi_{b}(z)$ in [74] differs from the definition in [71] and the one used
here.

$\lim_{b\to 0}\Phi_{b}(z)=\exp\left[\frac{{\rm
i}L}{\sqrt{\delta}}\mathrm{Li}_{2}(-e^{{\rm i}z})\right]+\mathcal{O}(L^{0}),$
(7.5)

and accordingly at the leading order $\mathcal{O}(L)$ the logarithmic overlap
is approximately

$\begin{split}&\lim_{\eta\to\epsilon/L}\log\langle\\{\vartheta\\}_{1}|\\{\vartheta\\}_{2}\rangle=C_{2}L+o(L),\\\
&C_{2}=\frac{1}{\sqrt{\delta}}\oint_{\mathcal{C}_{1\cup 2}}\frac{{\rm
d}\zeta}{2\pi(1+\delta\,\zeta^{2})}\mathrm{Li}_{2}\left[e^{{\rm
i}\big{(}\mathfrak{p}_{1}(\zeta)+\mathfrak{p}_{2}(\zeta)\big{)}}\right].\end{split}$
(7.6)

Similarly to the Gaudin norm, the dominant subleading correction is of form
$\log\langle\\{\vartheta\\}_{1}|\\{\vartheta\\}_{2}\rangle=C_{2}L+\frac{1}{2}\log
L+\mathcal{O}(L^{0})$.

For the coinciding sets of rapidities, this correctly reproduces the leading
order expression for the Gaudin norm,

$\lim_{\eta\to\epsilon/L}\log\langle\\{\vartheta\\}|\\{\vartheta\\}\rangle\simeq\frac{L}{\sqrt{\delta}}\oint_{\mathcal{C}}\frac{{\rm
d}\zeta}{\pi(1+\delta\,\zeta^{2})}\mathrm{Li}_{2}\left(e^{2{\rm
i}\mathfrak{p}(\zeta)}\right).$ (7.7)

which can be readily reconciled with Eq. (7.2) upon expressing the resolvent
density in terms of the quasi-momentum, ${\rm i}\pi\rho(\zeta)=\pi
n-\mathfrak{p}(\zeta)$, and performing the following integral

$\int_{0}^{\rho(\zeta)(1+\delta\zeta^{2})}\frac{{\rm
d}\xi}{1+\delta\xi^{2}}\log\big{(}2\sinh(\pi\xi)\big{)}=\frac{1}{2\pi}\mathrm{Li}_{2}\left(e^{2{\rm
i}\mathfrak{p}(\zeta)}\right)+\frac{\pi}{2}\rho^{2}(\zeta)(1+\delta\zeta^{2})^{2}-\frac{\pi}{12}.$
(7.8)

There is a practical limitation of Eq. (7.6) that concerns the placement of
integration contours, acknowledged previously in Ref. [33]. The requirement is
that the integration contours $\mathcal{C}$ must avoid crossing any branch cut
of the function in the integrand. This issue is presently further complicated
by additional logarithmic branch cuts due to the dilogarithm function. This
shortcoming makes the numerical verification a challenging task. We
nonetheless still managed to verify its validity in special case of overlaps
with a vacuum descendant (see Appendix D.2).

##  8 Correlation functions

We have thus far demonstrated that the knowledge of physical contours not only
proves useful in calculating the spectrum, overlaps and norms of semi-
classical Bethe eigenstates, but also facilitates the semi-classical
quantisation of the finite-gaps solutions. On the other hand, we have not yet
addressed the expectation values of physical observables. This section is
devoted to discuss some properties of correlation functions at the semi-
classical level.

Despite integrability, the task of computing exact expectation values of
physical observables, including their correlation functions, appears quite
challenging. There have already been numerous works on this subject, employing
either the form-factor expansion or bootstrap methods, e.g. [38, 75, 76, 77,
78]. Here, however, we are particularly interested in the semi-classical
regime where those methods are not directly applicable.

We are specifically interested whether the aforementioned classical–quantum
correspondence for correlations functions, established analytically in the
introductory section 2 on the toy example of the harmonic oscillator, holds on
more general grounds. A direct generalisation of this principle from a single-
particle paradigm is complicated by the fact that integrable field theories
governed by PDEs involve many degrees of freedom which, moreover, couple (i.e.
interact) in a non-trivial fashion. One viable empirical approach to obtain
correlation functions in classical regime (enabling a comparison with their
quantum counterparts) is to build on the semi-classical form-factor approach
developed by Smirnov [79]. Accordingly, the matrix elements would be
represented as integrals over $\gamma$ variables, see Eqs. (33)–(36) in [79].

Let us consider a periodic solution associated with one branch cut. It is
described by a single dynamic variable $\gamma$. We further replace the
integration over variable $\gamma$ with the integration over one period,
similar to the phase space averaging [80]. In addition, we compute numerically
(for small system sizes) the quantum correlators in the eigenstate that
corresponds to the cut.

(a)

(b)

Figure 11: (a) Quantum correlation functions $\langle\hat{\sigma}^{\rm
x}_{j}\hat{\sigma}^{\rm x}_{j+n}\rangle$ (shown in green circles for $L=15$
and black crosses for $L=18$) versus averaged classical correlation functions
(red dashed line), with anisotropy $\delta=1$ and classical period $\ell=1$.
(b) The same comparison for longitudinal correlators $\langle\hat{\sigma}^{\rm
z}_{j}\hat{\sigma}^{\rm z}_{j+n}\rangle$.

We focus only on static (i.e. equal-time) correlation functions, considering
one-point functions $\langle\hat{\sigma}^{\rm x}_{j}\rangle$ and
$\langle\hat{\sigma}^{\rm z}_{j}\rangle$ and two-point point functions
$\langle\hat{\sigma}^{\rm x}_{j}\hat{\sigma}^{\rm x}_{k}\rangle$ and
$\langle\hat{\sigma}^{\rm z}_{j}\hat{\sigma}^{\rm z}_{k}\rangle$.
Specifically, we set the filling fraction to $\nu=\frac{1}{3}$ and the mode
number to $n=1$. We present our results for system sizes $L=15$ and $18$ with
the with number of down-turned spins being $5$ and $6$, respectively. The
corresponding averages are denoted as $\langle\cdots\rangle_{5}$ and
$\langle\cdots\rangle_{6}$. The coordinate Bethe ansatz wavefunctions can be
represented in the local eigenstate basis of $\hat{\sigma}^{\rm z}$ by solving
Bethe ansatz equations (3.3). An important thing to keep in mind is that
quantum states (wavefunctions) are eigenfunctions of the momentum operator and
consequently expectations values of one-point observables do not have any
dependence on the spatial coordinate (lattice index). Two-point functions on
the other hand can only be a function of the distance. In contrast, classical
spin field configurations exhibit non-uniform dependence on the spatial
coordinate $x$ and can be thus compared to the quantum correlation functions
evaluated on semi-classical eigenstates after an appropriate phase-space
averaging. In particular, for a periodic classical spin-field configuration
with period $\ell$ and an operator (product of operators) $O[\mathcal{S}]$
that functionally depends on the spin configuration $\mathcal{S}$, we expect
the corresponding correlation functions to take the form

$\langle\\{\vartheta\\}|O[\mathcal{S}]|\\{\vartheta\\}\rangle\simeq\left\langle
O[\mathcal{S}]\right\rangle_{\rm cl}=\frac{1}{\ell}\int_{0}^{\ell}{\rm
d}x\,O[\mathcal{S}(x)],$ (8.1)

with $|\\{\vartheta\\}\rangle$ denoting the corresponding semi-classical
quantum eigenstate in the thermodynamic limit $L\to\infty$. The classical spin
configuration for $n=1$ and $\nu=\frac{1}{3}$ reads

$\begin{split}&\mathcal{S}^{\rm z}(x,t)=1-2\nu=\frac{1}{3},\\\
&\mathcal{S}^{\rm
x}=4(\nu-\nu^{2})\cos(kx+wt)=\frac{8}{9}\cos\left[\frac{2\pi}{\ell}x+\frac{1}{3}\left(\frac{4\pi^{2}}{\ell^{2}}+\delta\right)t\right].\end{split}$
(8.2)

Therefore for the one-point functions, we find

$\langle\hat{\sigma}^{\rm x}_{j}\rangle_{5}=\langle\hat{\sigma}^{\rm
x}_{k}\rangle_{6}=0=\frac{1}{\ell}\int_{0}^{\ell}{\rm d}x\,\mathcal{S}^{\rm
x}(x),\qquad\langle\hat{\sigma}^{\rm
z}_{j}\rangle_{5}=\langle\hat{\sigma}^{\rm
z}_{k}\rangle_{6}=\frac{1}{\ell}\int_{0}^{\ell}{\rm d}x\,\mathcal{S}^{\rm
z}(x)=\frac{1}{3},$ (8.3)

for $j=1,2,\ldots 15$ and $k=1,2,\ldots,18$, thus confirming the
correspondence.

Our results for two point functions are presented in Fig. 11. In this case we
notice some discrepancies between the classical expectation values,

$\langle\hat{\sigma}^{\rm x}_{j}\hat{\sigma}^{\rm
x}_{j+n}\rangle\simeq\frac{1}{\ell}\int_{0}^{\ell}{\rm d}x\,\mathcal{S}^{\rm
x}(x)\mathcal{S}^{\rm
x}\left(x+\frac{n\ell}{L}\right)=2(\nu-\nu^{2})\cos\left(\frac{2\pi
n}{L}\right),$ (8.4) $\langle\hat{\sigma}^{\rm z}_{j}\hat{\sigma}^{\rm
z}_{j+n}\rangle\simeq\frac{1}{\ell}\int_{0}^{\ell}{\rm d}x\,\mathcal{S}^{\rm
z}(x)\mathcal{S}^{\rm z}\left(x+\frac{n\ell}{L}\right)=(1-2\nu)^{2},$ (8.5)

and their quantum counterparts, which we attribute to the finite number of
spins. Indeed, instead of having a condensed contour of Bethe roots
representing the branch cut in the complex plane, we consider solutions with
at most 6 Bethe roots. We nonetheless find it plausible that with increasing
system sizes the deviations would gradually diminish and we therefore expect
to recover the classical result in the $L\to\infty$ limit. Finite-size effects
also depend on type of operators that appear in the correlator; in the case of
$\langle\hat{\sigma}^{\rm z}_{j}\hat{\sigma}^{\rm z}_{j+n}\rangle$ the
deviation from the asymptotic result is found to be larger than in the case of
$\langle\hat{\sigma}^{\rm x}_{j}\hat{\sigma}^{\rm x}_{j+n}\rangle$, see Fig.
11. Here we have demonstrated the first step in understanding such correlation
functions. A systematic and comprehensive numerical analysis of the
correspondence and finite-size corrections is postponed for future work.

##  9 Conclusion and outlook

We have studied the structure of the semi-classical spectrum of the