Scaling limit of the ground state Bethe roots for the inhomogeneous XXZ spin -
$\frac{1}{2}$ chain

Sascha Gehrmann1, Gleb A. Kotousov1 and Sergei L. Lukyanov2

1Institut f$\ddot{{\rm u}}$r Theoretische Physik, Leibniz Universit$\ddot{{\rm
a}}$t Hannover

Appelstraße 2, 30167 Hannover, Germany

2NHETC, Department of Physics and Astronomy

Rutgers University

Piscataway, NJ 08855-0849, USA

Abstract

It is known that for the Heisenberg XXZ spin - $\frac{1}{2}$ chain in the
critical regime, the scaling limit of the vacuum Bethe roots yields an
infinite set of numbers that coincide with the energy spectrum of the quantum
mechanical 3D anharmonic oscillator. The discovery of this curious relation,
among others, gave rise to the approach referred to as the ODE/IQFT (or
ODE/IM) correspondence. Here we consider a multiparametric generalization of
the Heisenberg spin chain, which is associated with the inhomogeneous six-
vertex model. When quasi-periodic boundary conditions are imposed the lattice
system may be explored within the Bethe Ansatz technique. We argue that for
the critical spin chain, with a properly formulated scaling limit, the scaled
Bethe roots for the ground state are described by second order differential
equations, which are multi-parametric generalizations of the Schrödinger
equation for the anharmonic oscillator.

###### Contents

  1. 1 Introduction and summary
  2. 2 The case odd $r$ and $A=\frac{r-1}{2}$
    1. 2.1 Scaling limit of the Bethe roots for the free fermion point
    2. 2.2 Differential equation
    3. 2.3 Relation between the RG invariants and the coefficients of the ODE
    4. 2.4 Some remarks about the scaling limit
    5. 2.5 Low energy spectrum in the scaling limit
      1. 2.5.1 Low energy spectrum for the ${\cal Z}_{r}$ invariant spin chain
      2. 2.5.2 Some comments on the low energy spectrum away from the ${\cal Z}_{r}$ invariant point
  3. 3 Basic constraints on ODEs for spin chain with anisotropy $0<\gamma<\pi$
    1. 3.1 The case $A=r-1$
    2. 3.2 General restrictions on integers $\mu$ and $j$
  4. 4 ODE with single coefficient $c_{\mu,j}$ for $A\mu>rj$ and $A=1,\ldots,r-2$
    1. 4.1 Bohr-Sommerfeld quantization condition
    2. 4.2 Specification of the scaling limit
      1. 4.2.1 Example: $r=5$, $A=1$ and $(\mu,j)=(2,0)$
      2. 4.2.2 Example: $r=5$, $A=1$ and $(\mu,j)=(1,0)$
      3. 4.2.3 Example: combining the cases $\mu=1$ and $\mu=2$ for $r=5$, $A=1$, $j=0$
      4. 4.2.4 Example: $r=7$, $A=1$ and $(\mu,j)=(1,0)$
      5. 4.2.5 The case $\frac{rj}{A}<\mu<\frac{r}{A+1}\,(j+\frac{1}{2})$
  5. 5 The case $A\mu=rj$ with $A=1,\ldots,r-2$
    1. 5.1 Special features of the scaling limit
    2. 5.2 Relation to $A=0$
  6. 6 General differential equation for $A=1,2,\ldots,r-2$
    1. 6.1 The case $A=1$
    2. 6.2 The case even $r$ and $A=\frac{r}{2}-1$
    3. 6.3 The case even $r$ and $A=\frac{r}{2}$
    4. 6.4 The case $A=r-2$
  7. 7  ODEs for $A=0$
  8. 8 Concluding remarks
  9. A Auxiliary functions appearing in the sum rules
  10. B Lagrange formula
  11. C Relations between $\mathfrak{a}_{2j+1}$ and $c_{2j+1}$ for $r=3,5,7$ and $A=\frac{r-1}{2}$
  12. D Scaling limit at the free fermion point with $\eta_{\ell}=(-1)^{r}\,\mbox{e}^{\frac{{\rm i}\pi}{r}\,(2\ell-1)}$
  13. E Scaling limit for $r$, $A$ even and $(\mu,j)=(\frac{r}{2},\frac{A}{2})$
  14. F Supplementary tables

##  1 Introduction and summary

Since the seminal work of Polyakov [1], it is believed that the isotropic and
homogeneous scaling behaviour of a critical statistical system can be
adequately described within the framework of Conformal Field Theory (CFT).
While this fundamental insight remains unproven in general, there exists an
enormous amount of supporting evidence. Much of it comes from the study of 2D
classical statistical models. In many cases, such models admit an equivalent
description in terms of 1D quantum spin chains. Then the problem of
identification of the underlying CFT can be approached by studying the scaling
limit of the spectrum of a spin chain Hamiltonian [2]. This provides a
convenient way for a numerical investigation of the critical behaviour. For
the so-called integrable statistical systems, a detailed analytical study is
also possible. By integrable, what is usually meant is that the model can be
explored via the Bethe Ansatz approach or some variant thereof [3].

In the case of integrable systems at criticality, apart from the conformal
structure of the underlying continuous theory, it is natural to investigate
also the integrable structure of the corresponding CFT [4, 5, 6]. The latter
is encoded in the scaling limit of the Bethe roots – the solutions of the
Bethe Ansatz equations. The most powerful approach for studying the integrable
structure of CFT or, more generally, QFT is the so-called ODE/IM
correspondence [7, 8, 9, 10]. The abbreviation IM stands for Integrable Model.
For the most interesting examples these are 2D Integrable Quantum Field
Theories (IQFT) and, as such, we prefer to use the term ODE/IQFT
correspondence rather than ODE/IM.

One of the simplest manifestations of the ODE/IQFT correspondence, which is
directly relevant to this work, occurs in the scaling limit of the XXZ spin -
$\frac{1}{2}$ chain. The latter is the archetypical example of an integrable
model and, in fact, the Bethe Ansatz approach was originally introduced for
its special (isotropic) case [11]. The Hamiltonian of the spin-$\frac{1}{2}$
XXZ chain of length $N$ is given by

$\displaystyle\mathbb{H}_{\rm
XXZ}=-\sum_{m=1}^{N}\Big{(}\sigma_{m}^{x}\sigma_{m+1}^{x}+\sigma_{m}^{y}\sigma_{m+1}^{y}+\Delta\,\sigma_{m}^{z}\sigma_{m+1}^{z}\Big{)}\
.$ (1.1)

With quasi-periodic boundary conditions imposed [14],

$\sigma^{x}_{N+m}\pm{\rm i}\,\sigma^{y}_{N+m}=\omega^{\pm
2}\,\big{(}\,\sigma^{x}_{m}\pm{\rm
i}\,\sigma^{y}_{m}\,\big{)}\,,\qquad\qquad\sigma_{N+m}^{z}=\sigma^{z}_{m}\,,$
(1.2)

$\mathbb{H}_{\rm XXZ}$ commutes with the $z$ \- component of the total spin
operator and the eigenstates are labelled by the (half-)integer numbers
$S^{z}$. The energy spectrum is described in terms of solutions of the Bethe
Ansatz equations:

$\bigg{(}\frac{1+q^{+1}\,\,\zeta_{j}}{1+q^{-1}\,\zeta_{j}}\bigg{)}^{N}=-\omega^{2}\,q^{2S^{z}}\,\prod_{i=1}^{N/2-S^{z}}\,\frac{\zeta_{i}-q^{+2}\,\zeta_{j}}{\zeta_{i}-q^{-2}\,\zeta_{j}}\
,$ (1.3)

where $q$ parameterizes the anisotropy as $\Delta=\frac{1}{2}(q+q^{-1})$. It
turns out that the system is critical when $-1\leq\Delta<1$ [12, 13, 14] or,
equivalently, $q$ is a unimodular number, i.e.,

$q=\mbox{e}^{{\rm i}\gamma}\ .$ (1.4)

In this case, for an even number of sites, the ground state of
$\mathbb{H}_{\rm XXZ}$, i.e., the state with the lowest energy, is non-
degenerate and occurs in the sector $S^{z}=0$. It was proven in the work of
Yang and Yang [15] that for $\gamma\in(0,\pi]$ and periodic boundary
conditions $(\omega=1)$, the solution set to the Bethe Ansatz equations
corresponding to the ground state has $\zeta_{i}$ all real, distinct positive
numbers. The proof may be extended to quasi-periodic boundary conditions when

$\omega=\mbox{e}^{{\rm i}\pi{\tt k}}$ (1.5)

with sufficiently small real ${\tt k}$, namely,

$|{\tt k}|<\frac{\gamma}{2\pi}\ .$ (1.6)

The vacuum Bethe roots can be naturally ordered as

$0<\zeta_{1}<\zeta_{2}<\ldots<\zeta_{N/2}\ .$ (1.7)

In the limit when the number of sites $N$ goes to infinity, most of the roots
become densely distributed such that $(\zeta_{j}-\zeta_{j+1})/\zeta_{j}\sim
1/N$. However, at the edges of the distribution the Bethe roots develop a a
certain scaling behaviour. In particular, keeping $j$ fixed as $N\to\infty$
the following limits exist

$\lim_{N\to\infty\atop j-{\rm fixed}}N^{2(1-\gamma/\pi)}\zeta_{j}$ (1.8)

and form an infinite set of non-vanishing numbers. These turn out to admit a
remarkable description in terms of the Schrödinger equation for the 3D
anharmonic oscillator [9],

$\displaystyle\bigg{(}-{\frac{{\mbox{d}}^{2}}{{\mbox{d}}x^{2}}}+\frac{l(l+1)}{x^{2}}+x^{2\alpha}-E\bigg{)}\,\Psi=0\
.$ (1.9)

Namely, the differential equation possesses only a discrete component of the
spectrum and the corresponding energy levels $E_{j}$, coincide with the limits
(1.8) up to an overall factor provided the parameters of the ODE are taken to
be

$\displaystyle\alpha=\frac{\pi}{\gamma}-1\ ,\ \ \ \ \ \ \ \ \ \ \ \
l+\frac{1}{2}=\frac{\pi{\tt k}}{\gamma}\ .$ (1.10)

A similar relation between the scaling behaviour of the Bethe roots of not
only the ground state, but also the class of so-called low energy excited
states, holds true (see, e.g., ref.[16]). Moreover such a phenomenon, i.e.,
the ODE/IQFT correspondence, has been observed for a number of other
integrable, critical statistical models as well [17, 18, 19, 20, 21].

The XXZ spin chain belongs to the integrability class of a 2D classical
statistical system known as the six-vertex model. As was shown by Baxter in
the work [22], the latter admits a multiparametric generalization, which is
solvable within the Bethe Ansatz approach. In this case, the equations (1.3)
are modified to

$\prod_{J=1}^{N}\frac{\eta_{J}+q^{+1}\,\,\zeta_{j}}{\eta_{J}+q^{-1}\,\zeta_{j}}=-\omega^{2}\,q^{2S^{z}}\,\prod_{i=1}^{N/2-S^{z}}\,\frac{\zeta_{i}-q^{+2}\,\zeta_{j}}{\zeta_{i}-q^{-2}\,\zeta_{j}}\,\qquad\qquad(\,i=1,2,\ldots,\tfrac{N}{2}-S^{z}\,)$
(1.11)

and involve the complex parameters $\\{\eta_{J}\\}_{J=1}^{N}$, the so-called
inhomogeneities. The corresponding statistical system is known as the
inhomogeneous six-vertex model. With the parameters $q$ and $\eta_{J}$ obeying
certain conditions, the lattice system develops critical behaviour. We’ll
assume that the number of sites $N$ is divisible by $r$ and the
inhomogeneities obey the $r$-site periodicity condition

$\eta_{J+r}=\eta_{J}\qquad\qquad\qquad(\,J=1,2,\ldots,N-r\,;\
N/r\in\mathbb{Z}\,)\,.$ (1.12)

In the scaling limit, $N\to\infty$ while the integer $r$ is kept fixed. This
ensures the presence of translational invariance for the continuous theory.

The critical behaviour of the $r$ \- site periodic inhomogeneous six-vertex
model can be explored by studying the spectrum of the Hamiltonian $\mathbb{H}$
of a certain spin chain. For the case $r=1$ it coincides, up to an overall
multiplicative factor and additive constant, with $\mathbb{H}_{\rm XXZ}$
(1.1). The explicit form of $\mathbb{H}$ with $r\geq 2$ is rather cumbersome
(see, e.g., [23, 24] and the work [25], where the same notation as in this
paper is used). What is important is that $\mathbb{H}$ is given by a sum over
terms that describe the interactions of up to $r+1$ adjacent $\frac{1}{2}$
spins, i.e., the corresponding Hamiltonian density is local. The quasi-
periodic boundary conditions (1.2) are still assumed and, similar to the
homogeneous case, the Hamiltonians commute with the $z$ \- projection of the
total spin operator.

A numerical investigation shows that the spin chain governed by the
Hamiltonian $\mathbb{H}$ exhibits critical behaviour when the anisotropy
parameter $q$ is a unimodular number, $q=\mbox{e}^{{\rm i}\gamma}$. It turns
out that the phase can be restricted to the domain $\gamma\in(0,\pi)$ and the
critical behaviour is described differently as
$\frac{\pi}{r}\,A<\gamma<\frac{\pi}{r}\,(A+1)$. Inside each segment labelled
by the integer $A$ we will parameterize $\gamma$ by real positive $n$ such
that

$\gamma=\frac{\pi}{r}\,A+\frac{\pi}{n+r}\
\qquad\qquad\qquad\qquad(A=0,1,\ldots,r-1\,,\;n>0)\ .$ (1.13)

Mapping out the critical surfaces in the space of the parameters
$\big{(}\gamma,\\{\eta_{\ell}\\}\big{)}$ is a complicated problem that, apart
from some special cases, has not been fully achieved. A useful starting point
is when the inhomogeneities are given by

$\eta_{\ell}=(-1)^{r}\,\mbox{e}^{\frac{{\rm i}\pi}{r}\,(2\ell-1)}\ \ \
\qquad\qquad(\ell=1,\ldots,r)$ (1.14)

for which the model possesses ${\cal Z}_{r}$ invariance (see, e.g., sec. 7 in
ref.[25]). In this case the description of the scaling limit for the ground
state is especially simple. For $N$ divisible by $2r$ and $S^{z}=0$ the
corresponding Bethe roots are split into $r$ groups of equal size. The Bethe
roots within each group lie on the rays

$\displaystyle\arg(\zeta)=-\frac{\pi}{r}\,A,\ \frac{\pi}{r}\,(2-A),\
\ldots,\,\frac{\pi}{r}\,\big{(}2\,(r-1)-A\big{)}\ \ \ \ \ \ \ \ \ \ \ ({\rm
mod}\,2\pi)\ .$ (1.15)

These can be labelled by the integer $a=1,2,\ldots,r$ and we denote by
$\zeta_{m}^{(a)}$ the roots lying on the $a^{\rm th}$ ray ordered by their
absolute value. Then the real positive numbers
$\big{|}\zeta_{m}^{(a)}\big{|}^{r}$ with $m=1,2,\ldots,N/(2r)$ are the same
for any $a$. They turn out to coincide with the solution to (1.3) with
$S^{z}=0$ and the substitutions $q\mapsto\mbox{e}^{\frac{{\rm i}\pi r}{n+r}}$,
$N\mapsto N/r$ which corresponds to the ground state of the XXZ spin chain. It
follows that the limits

$E_{m}^{(a)}=\lim_{N\to\infty\atop m-{\rm
fixed}}\,\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{\frac{2n}{r(n+r)}}\zeta_{m}^{(a)}$
(1.16)

exist and are non-vanishing. Here for later convenience, an extra $n$
dependent factor has been introduced with

$N_{0}=\frac{\sqrt{\pi}\Gamma(1+\frac{r}{2n})}{r\Gamma(\frac{3}{2}+\frac{r}{2n})}\
.$ (1.17)

The numbers $E_{m}^{(a)}$ would be expressed in terms of the spectrum of the
ODE (1.9). For the purpose of this work, it is convenient to re-write the
differential equation using the variables
$y=\frac{2}{r}\log(x)+\frac{2}{n+r}\log(\frac{r}{2})$ and
$\psi=x^{-\frac{1}{2}}\,\Psi$, so that it becomes111With some abuse of
notation we use the same symbol $E$ in eqs. (1.9) and (1.18). In fact, $E$
from (1.9) coincides with the combination
$(-1)^{A}\,(2/r)^{\frac{2n}{n+r}}E^{r}$ with $E$ from (1.18). Also,
$\alpha=\frac{n}{r}$ and $l+\frac{1}{2}=\frac{2p}{r}$.

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}\bigg{]}\,\psi=0$ (1.18)

with

$\displaystyle p=\frac{n+r}{2}\ {\tt k}\ .$ (1.19)

The above differential equation admits the Jost solution $\psi_{p}(y)$, which
is uniquely defined through the asymptotic

$\psi_{p}(y)=\mbox{e}^{py}\qquad\qquad{\rm as}\qquad\qquad
y\to-\infty\qquad\qquad\qquad\qquad\big{(}\Re e(p)\geq 0\big{)}\ .$ (1.20)

For generic $E$ the solution grows unboundedly as $y\to+\infty$. However, at
the special values $E=E_{m}^{(a)}$ (1.16), one has

$\psi_{p}(y)\big{|}_{E=E_{m}^{(a)}}\to 0\qquad{\rm as}\qquad y\to+\infty\ .$
(1.21)

In this work we consider the spin chain, where the ${\cal Z}_{r}$ symmetry is
“softly” broken for finite system size. This is carried out by treating the
inhomogeneities $\eta_{\ell}$ as $N$ \- dependent bare couplings of the spin
chain Hamiltonian which tend to the values (1.14) as $N\to\infty$. The
quantities

$\displaystyle{\mathfrak{a}}_{s}=\frac{1}{s}\
\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{d_{s}}\ \frac{1}{r}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-s}$ (1.22)

will be taken to be independent of $N$, i.e., they play the rôle of the RG
invariants. With their proper specialization (which, among other things,
involves fixing the exponent $d_{s}>0$) the Bethe roots possess a scaling
behaviour similar to that as for the ${\cal Z}_{r}$ invariant case. In
particular, for the ground state in the sector with $S^{z}=0$, the roots may
still be split into $r$ groups with roughly equal phases. Then there exist the
limits (1.16), where the integer $m$ labels the roots according to their
absolute values. We will argue that $E_{m}^{(a)}$ can be described as above,
but with the differential equation (1.18) replaced by the ODE of the form

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-\sum_{(\mu,j)\in\,\bm{\Xi}_{r,A}}c_{\mu,j}\,E^{\mu}\
\mbox{e}^{\big{(}(A\mu-rj)\,\frac{n+r}{r}+\mu\big{)}y}\,\bigg{]}\,\psi=0\ .$
(1.23)

The coefficients $c_{\mu,j}$ are related to the values of the RG invariants
${\mathfrak{a}}_{s}$ from eq. (1.22), while the summation indices $\mu$ and
$j$ run over a certain integer set. The latter is given by222The number of
admissible pairs $(\mu,j)$ (1.24) coincides with the total number of solutions
of the Diophantine equations $A\mu-rj=L$ subject to the constraints
$1\leq\mu<r-L\,\&\ j\geq 0$ for $L=1,2,\ldots$ .

$\displaystyle\bm{\Xi}_{r,A}$ $\displaystyle=$
$\displaystyle\big{\\{}(\mu,j)\,:\ \ \tfrac{rj}{A}<\mu<\tfrac{r}{A+1}\ (j+1)\
\ \&\ \ j\geq 0\big{\\}}\qquad{\rm for}\qquad A=1,2,\ldots,r-2\,,$ (1.24)

while for $A=0$,

$\displaystyle\bm{\Xi}_{r,0}$ $\displaystyle=$
$\displaystyle\big{\\{}(\mu,j)\,:\ \ \mu=1,2,\ldots,r-1\ \ \&\ \
j=0\big{\\}}\,.$ (1.25)

As $0<A<r-1$ we expect that it is possible to organize the scaling limit in
such a way that the differential equation (1.23) appears for any given values
of the coefficients $c_{\mu,j}$ (some restrictions on the domain of $n$ may be
required). In particular, the specification of the RG invariants such that
there is only one non-vanishing term occurring in the sum in (1.23) can be
found in Appendix F for $r=3,4,\ldots,10$. When $A=0$ only the coefficient
$c_{\mu,0}$ with $\mu\geq\frac{r}{2}$ can be chosen at will, while the others
turn out to be not independent. The case $A=r-1$ was already discussed in
ref.[21]. Its peculiar feature is that all of the exponents $d_{s}$ in (1.22)
turn out to be zero. The scaling of the Bethe roots for the ground state is
still described by means of the differential equation of the form (1.23) with

$\displaystyle\bm{\Xi}_{r,r-1}=\big{\\{}(\mu,j)\,:\ \ \mu=j+1\ \ \&\ \
j=0,1,\ldots,r-2\big{\\}}\qquad\qquad\qquad(A=r-1)\,.$ (1.26)

##  2 The case odd $r$ and $A=\frac{r-1}{2}$

An important source of intuition in the study of the spin - $\frac{1}{2}$
Heisenberg chain comes from the XX spin chain — the special case when the
anisotropy parameter $\Delta$ in (1.1) is set to zero or, equivalently,
$q={\rm i}$. By means of the Jordan-Wigner transformation, the lattice model
can then be mapped to a system of free fermions, which dramatically simplifies
the analysis. A similar phenomenon occurs for the spin chain associated with
the $r$-site periodic inhomogeneous six-vertex model for any odd $r$,

$\displaystyle r=3,5,7,\ldots\,,$ (2.1)

when the integer $A$ in formula (1.13) for the argument of $q$ is chosen to be

$\displaystyle A=\frac{r-1}{2}\ .$ (2.2)

In this case

$\displaystyle q={\rm i}\,\mbox{e}^{\frac{{\rm i}\pi(r-n)}{2r(n+r)}}$ (2.3)

so that as $n=r$ one has $q={\rm i}$. Below we discuss the scaling limit of
the ground state for the spin chain in the regime (2.2) ($r$ odd) using the
“free fermion” point as a launch pad.

###  2.1 Scaling limit of the Bethe roots for the free fermion point

At $q={\rm i}$ the Bethe Ansatz equations (1.11) for the $r$-site periodic
model (1.12) become

$\bigg{(}\prod_{\ell=1}^{r}\frac{\eta_{\ell}+{\rm
i}\zeta_{j}}{\eta_{\ell}-{\rm
i}\zeta_{j}}\bigg{)}^{\frac{N}{r}}=-\mbox{e}^{{\rm i}\pi(2{\tt k}+S^{z})}\ .$
(2.4)

Taking the logarithm of both sides, one obtains

$\frac{N}{{\rm i}\pi r}\ \log\Bigg{(}\prod_{\ell=1}^{r}\frac{1+{\rm
i}\zeta_{j}/\eta_{\ell}}{1-{\rm
i}\zeta_{j}/\eta_{\ell}}\Bigg{)}=2n_{j}-1+2{\tt k}+S^{z}$ (2.5)

with $n_{j}$ being certain integers. For the vacuum state of the ${\cal
Z}_{r}$ invariant case, where the inhomogeneities are given by (1.14), the
Bethe roots are split into the groups $\zeta_{m}^{(a)}$ with $a=1,\ldots,r$.
Then it turns out that with a proper choice of the branch of the logarithm,

$\frac{N}{{\rm i}\pi r}\ \log\Bigg{(}\prod_{\ell=1}^{r}\frac{1+{\rm
i}\zeta_{m}^{(a)}/\eta_{\ell}}{1-{\rm
i}\zeta_{m}^{(a)}/\eta_{\ell}}\Bigg{)}=2m-1+2{\tt
k}\qquad\qquad\big{(}m=1,2,\ldots,N/(2r)\big{)}\ ,$ (2.6)

where $S^{z}$ has been set to zero. One may expect that the above formula
remains true when the inhomogeneities $\eta_{\ell}$ differ slightly from the
values $(-1)^{r}\,\mbox{e}^{\frac{{\rm i}\pi}{r}\,(2\ell-1)}$.

To perform the scaling limit, we write $\zeta$ in the form

$\displaystyle\zeta=\bigg{(}\frac{\pi}{2N}\bigg{)}^{\frac{1}{r}}\ E$ (2.7)

and keep $E$ fixed as $N\to\infty$. Further, the inhomogeneities are assigned
an $N$ \- dependence in such a way that the quantities

$\displaystyle{\mathfrak{a}}_{2j+1}$ $\displaystyle=$
$\displaystyle\frac{1}{2j+1}\
\bigg{(}\frac{2N}{\pi}\bigg{)}^{1-\frac{2j+1}{r}}\ \frac{1}{r}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-2j-1}\ \ \ \ \
\qquad\big{(}j=1,\ldots,\tfrac{r-1}{2}\big{)}$ (2.8)

remain independent of the length of the spin chain. For $j=A=\frac{r-1}{2}$,
the exponent of the $N$ \- dependent factor in the r.h.s. vanishes and,
without loss of generality, ${\mathfrak{a}}_{r}$ can be set to be any positive
number. We take its value to be the same as for the ${\cal Z}_{r}$ invariant
case (1.14), i.e., ${\mathfrak{a}}_{r}=\frac{1}{r}$. This way, in the scaling
limit, the Bethe Ansatz equations (2.6) simplify to

$\displaystyle\frac{(-1)^{A}}{r}\,E^{r}+\sum_{j=0}^{A-1}(-1)^{j}\,{\mathfrak{a}}_{2j+1}\
E^{2j+1}=2m-1+2{\tt k}\ \ \ \qquad\ \ \ \ (m=1,2,\ldots)\ .$ (2.9)

For given $m$, the above algebraic equation possesses exactly $r$ complex
solutions, which we denote as $E_{m}^{(a)}$. Assuming that the absolute value
of ${\mathfrak{a}}_{2j+1}$ is sufficiently small, the roots may be specified
by the condition that

$E_{m}^{(a)}\sim\mbox{e}^{\frac{{\rm i}\pi}{r}(2a-A)}\,\big{(}(2m-1+2{\tt
k})\,r\big{)}^{\frac{1}{r}}\qquad\qquad{\rm for}\qquad\qquad m\gg 1$ (2.10)

and $E_{m}^{(a)}$ are split into $r$ groups of roots which have approximately
the same phase

$\arg\big{(}E_{m}^{(a)}\big{)}\approx\frac{\pi}{r}\ \big{(}2a-A\big{)}\ \ \ \
\ \ \qquad(a=1,\ldots,r)\ .$ (2.11)

As such, the scaling limit of the Bethe roots for the ground state is
described by the formula

$E_{m}^{(a)}=\lim_{N\to\infty\atop m-{\rm
fixed}}\,\bigg{(}\frac{2N}{\pi}\bigg{)}^{\frac{1}{r}}\zeta_{m}^{(a)}\,.$
(2.12)

This provides an illustration of the scaling relation (1.16) discussed in the
Introduction.

It is clear how to come up with a differential equation, whose spectrum is
given by $E_{m}^{(a)}$. One just takes (1.18) with $n=r$ and replaces
$(-1)^{A}E^{r}$ by $r\lambda(E)$, where

$\lambda(E)=\frac{(-1)^{A}}{r}E^{r}+\sum_{j=0}^{\frac{r-3}{2}}(-1)^{j}\,{\mathfrak{a}}_{2j+1}\
E^{2j+1}\ .$ (2.13)

The ODE becomes

$\displaystyle\Big{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{2ry}-r\lambda(E)\
\mbox{e}^{ry}\Big{]}\ \psi=0\ ,$ (2.14)

which is a form of the confluent hypergeometric equation. The solution defined
via the asymptotic condition (1.20) reads explicitly as

$\psi_{p}(y)=\mbox{e}^{py}\ \exp\Big{(}-\tfrac{1}{r}\,\mbox{e}^{ry}\Big{)}\
{}_{1}F_{1}\Big{(}\tfrac{1}{2}+\tfrac{p}{r}-\tfrac{\lambda}{2},1+\tfrac{2p}{r},\tfrac{2}{r}\,\mbox{e}^{ry}\Big{)}\,.$
(2.15)

For $y\to+\infty$, it develops the asymptotic behaviour

$\psi_{p}(y)\asymp(2r)^{\frac{p}{r}+\frac{1}{2}}\
\frac{\Gamma(1+\tfrac{p}{r})}{2\sqrt{\pi}}\ \ D_{+}\
\Big{(}\frac{r}{2}\Big{)}^{\frac{1}{2}\lambda}\exp\Big{(}\tfrac{1}{r}\,\mbox{e}^{ry}-\big{(}\lambda+1\big{)}\,\tfrac{ry}{2}+o(1)\Big{)}\
,$ (2.16)

where

$D_{+}=\frac{\Gamma(\frac{1}{2}+\tfrac{p}{r})}{\Gamma\big{(}\tfrac{1}{2}+\tfrac{p}{r}-\tfrac{1}{2}\lambda(E)\big{)}}\
.$ (2.17)

The Jost solution $\psi_{p}(v)$ vanishes at large $v$ when $E$ is a zero of
$D_{+}=D_{+}(E)$. Hence, the spectrum is defined through the equation
$\lambda(E)=2m-1+\frac{2p}{r}$ with $m=1,2,3,\ldots\ $. This coincides with
(2.9), provided $\frac{p}{r}={\tt k}$ as in (1.19) with $n=r$.

###  2.2 Differential equation

For the case $q={\rm i}\,\mbox{e}^{\frac{{\rm i}\pi(r-n)}{2r(n+r)}}$ with
arbitrary $n>0$ we propose that the differential equation describing the
scaling limit of the Bethe roots for the ground state is given by

$\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-\sum_{j=0}^{A-1}c_{2j+1}\,E^{2j+1}\
\mbox{e}^{\big{(}\frac{n+r}{2}-\frac{n-r}{2r}\,(2j+1)\big{)}y}\,\bigg{]}\,\psi=0\
\ \ \ \ \big{(}A=\tfrac{r-1}{2}\big{)}\,.$ (2.18)

When all the coefficients $c_{2j+1}$ vanish, the ODE becomes that for the
${\cal Z}_{r}$ invariant case (1.18). Also, as $n=r$, it reduces to the
equation (2.14) where $\lambda(E)$ is given by (2.13), provided $c_{2j+1}$ are
set to $(-1)^{j}\,{\mathfrak{a}}_{2j+1}$. In general, these coefficients would
be some functions of the RG invariants of the form (1.22).

The ODE may be advocated for via the same types of arguments that were
originally developed in ref.[9] for the Schrödinger equation for the
anharmonic oscillator (which is equivalent to (1.18)). They are based on the
invariance of (2.18) w.r.t. the transformation

$\hat{\Omega}\ :\ y\mapsto y+\frac{2\pi{\rm i}}{n+r}\,,\qquad\qquad E\mapsto
q^{-2}\,E$ (2.19)

that holds true not only for $c_{2j+1}=0$, but any values of these
coefficients. Using this symmetry, one can derive the set of functional
relations for the connection coefficients of the ODE. The most fundamental one
is

$\mbox{e}^{\frac{2\pi{\rm
i}p}{n+r}}\,D_{+}(qE)\,D_{-}(q^{-1}E)-\mbox{e}^{-\frac{2\pi{\rm
i}p}{n+r}}\,D_{+}(q^{-1}E)\,D_{-}(qE)=2{\rm i}\sin\big{(}\tfrac{2\pi
p}{n+r}\big{)}\ ,$ (2.20)

which is obeyed by the spectral determinants $D_{\pm}(E)$. In order to define
the latter, apart from $\psi_{p}(y)$ (1.20), one should consider other
solutions of the differential equation, $\psi_{-p}(y)$ and $\chi(y)$.

For large negative $y$, the potential in (2.18) approaches $p^{2}$. As a
result, the ODE possesses a solution satisfying the asymptotic condition

$\displaystyle\psi_{-p}(y)\to\mbox{e}^{-py}\ \ \ \ \ \ \qquad{\rm as}\ \ \
\qquad y\to-\infty\qquad\qquad\big{(}\Re e(p)\leq 0\big{)}\,.$ (2.21)

It turns out that the analytic continuation of $\psi_{-p}(y)$ to the full
complex $p$ \- plane results in a meromorphic function of $p$. This allows one
to specify the solution $\psi_{-p}(y)$ for any complex $p$, except for the
discrete set of real, positive values
$p=\frac{1}{2}\,(n+r),\,(n+r),\,\frac{3}{2}\,(n+r),\ldots\ $. Similarly
$\psi_{p}(y)$, defined in the Introduction in formula (1.20), makes sense for
any $p\neq-\frac{1}{2}\,(n+r),\,-(n+r),\ldots\ $. The Wronskian
$W[\psi_{-p},\psi_{p}]\equiv\psi_{-p}\,\partial_{y}\psi_{p}-\psi_{p}\,\partial_{y}\psi_{-p}$
is given by

$W[\psi_{-p},\psi_{p}]=2p$ (2.22)

so that the two functions form a basis in the space of solutions of the ODE
(2.18) if $\frac{2p}{n+r}\notin\mathbb{Z}$.

The solution $\chi(y)$ is unambiguously defined through its asymptotic
behaviour

$\displaystyle\chi(y)\asymp\exp\Big{(}-\tfrac{n+r}{4}\ y-\tfrac{2}{n+r}\
\mbox{e}^{\frac{n+r}{2}y}+o(1)\Big{)}\qquad\qquad{\rm as}\qquad\qquad
y\to+\infty\qquad\qquad(n>r)\ .$ (2.23)

Note that, for now, we make the technical assumption $n>r$, while the case
$0<n<r$ will be discussed below. The spectral determinants are expressed via
the Wronskians of $\psi_{\pm p}(y)$ and $\chi(y)$ as

$\displaystyle D_{\pm}(E)=\frac{\sqrt{\pi}}{\Gamma(1\pm\frac{2p}{n+r})}\
(n+r)^{-\frac{1}{2}\mp\frac{2p}{n+r}}\ W[\chi,\psi_{\pm p}]\ .$ (2.24)

They are entire functions of $E$ and the overall multiplicative factor has
been chosen in such a way that

$D_{\pm}(0)=1\ .$ (2.25)

For the proof of the relation (2.20) one should expand $\chi$ in the basis of
solutions $\psi_{\pm p}$. It follows from (2.24) that

$\chi(y)=\frac{1}{\sqrt{\pi(n+r)}}\
\Big{(}\Gamma\big{(}-\tfrac{2p}{n+r}\big{)}\
(n+r)^{-\frac{2p}{n+r}}\,D_{-}(E)\,\psi_{+p}(y)+\Gamma\big{(}\tfrac{2p}{n+r}\big{)}\
(n+r)^{\frac{2p}{n+r}}\,D_{+}(E)\,\psi_{-p}(y)\Big{)}\,.$ (2.26)

Then one applies the symmetry transformation $\hat{\Omega}$ (2.19) to both
sides of this equation and computes the Wronskian of $\chi(y)$ and
$\hat{\Omega}\chi(y)$. Taking into account that
$W[\hat{\Omega}\chi,\chi]=2{\rm i}$, which is a consequence of (2.23), one
easily arrives at (2.20).

To make a link to the lattice system, recall that the inhomogeneous six-vertex
model possesses a commuting family of operators, a prominent member of which
is the Baxter $Q$ \- operator [26]. In fact, there are two of them [5, 6]. We
will use the notation $\mathbb{A}_{\pm}(\zeta)$ borrowed from ref.[25] and the
reader may consult sec. 3 of that work for their construction. Important
properties of $\mathbb{A}_{\pm}(\zeta)$ is that they commute amongst
themselves for different values of the spectral parameter,

$\displaystyle\big{[}\mathbb{A}_{\pm}(\zeta_{1}),\,\mathbb{A}_{\pm}(\zeta_{2})\big{]}=\big{[}\mathbb{A}_{\pm}(\zeta_{1}),\,\mathbb{A}_{\mp}(\zeta_{2})\big{]}=0\
,$ (2.27)

and satisfy the quantum Wronskian relation

$q^{+2\mathbb{P}}\
\mathbb{A}_{+}\big{(}q^{+1}\zeta\big{)}\,\mathbb{A}_{-}\big{(}q^{-1}\zeta\big{)}-q^{-2\mathbb{P}}\
\mathbb{A}_{-}\big{(}q^{+1}\zeta\big{)}\,\mathbb{A}_{+}\big{(}q^{-1}\zeta\big{)}=\big{(}q^{+2\mathbb{P}}-q^{-2\mathbb{P}}\,\big{)}\,f(\zeta)\
.$ (2.28)

Here

$f(\zeta)=\prod_{J=1}^{N}\big{(}1+\zeta/\eta_{J}\big{)}\,,$ (2.29)

while

$q^{2{\mathbb{P}}}=\mbox{e}^{{\rm i}\pi{\tt k}}\,q^{\mathbb{S}^{z}}$ (2.30)

with $\mathbb{S}^{z}$ being the $z$ \- projection of the total spin operator.
In addition, in the sector with given value of $S^{z}\geq 0$, the eigenvalues
of $\mathbb{A}_{+}(\zeta)$ are polynomials in $\zeta$ of order $N/2-S^{z}$.
The above immediately implies that the zeroes of these polynomials coincide
with the roots $\zeta_{j}$ of the Bethe Ansatz equations (1.11). The
eigenvalues of $\mathbb{A}_{-}(\zeta)$ are likewise polynomials in $\zeta$ of
order $N/2+S^{z}$, whose zeroes obey the equations similar to (1.11). The
relation (2.28) holds true for arbitrary inhomogeneities, though we focus on
the case where $\eta_{J}$ obey the $r$ \- site periodicity condition (1.14).

The Hamiltonian $\mathbb{H}$ associated with the inhomogeneous six-vertex
model belongs to the commuting family, i.e.,

$\big{[}\mathbb{H},\mathbb{A}_{\pm}(\zeta)\big{]}=0\,.$ (2.31)

Let $A_{\pm}(\zeta)$ be the eigenvalues of $\mathbb{A}_{\pm}(\zeta)$
corresponding to the ground state of the spin chain. In the ${\cal Z}_{r}$
invariant case the scaling relation for the Bethe roots (1.16) can be
equivalently written as

$\lim_{N\to\infty}A_{\pm}\Big{(}\big{(}N/(rN_{0})\big{)}^{-\frac{2n}{r(n+r)}}\,E\Big{)}=D_{\pm}(E)\qquad\qquad\qquad\qquad(n>r)$
(2.32)

with the constant $N_{0}$ defined in (1.17) and $D_{\pm}(E)$ are the spectral
determinants for the ODE (1.18). In the scaling limit, the quantum Wronskian
relation for $A_{\pm}(\zeta)$, which follows from (2.28), becomes the
functional relation (2.20). Motivated by the analysis of the free fermion
point, we propose to extend (2.32) to the case when $D_{\pm}$ are the spectral
determinants for the ODE (2.18). In taking the limit, the inhomogeneities
become $N$ \- dependent such that the values of

$\displaystyle{\mathfrak{a}}_{2j+1}=\frac{1}{2j+1}\
\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{1-\frac{2j+1}{r}}\ \frac{1}{r}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-2j-1}\ \ \ \ \
\qquad\big{(}j=0,\ldots,\tfrac{r-3}{2}\big{)}$ (2.33)

are kept fixed. Like in the free fermion case, ${\mathfrak{a}}_{r}$ can be
specified to be $\frac{1}{r}$ or, equivalently,333In this section, $r$ is
always taken to be odd so that the sign factor in the r.h.s. of (2.34) is one.
Nevertheless, keeping it there makes the formula applicable for any integer
$r$, which will be explored later.

$\frac{1}{r}\sum_{\ell=1}^{r}(\eta_{\ell})^{-r}=(-1)^{r-1}\,.$ (2.34)

As for the quantities ${\mathfrak{a}}_{s}$ (1.22) with $s$ even, we set them
to be zero:

$\sum_{\ell=1}^{r}(\eta_{\ell})^{-2j}=0\qquad\qquad\qquad\big{(}j=1,\ldots,\tfrac{r-1}{2}\big{)}\,.$
(2.35)

This can be thought of as taking the exponent $d_{s}$ in (1.22) to be
$+\infty$ for even $s$.

Our numerical work in support of the relation (2.32) went along the following
line. First of all we note that for given values of ${\mathfrak{a}}_{2j+1}$,
formulae (2.33)-(2.35) may be regarded as a system of equations for the
determination of the inhomogeneities. It is easy to show that
$(\eta_{\ell})^{-1}$ are zeroes of a single polynomial equation of order $r$,
whose coefficients are expressed in terms of ${\mathfrak{a}}_{2j+1}$. This
defines the set $\\{\eta_{\ell}\\}_{\ell=1}^{r}$ up to permutations, which has
no effect on the Bethe Ansatz equations. The roots corresponding to the ground
state of the Hamiltonian $\mathbb{H}$ can be obtained by continuously
deforming the solution to the Bethe Ansatz equations from the ${\cal Z}_{r}$
invariant case with ${\mathfrak{a}}_{2j+1}=0$ to the given values of the RG
invariants. Then one may consider the sums

$\displaystyle
h_{s}^{(N)}=s^{-1}\sum_{j=1}^{\frac{N}{2}}\big{(}\zeta_{j}\big{)}^{-s}\qquad\qquad(s=1,2,\ldots)\
.$ (2.36)

Taking into account that

$A_{+}(\zeta)=\prod_{j=1}^{\frac{N}{2}}\bigg{(}1-\frac{\zeta}{\zeta_{j}}\bigg{)}\,,$
(2.37)

the convergence of the l.h.s. of (2.32) implies the existence of the limits

$h_{s}^{(\infty)}=\lim_{N\to\infty}\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{\frac{2sn}{r(n+r)}}h_{s}^{(N)}\,.$
(2.38)

For $n>r$ this was indeed observed numerically. The spectral determinant
$D_{+}(E)$ occurring in the scaling limit of $A_{+}(\zeta)$ is an entire
function of $E$ which satisfies the normalization condition $D_{+}(0)=1$.
Hence the series

$\displaystyle\log D_{+}(E)=-\sum_{s=1}^{\infty}J_{s}\ E^{s}$ (2.39)

has a finite radius of convergence. This way, the scaling relation (2.32) is
equivalent to the infinite set of so-called sum rules:

$h_{s}^{(\infty)}=J_{s}\qquad\qquad\qquad(s=1,2,\ldots)\,.$ (2.40)

They can be used to provide numerical support for the ODE/IQFT correspondence
as was originally done in ref.[19] in the study of the scaling limit of the
Fateev-Zamolodchikov spin chain [27, 28, 29].

It is possible to calculate the first few coefficients $J_{s}$ analytically by
applying perturbation theory in $E$ to the ODE (2.18). From the results of the
work [30] one can show that

$\displaystyle J_{1}$ $\displaystyle=$ $\displaystyle c_{1}\rho_{1}\
f_{1}(h,g_{0})\ $ $\displaystyle J_{2}$ $\displaystyle=$
$\displaystyle\big{(}c_{1}\rho_{1}\big{)}^{2}\ f_{2}(h,g_{0})$ (2.41)
$\displaystyle J_{3}$ $\displaystyle=$
$\displaystyle\big{(}c_{1}\rho_{1}\big{)}^{3}\,f_{3}(h,g_{0})+c_{3}\rho_{3}\,f_{1}(h,g_{1})\
.$

The functions $f_{s}$ are given in the Appendix A, while

$\displaystyle\rho_{2j+1}=\frac{(n+r)^{2g_{j}-2}}{\Gamma^{2}(1-g_{j})}$ (2.42)

and

$\displaystyle h=\frac{p}{n+r}\ ,\ \ \ \ \ \ \ \ \ \ \
g_{j}=\frac{1}{2}-\frac{(2j+1)(n-r)}{2r(n+r)}\ .$ (2.43)

Since $h=\frac{\tt k}{2}$, see (1.19), formula (2.40) yields a prediction for
the dependence of $h_{s}^{(\infty)}$ on the twist parameter. This can be
checked even without the explicit knowledge of the relation between the RG
invariants ${\mathfrak{a}}_{2j+1}$ and the coefficients $c_{2j+1}$ of the
differential equation. Note that the expressions for $J_{s}$ for larger $s$
are not available in analytical form. Nevertheless, it may be verified that
the ${\tt k}$ \- dependence of both sides of (2.40) agrees for $s\geq 4$ by
computing $J_{s}$ via a numerical integration of the ODE (2.18).

Through the study of the cases $r=3,5,7$ we found that the relation between
the coefficients of the differential equation and the RG invariants takes the
general form

$\displaystyle c_{2j+1}=C_{j}^{(j)}\,{\mathfrak{a}}_{2j+1}+\sum_{k=2}^{A-j}\
\sum_{j_{1}+\ldots+j_{k}=j+(k-1)A\atop 0\leq j_{1},\ldots,j_{k}\leq A-1}\
C^{(j)}_{j_{1}\cdots j_{k}}\
{\mathfrak{a}}_{2j_{1}+1}\cdots{\mathfrak{a}}_{2j_{k}+1}\ .$ (2.44)

Here, for given $r$, $C^{(j)}_{j_{1}\cdots j_{k}}$ are some $n$ \- dependent
constants. In the case $r=n$, all $C^{(j)}_{j_{1}\cdots j_{k}}$ with $k\geq 2$
vanish, while $C_{j}^{(j)}=(-1)^{j}r$. In what follows, we present a
conjectured formula that expresses ${\mathfrak{a}}_{2j+1}$ in terms of
$c_{2j+1}$, i.e., the inversion of (2.44). However, this first requires a
discussion of the case $0<n<r$.

###  2.3 Relation between the RG invariants and the coefficients of the ODE

For $0<n\leq r$ the asymptotic condition (2.23), which is used to specify the
subdominant as $y\to+\infty$ solution of the differential equation, is not
literally applicable. The function $\chi(y)$ can be instead defined by means
of its WKB asymptotic, and the spectral determinants are then calculated, as
before, by means of formula (2.24). This out to be most efficient for the
numerical computation of $D_{\pm}(E)$. However, for our purposes, it would be
useful to introduce them for generic $n>0$ via analytic continuation. Assuming
$n>r$, we first derive the large - $E$ asymptotic behaviour of $D_{+}(E)$. The
latter has zeroes at $E_{m}^{(a)}$ defined through the condition (1.21):

$D_{+}\big{(}E_{m}^{(a)}\big{)}=0\,.$ (2.45)

These accumulate along the rays

$\displaystyle\arg(E)=-\frac{\pi}{r}\,A,\ \frac{\pi}{r}\,(2-A),\ \ldots,\
\frac{\pi}{r}\,\big{(}2\,(r-1)-A\big{)}\ \ \ \ \ \ \ \ \ \ \ ({\rm
mod}\,2\pi)\,.$ (2.46)

The rays are the Stokes lines that divide the complex $E$ \- plane into $r$
wedges, in which the large - $E$ asymptotic behaviour of the spectral
determinant is described differently. Inside each wedge, labelled by
$a=1,2,\ldots,r$, we can parameterize $E$ in terms of $\theta$ as

$\displaystyle E=(-1)^{r-1}\ \mbox{e}^{\frac{{\rm i}\pi}{r}(2a-1-A)}\
\mbox{e}^{\frac{2n\theta}{r(n+r)}}$ (2.47)

with444For odd $r$, which we consider in this section, one can neglect the
factor $(-1)^{r-1}$. It is included in eq. (2.47) since it will be applied for
any $r$, $A=0,1,\ldots,r-2$ and to make the notation consistent with that of
the work [31], devoted to the case $A=0$. The remark also carries over to eq.
(2.62), below.

$\big{|}\Im m(\theta)\big{|}<\frac{\pi(n+r)}{2n}\,.$ (2.48)

Then the standard WKB analysis yields

$\log D_{+}\asymp\frac{N_{0}\,\mbox{e}^{\theta}}{\cos(\frac{\pi
r}{2n})}-\sum_{j=0}^{A-1}\
\frac{D_{2j+1}}{\sin\big{(}\frac{\pi(2j+1)(n-r)}{2nr}\big{)}}\
\mbox{e}^{\frac{(2j+1)}{r}\,\big{(}\theta+{\rm
i}\pi(2a-1-A)\big{)}}-\frac{2np\,\theta}{r(n+r)}+\log C_{p}+o(1)$ (2.49)

as $\Re e(\theta)\to+\infty$. Here

$\displaystyle D_{2j+1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\pi}\Gamma\big{(}\tfrac{1}{2}-\tfrac{n-r}{2nr}(2j+1)\big{)}}{2n\Gamma\big{(}1-\tfrac{n-r}{2nr}(2j+1)\big{)}}\sum_{k=1}^{A}\mbox{e}^{-\frac{{\rm
i}\pi}{2}(r+1)(k-1)}\
\frac{\Gamma\big{(}m-1+\tfrac{n-r}{2nr}(2j+1)\big{)}}{m!\,\Gamma\big{(}\tfrac{n-r}{2nr}(2j+1)\big{)}}$
(2.50) $\displaystyle\times$
$\displaystyle\sum_{j_{1}+\ldots+j_{k}=j+(k-1)A\atop j_{1},\ldots,j_{k}\geq
0}c_{2j_{1}+1}\cdots c_{2j_{k}+1}$

and

$\displaystyle C_{p}=\sqrt{\frac{r}{n+r}}\
r^{\frac{2p}{r}}\,(n+r)^{-\frac{2p}{n+r}}\,\frac{\Gamma(1+\frac{2p}{r})}{\Gamma(1+\frac{2p}{n+r})}\,.$
(2.51)

For $0<n<r$ the spectral determinant can be defined to be the entire function
of $E$, whose zeroes coincide with $E_{m}^{(a)}$ (2.45) and which obeys the
asymptotic formula (2.49). This specifies $D_{+}(E)$ for any $n>0$ except for
a discrete set of values when the coefficients in the expansion (2.49) possess
simple poles:

$n=\frac{r\,(2j+1)}{2j+1+2wr}\qquad\qquad{\rm with}\qquad\qquad
j=0,\ldots,A-1;\ w=0,1,2,\ldots\ .$ (2.52)

The spectral determinant $D_{-}(E)$ may be introduced analogously. Together,
$D_{\pm}(E)$ satisfy the quantum Wronskian relation (2.20) for any $n>0$
except (2.52).555If we require the spectral determinants to be entire
functions of $E$, the quantum Wronskian relation (2.20) must be modified when
$n$ takes the values (2.52). The simplest illustration is provided by the case
$r=n$, where $D_{+}(E)$ is given by (2.17), while $D_{-}(E)$ is obtained from
the former via the substitution $p\mapsto-p$. We will exclude such exceptional
$n$ from our analysis to avoid making the discussion too technical (see, e.g.,
ref.[39] for a treatment of the cases $r=1,2$).

The scaling limit of the individual Bethe roots,

$\lim_{N\to\infty\atop m-{\rm
fixed}}\,\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{\frac{2n}{r(n+r)}}\zeta_{m}^{(a)}\,,$
(2.53)

exists for any positive $n$. However as $0<n\leq r$, it turns out there are
issues with the existence of the limit (2.38). To properly define
$h_{s}^{(\infty)}$ certain subtractions need to be made so that the large -
$N$ limit can be taken. Namely, introduce the “regularized” version of the sum
$h_{s}^{(N)}$:

$\displaystyle h_{s}^{(N,{\rm
reg})}=s^{-1}\sum_{m=1}^{M}(\zeta_{m})^{-s}+\frac{(-1)^{s-1}N}{2sr\cos(s\gamma)}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-s}\ ,$ (2.54)

where $M$ stands for the total number of Bethe roots and in the case of the
ground state $M=\frac{N}{2}$. Without going into details, we mention that the
form of the counterterm is motivated by the Bethe Ansatz equations (1.11) with
the inhomogeneities subject to the $r$ \- site periodicity condition
$\eta_{J+r}=\eta_{J}$. Taking into account the definition of
${\mathfrak{a}}_{s}$ (1.22), the above formula can be re-written as

$\displaystyle h_{s}^{(N,{\rm
reg})}=s^{-1}\sum_{m=1}^{M}(\zeta_{m})^{-s}+\frac{(-1)^{s-1}\
rN_{0}{\mathfrak{a}}_{s}}{2\cos(s\gamma)}\
\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{1-d_{s}}\ .$ (2.55)

This will be used to define $h_{s}^{(N,{\rm reg})}$ for any value of
$A=0,1,2,\ldots,r-2$. In the case we are currently considering, where $r$ is
odd and $A=\frac{r-1}{2}$, all ${\mathfrak{a}}_{2j}$ with $j=1,2,\ldots,A$ are
set to zero, so that

$h_{2j}^{(N,{\rm reg})}=h_{2j}^{(N)}\,.$ (2.56)

At the same time, the values of ${\mathfrak{a}}_{2j+1}$ with
$j=0,1,\ldots,A-1$ are generic complex numbers, which are held fixed in taking
the scaling limit, and the exponent $d_{2j+1}$ in (2.55) is given by

$\displaystyle d_{2j+1}=1-\frac{2j+1}{r}\ \ \ \ \ \ \qquad\ \ \ \
(j=0,1,\ldots,A-1)\,.$ (2.57)

Hence,

$\displaystyle\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{\frac{2(2j+1)n}{r(n+r)}}\,h_{2j+1}^{(N,{\rm
reg})}=\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{\frac{2(2j+1)n}{r(n+r)}}\,h_{2j+1}^{(N)}+\frac{(-1)^{j}\,rN_{0}\,\mathfrak{a}_{2j+1}}{2\sin\big{(}\frac{\pi(2j+1)(n-r)}{2r(n+r)}\big{)}}\
\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{\frac{(2j+1)(n-r)}{r(n+r)}}\ .$ (2.58)

For $n>r$ the second term in the r.h.s. tends to zero as $N\to\infty$ and can
be neglected. However, as $0<n<r$ it grows and cancels the divergent behaviour
of the first term, ensuring the existence of the limit

$h_{s}^{(\infty)}=\lim_{N\to\infty}\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{\frac{2sn}{r(n+r)}}h_{s}^{(N,{\rm
reg})}\,.$ (2.59)

Note that even when $n>r$, the last formula is more efficient than (2.38) for
numerical purposes since the inclusion of the counterterm greatly improves the
convergence. For instance, we observed

$\displaystyle\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{\frac{2n}{r(n+r)}}\,h_{1}^{(N,{\rm
reg})}=h_{1}^{(\infty)}+O\big{(}N^{-2}\big{)}\ \quad\ \ \ (n>r-2)\,,$ (2.60)

while for $0<n<r-2$ the correction term goes to zero as $N\to\infty$ with a
certain exponent which depends on $n$.

Closely related to the divergent behaviour of $h_{s}^{(N)}$ for $0<n\leq r$ is
that some of the coefficients in the expansion of the spectral determinant
(2.39) possess a simple pole at $n=r$. It is worth discussing this in some
detail since it turns out to give a hint of how to obtain the relations
between the RG invariants and the coefficients of the differential equation.
Since $E_{m}^{(a)}$ are zeroes of $D_{+}(E)$, it follows that $J_{s}$ can be
written in the form

$\displaystyle J_{s}=\frac{1}{s}\
\sum_{m=1}^{\Lambda}\sum_{a=1}^{r}\big{(}E_{m}^{(a)}\big{)}^{-s}+\Xi_{s}(\Lambda)\,.$
(2.61)

Here the remainder $\Xi_{s}(\Lambda)\to 0$ as the cutoff $\Lambda\to\infty$
for $n>r$. Obtaining the large - $\Lambda$ behaviour of $\Xi_{s}(\Lambda)$
requires knowledge of the asymptotics of the zeroes $E_{m}^{(a)}$ as $m\gg 1$.
The latter can be extracted from the ODE (2.18) using the WKB approximation.

The asymptotic formula (2.49) for $D_{\pm}(E)$ holds true inside the wedge in
the complex $E$ \- plane described by eqs. (2.47), (2.48) and can not be
applied along its boundary. At the boundary of the $a^{\rm th}$ and
$(a+1)^{\rm th}$ (mod $r$) wedges, $E$ can be parameterized by real positive
$\theta$ as

$\displaystyle E=(-1)^{r-1}\ \mbox{e}^{\frac{{\rm i}\pi}{r}(2a-A)}\
\mbox{e}^{\frac{2n\theta}{r(n+r)}}\ \ \ \ \qquad\qquad(\theta>0)\,.$ (2.62)

Along this ray, the asymptotic behaviour of $D_{+}(E)$ is given by

$\displaystyle D_{+}\asymp 2C_{p}\ \mbox{e}^{-\frac{2np\theta}{r(n+r)}}\
\mbox{e}^{-\chi_{a}(\theta)}\cos\big{(}\phi_{a}(\theta)\big{)}\ \ \ \ \ \ \ \
\qquad(\theta\to+\infty)\ ,$ (2.63)

where

$\displaystyle\chi_{a}(\theta)=N_{0}\tan\big{(}\tfrac{\pi
r}{2n}\big{)}\,\mbox{e}^{\theta}+\sum_{j=0}^{A-1}\ D_{2j+1}\
\cot\big{(}\tfrac{\pi(2j+1)(n-r)}{2nr}\big{)}\
\mbox{e}^{\frac{(2j+1)}{r}\,\big{(}\theta+{\rm i}\pi(2a-A)\big{)}}+o(1)$
(2.64)

and

$\displaystyle\phi_{a}(\theta)=N_{0}\mbox{e}^{\theta}+\sum_{j=0}^{A-1}\mbox{e}^{\frac{{\rm
i}\pi}{r}(2a-A)(2j+1)}\ D_{2j+1}\ \mbox{e}^{\frac{2j+1}{r}\,\theta}-\tfrac{\pi
p}{r}+o(1)\,.$ (2.65)

In view of eq. (2.62) we write the zeroes $E_{m}^{(a)}$ in the form

$E_{m}^{(a)}=(-1)^{r-1}\ \mbox{e}^{\frac{{\rm i}\pi}{r}(2a-A)}\
\exp\Big{(}\tfrac{2n}{r(n+r)}\,\theta_{m}^{(a)}\Big{)}\ .$ (2.66)

Then $\theta_{m}^{(a)}$ solves the equation

$\displaystyle\phi_{a}\big{(}\theta_{m}^{(a)}\big{)}=\pi\big{(}m-\tfrac{1}{2}\big{)}\
\ \ \ \ \ \ \ \ \ (m=1,2,\ldots)\ ,$ (2.67)

which is similar to the usual Bohr-Sommerfeld quantization condition. It is
convenient to re-write the latter in a different way. Introduce $X$ and $Y$
such that

$\displaystyle X=-{\rm i}\,\mbox{e}^{\frac{{\rm i}\pi}{r}(2a-A)}\
\mbox{e}^{\frac{1}{r}\,\theta_{m}^{(a)}}\ ,\ \ \ \ Y=-{\rm
i}\,\mbox{e}^{\frac{{\rm i}\pi}{r}(2a-A)}\ \bigg{(}\frac{\pi}{N_{0}}\
\Big{(}m-\frac{1}{2}+\frac{p}{r}\,\Big{)}\bigg{)}^{\frac{1}{r}}$ (2.68)

as well as

$\displaystyle G_{r-2j-1}=(-1)^{j}\ \frac{D_{2j+1}}{N_{0}}\ .$ (2.69)

Formula (2.67), upon dropping the term $o(1)$ in (2.65), then takes the form

$\displaystyle X^{r}+\sum_{m=1}^{A}G_{2m}X^{r-2m}={Y}^{r}\ ,$ (2.70)

which can be considered as a polynomial equation to determine $X$ given $Y$.
Since $m$ is assumed to be large it follows from the definition (2.68) that
$|Y|\gg 1$. Moreover, eq. (2.70) has $r$ roots and we should focus on the one,
$X(Y)$, such that

$\displaystyle X(Y)=Y+o(1)\ \ \ \ \ {\rm as}\ \ \ \ Y\to\infty\ .$ (2.71)

The Lagrange formula [32, 33] (see also Appendix B) gives this root as a
series in inverse powers of $Y$:

$\displaystyle X(Y)=Y+\sum_{j=0\atop 2j+1\not=0\ {\rm
mod}\,r}^{\infty}R_{2j+1}\ Y^{-2j-1}$ (2.72)

with

$R_{2j+1}=\frac{1}{r}\sum_{\alpha_{1}+2\alpha_{2}+\ldots+(j+1)\alpha_{j+1}=j+1\atop\alpha_{1},\ldots,\alpha_{j+1}\geq
0}\frac{(-1)^{\alpha_{1}+\ldots+\alpha_{j+1}}}{\alpha_{1}!\alpha_{2}!\ldots\alpha_{j+1}!}\
\frac{\Gamma(\alpha_{1}+\ldots+\alpha_{j+1}-\frac{2j+1}{r})}{\Gamma(1-\frac{2j+1}{r})}\
G_{2}^{\alpha_{1}}G_{4}^{\alpha_{2}}\dots G_{2j+2}^{\alpha_{j+1}}$ (2.73)

In the case at hand, as it follows from (2.69) and (2.50), $G_{2m}$ is given
by

$\displaystyle G_{2m}$ $\displaystyle=$ $\displaystyle(-1)^{A-m}\
\frac{\Gamma\big{(}\frac{3}{2}+\frac{r}{2n}\big{)}\Gamma\big{(}\frac{1}{2}-\frac{(r-2m)(n-r)}{2nr}\big{)}}{\Gamma\big{(}\frac{r}{2n}\big{)}\Gamma\big{(}1-\frac{(r-2m)(n-r)}{2nr}\big{)}}\
\sum_{k=1}^{m}(-1)^{(A+1)(k-1)}\
\frac{\Gamma\big{(}k-1+\tfrac{n-r}{2nr}(r-2m)\big{)}}{m!\,\Gamma\big{(}\tfrac{n-r}{2nr}(r-2m)\big{)}}$
(2.74) $\displaystyle\times$ $\displaystyle\sum_{j_{1}+\ldots+j_{k}=kA-m\atop
j_{1},\ldots,j_{k}\geq 0}c_{2j_{1}+1}\cdots c_{2j_{k}+1}\ .$

Combining the above with (2.66) yields the asymptotic formula for the zeroes:

$E_{m}^{(a)}\asymp\mbox{e}^{\frac{{\rm i}\pi}{r}(2a-A)}\
|Y_{m}|^{\frac{2n}{n+r}}\ \bigg{(}1+\sum_{k=1}^{A}R_{2k-1}\
\mbox{e}^{-\frac{{\rm i}\pi}{r}(4a+1)k}\
|Y_{m}|^{-2k}+O\big{(}m^{-1}\big{)}\bigg{)}^{\frac{2n}{n+r}}\qquad(m\gg 1)\ ,$
(2.75)

where

$|Y_{m}|=\bigg{(}\frac{\pi}{N_{0}}\,\Big{(}m-\frac{1}{2}+\frac{p}{r}\Big{)}\bigg{)}^{\frac{1}{r}}\,.$
(2.76)

The remainder $\Xi_{s}(\Lambda)$ in (2.61) involves a sum over negative
integer powers of $E_{m}^{(a)}$. They are expressed through $Q_{2k-1}(z)$,
which are the coefficients in the formal expansion

$\displaystyle\bigg{(}1+\sum_{k=1\atop 2k-1\not=0\ {\rm
mod}\,r}^{\infty}R_{2k-1}\
Y^{-2k}\bigg{)}^{z}=1+\sum_{k=1}^{\infty}Q_{2k-1}(z)\ Y^{-2k}\ .$ (2.77)

Then

$\displaystyle\big{(}E_{m}^{(a)}\big{)}^{-s}\asymp\mbox{e}^{-\frac{{\rm
i}\pi}{r}(2a-2-A)s}\ |Y_{m}|^{-s\nu}\ \bigg{(}1+\sum_{k=1}^{A}Q_{2k-1}(-s\nu)\
\mbox{e}^{-\frac{{\rm i}\pi}{r}(4a+1)k}\
|Y_{m}|^{-2k}+O\big{(}m^{-1}\big{)}\bigg{)}\,,$ (2.78)

where we use

$\displaystyle\nu=\frac{2n}{n+r}\ .$ (2.79)

Taking $s$ to be odd and performing the sum over $a$ leads to

$\sum_{a=1}^{r}\big{(}E_{m}^{(a)}\big{)}^{-s}\asymp r\,\mbox{e}^{\frac{{\rm
i}\pi}{2}(s-1)}\ Q_{r-s-1}(-s\nu)\
|Y_{m}|^{-r+s(1-\nu)}+O\Big{(}m^{-2+(1-\nu)\frac{s}{r}}\Big{)}\qquad(s=2j+1)\
.$ (2.80)

It follows that the large $\Lambda$ behaviour for the remainder in (2.61) is
given by

$\displaystyle\Xi_{2j+1}(\Lambda)=\frac{(-1)^{j}\,rN_{0}}{\big{(}\frac{\pi(2j+1)(n-r)}{r(n+r)}\big{)}}\
\,\frac{Q_{r-2j-2}\big{(}-(2j+1)\nu\big{)}}{2j+1}\
\bigg{(}\frac{N_{0}}{\pi\Lambda}\bigg{)}^{\frac{(2j+1)(n-r)}{r(n+r)}}+O\bigg{(}\Lambda^{-1-\frac{(2j+1)(n-r)}{r(n+r)}}\bigg{)}\,.$
(2.81)

Let’s compare the residues of (2.81) and (2.58) at the pole $n=r$. It leads to
the relation

$\mathfrak{a}_{2j+1}=\frac{1}{2j+1}\ Q_{r-2j-2}\big{(}-(2j+1)\nu\big{)}\ ,$
(2.82)

or, explicitly, taking into account the definition (2.77),

$\displaystyle\mathfrak{a}_{2j+1}=\frac{1}{2j+1}\ \sum_{k=1}^{A-j}(-1)^{k}\
\frac{\Gamma\big{(}(2j+1)\nu+k\big{)}}{k!\Gamma\big{(}(2j+1)\nu\big{)}}\
\sum_{j_{1}+j_{2}\ldots+j_{k}=A-j-k\atop j_{1},\ldots,j_{k}\geq
0}R_{2j_{1}+1}\ldots R_{2j_{k}+1}\ ,$ (2.83)

where $R_{2j+1}$ is given by eqs. (2.73) and (2.74). Of course, the above is
expected to be valid only for $n=r$. Nevertheless, we found from the numerical
work that it holds true for any $n>0$. The relation can be inverted, whereupon
it takes the form (2.44). The coefficients $C^{(j)}_{j_{1}\cdots j_{k}}$ are
rather cumbersome in general and the simplest of them is given by

$C_{j}^{(j)}=(-1)^{j}\,\frac{r\Gamma\big{(}\frac{r}{2n}\big{)}\Gamma\big{(}1-\frac{(2j+1)(n-r)}{2rn}\big{)}}{\Gamma\big{(}\frac{1}{2}+\frac{r}{2n}\big{)}\Gamma\big{(}\frac{1}{2}-\frac{(2j+1)(n-r)}{2rn}\big{)}}\
.$ (2.84)

Formula (2.44) is written out explicitly in the Appendix C for the cases
$r=3,5,7$. Numerical work confirms its validity.

###  2.4 Some remarks about the scaling limit

In our discussion we focused on describing the large - $N$ behaviour of the
Bethe roots in the vicinity of $\zeta=0$. A similar analysis can be performed
for the Bethe roots located near $\zeta=\infty$. In this case, (1.16) is
replaced by

$\bar{E}_{m}^{(a)}=\lim_{N\to\infty\atop m-{\rm
fixed}}\,\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{\frac{2n}{r(n+r)}}\Big{(}\zeta_{M_{a}-m}^{(a)}\Big{)}^{-1}\,,$
(2.85)

where for the ground state in the sector $S^{z}=0$ all $M_{a}$ are equal and
coincide with $N/(2r)$. If the scaling limit is taken with
$\mathfrak{a}_{2j+1}$ (2.33) kept fixed and the conditions (2.34), (2.35) are
imposed on the inhomogeneities, it turns out that $\bar{E}_{m}^{(a)}$ are
described by the differential equation corresponding to the ${\cal Z}_{r}$
invariant case similar to (1.18). This can be observed by plotting the Bethe
roots in the complex $\beta$ \- plane with $\beta=-\frac{1}{2}\log(\zeta)$,
which is done in the left panel of fig. 1 for the case $r=3$. One sees that
for the roots on the right side of the distribution, for which $|\zeta_{j}|\ll
1$, there is significant deviation from the lines $\Im
m(\beta)=\pm\frac{\pi}{6},\,\frac{\pi}{2}$ so that their scaling behaviour is
described via the differential equation (2.18) with $r=3$. The Bethe roots
depicted on the left side of the figure have $|\zeta_{j}|\gg 1$ and are
visibly located along the lines with $\Im
m(\beta)=\pm\frac{\pi}{6},\,\frac{\pi}{2}$. Their scaling limit turns out to
be same as that of the Bethe roots for the ${\cal Z}_{3}$ invariant case.

$\beta$$\beta$$\Im m(\beta)=\frac{\pi}{2}$$\Im m(\beta)=+\frac{\pi}{6}$$\Im
m(\beta)=-\frac{\pi}{6}$$\Im m(\beta)=\frac{\pi}{2}$$\Im
m(\beta)=+\frac{\pi}{6}$$\Im m(\beta)=-\frac{\pi}{6}$ Figure 1: Depicted are
plots of the Bethe roots in the complex plane $\beta=-\frac{1}{2}\log(\zeta)$
for $r=3$, $N=120$, $n=5$ and ${\tt k}=0.05$. For the left panel, the
inhomogeneities are taken to satisfy the conditions (2.33) - (2.35) with
$\mathfrak{a}_{1}=0.4$. The right panel corresponds to the case when the spin
chain possesses ${\cal C}{\cal P}$ and ${\cal T}$ symmetry. The
inhomogeneities are specified according to (2.96) - (2.98) with
$\mathfrak{b}_{1}=0.4$.

Keeping in mind the invariance of the Bethe Ansatz equations (1.11) w.r.t. the
simultaneous inversion $\zeta_{j}\mapsto\zeta_{j}^{-1}$,
$\eta_{\ell}\mapsto(\eta_{\ell})^{-1}$ and $\omega\mapsto\omega^{-1}$ (${\tt
k}\mapsto-{\tt k}$), it is clear how to organize the scaling limit such that
the scaling behaviour of the roots near $\zeta=0$ and $\zeta=\infty$ is
reversed. One should treat

$\bar{{\mathfrak{a}}}_{s}=\frac{1}{s}\
\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{\bar{d}_{s}}\ \frac{1}{r}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{s}$ (2.86)

with $s=2j+1$ and $\bar{d}_{2j+1}=1-\frac{2j+1}{r}$ as the RG invariants and
impose the additional restrictions on the inhomogeneities:

$\sum_{\ell=1}^{r}(\eta_{\ell})^{2j}=0\qquad\qquad(j=1,2,\ldots,A)$ (2.87)

and

$\frac{1}{r}\sum_{\ell=1}^{r}(\eta_{\ell})^{r}=(-1)^{r-1}\,.$ (2.88)

The numbers $\bar{E}_{m}^{(a)}$ (2.85) would then be described by the “barred”
version of the ODE (2.18):

$\bigg{[}-\partial_{\bar{y}}^{2}+\bar{p}^{2}+\mbox{e}^{(n+r)\bar{y}}-(-1)^{A}\,\bar{E}^{r}\
\mbox{e}^{r\bar{y}}-\sum_{j=0}^{A-1}\bar{c}_{2j+1}\,\bar{E}^{2j+1}\
\mbox{e}^{\big{(}\frac{n+r}{2}-\frac{n-r}{2r}\,(2j+1)\big{)}\bar{y}}\,\bigg{]}\,\bar{\psi}=0\
\ \ \ \ \big{(}A=\tfrac{r-1}{2}\big{)}\,.$ (2.89)

Here $\bar{p}=-p=-\frac{n+r}{2}\ {\tt k}$, while formula (2.44) with the
replacements
$(\mathfrak{a}_{2j+1},c_{2j+1})\mapsto(\bar{\mathfrak{a}}_{2j+1},\bar{c}_{2j+1})$
gives the relation between the RG invariants and the coefficients of the
differential equation.

In view of applications to QFT, of special interest are the spin chains
admitting extra global symmetries — charge conjugation ${\cal C}$, parity
${\cal P}$ and time reversal ${\cal T}.$ If the inhomogeneities obey the
condition

$\eta_{\ell}=\big{(}\eta_{r+1-\ell}\big{)}^{-1}\qquad\qquad(\ell=1,2,\ldots,r)$
(2.90)

then the system possesses ${\cal CP}$ invariance, i.e.,

$\hat{{\cal C}}\hat{{\cal P}}\,\mathbb{H}\,\hat{{\cal C}}\hat{{\cal
P}}=\mathbb{H}\,.$ (2.91)

For the action of the involutory operators $\hat{{\cal C}}$ and $\hat{{\cal
P}}$ on the space of states of the spin chain we refer the reader to sec. 4 of
the paper [25]. In the latter, it is also described how the ${\cal C}{\cal P}$
transformation acts on the Baxter $Q$ \- operators:

$\displaystyle{\hat{\cal C}}{\hat{\cal
P}}\,{\mathbb{A}}_{\pm}(\zeta)\,{\hat{\cal C}}{\hat{\cal P}}$ $\displaystyle=$
$\displaystyle\zeta^{\frac{N}{2}\mp\,\mathbb{S}^{z}}\
{\mathbb{A}}_{\mp}\big{(}\zeta^{-1}\big{)}\,\big{[}\,\mathbb{A}^{(\infty)}_{\pm}\,\big{]}^{-1}\,,$
(2.92)

where

${\mathbb{A}}^{(\infty)}_{\pm}=\lim_{\zeta\to\infty}\zeta^{-\frac{N}{2}\pm\mathbb{S}^{z}}\mathbb{A}_{\pm}(\zeta)\,.$
(2.93)

Note that $\hat{\cal C}\hat{\cal P}$ intertwines the sectors with $S^{z}$ and
$-S^{z}$. If, in addition to (2.90), the inhomogeneities obey the reality
condition

$\eta_{\ell}^{*}=\eta_{\ell}^{-1}\,,$ (2.94)

then the spin chain possesses time reversal symmetry ${\cal T}$ (again see
ref.[25] for details). This way, the requirements (2.90) and (2.94) will
guarantee that the field theory underlying the critical behaviour of the spin
chain is ${\cal CP}$ and ${\cal T}$ invariant. Also note that, without loss of
generality, the inhomogeneities can be chosen in such a way to satisfy the
condition

$\prod_{\ell=1}^{r}\eta_{\ell}=1$ (2.95)

even for the spin chain, where no other constraints on $\eta_{\ell}$ are being
assumed.

For the ${\cal Z}_{r}$ invariant case,
$\eta_{\ell}=(-1)^{r}\,\mbox{e}^{\frac{{\rm i}\pi}{r}\,(2\ell-1)}$ and the
conditions (2.90), (2.94) and (2.95) are satisfied. For the model, which
possesses ${\cal C}{\cal P}$ and ${\cal T}$ symmetries but ${\cal Z}_{r}$
invariance is broken, the inhomogeneities can be parameterized through the
phases $\delta_{\ell}$ as

$\eta_{\ell}=(-1)^{r}\mbox{e}^{\frac{{\rm i}\pi}{r}(2\ell-1)+{\rm
i}\delta_{\ell}}$ (2.96)

with

$\delta_{\ell}=\delta_{\ell}^{*}=-\delta_{r+1-\ell}\,,\qquad\qquad\delta_{1}+\ldots+\delta_{r}=0\,.$
(2.97)

In the case when $r$ is odd and $A=\frac{r-1}{2}$, one can consider the limit
with $\delta_{\ell}\to 0$ as $N\to\infty$ such that

$\displaystyle\delta_{\ell}=-\delta_{r+1-\ell}=2\,\sum_{j=0}^{A-1}\
\sin\big{(}\tfrac{\pi}{r}(2\ell-1)(2j+1)\big{)}\ {\mathfrak{b}}_{2j+1}\
\bigg{(}\frac{rN_{0}}{N}\bigg{)}^{1-\frac{2j+1}{r}}\ \ \ (\ell=1,\ldots,A+1)$
(2.98)

keeping the real numbers ${\mathfrak{b}}_{2j+1}$ fixed. The Bethe roots for
the ground state of the spin chain with $N=120$ and $\mathfrak{b}_{1}=0.4$ are
depicted in the right panel of fig. 1, which may be compared with the plot in
the left panel. Our numerical work for the cases $r=3,5,7$ show that the Bethe
roots still possess the scaling behaviour (1.16) and (2.85) with $E_{m}^{(a)}$
and $\bar{E}_{m}^{(a)}$ described in terms of the differential equations
(2.18) and (2.89), respectively. If we assume that the scaling limit exists
then formula (2.92) implies that the coefficients $c_{2j+1}$ must coincide
with their barred counterparts:

$\displaystyle\bar{c}_{2j+1}={c}_{2j+1}\ \ \ \ \qquad\ (j=0,\ldots,A-1)\ .$
(2.99)

The relation between $\mathfrak{b}_{2j+1}$ and $c_{2j+1}$ turns out to take
the form

$\displaystyle c_{2j+1}={C}_{j}^{(j)}\,{\mathfrak{b}}_{2j+1}+\ldots\,,$
(2.100)

where the ellipses involve a sum over monomials in
${\mathfrak{b}}_{2j^{\prime}+1}$ with $j^{\prime}<j$, each of which is
quadratic or of higher order in the RG invariants. The coefficient for the
linear term coincides with ${C}_{j}^{(j)}$ from (2.44) and is given by (2.84).
It is worthing noting that though the $N\to\infty$ limit (2.59) involving the
regularized sum over the Bethe roots $h_{s}^{(N,{\rm reg})}$ still exists and
is described by (2.40), the rate of convergence is considerably slower.
Moreover, the definition of $h_{s}^{(N,{\rm reg})}$ requires different
subtractions, depending on the RG invariant $\mathfrak{b}_{2j+1}$ but not on
${\tt k}$, in order to ensure the existence of the limit for any $n>0$. In
particular, when the scaling limit was defined as in sec. 2.2,
$h_{2j}^{(N,{\rm reg})}=h_{2j}^{(N)}$, see (2.56), while for the case (2.98)
this is no longer true.

Thus we see that there is freedom in the way in which one takes the scaling
limit — different prescriptions give rise to the ODEs of the same form (2.18)
and (2.89) which describe the scaling of the Bethe roots for the vacua. The
coefficients of the differential equations depend on the set of RG invariants
that specify the scaling procedure. Note that the number of non-vanishing RG
invariants may not necessarily coincide with the number of coefficients
$c_{2j+1}$. For instance, we discussed the scaling limit, where
$\sum_{\ell=1}^{r}(\eta_{\ell})^{-2j}$ were set to zero for $j=1,2,\ldots,A$.
As was mentioned before, this formally corresponds to taking the exponent
$d_{2j}$ in (1.22) to be $+\infty$. However, if the condition is relaxed and
one considers the limit with non-vanishing RG invariants $\mathfrak{a}_{2j}$
but $d_{2j}$ sufficiently large positive numbers, the ODEs describing the
scaling limit remain the same. In particular for the cases $r=3,5$ we found
that if the exponent $d_{2j}>\frac{1}{r}$, the presence of the non-zero RG
invariants $\mathfrak{a}_{2j}$ do not essentially affect the scaling
behaviour.

###  2.5 Low energy spectrum in the scaling limit

In this paper, we focus on describing the scaling of the Bethe roots for the
ground state of the spin chain, which may be regarded as a rather formal
mathematical problem. From the physical point of view the most interesting
question would be the study of the field theory underlying the critical
behaviour of the lattice system. While this goes beyond the scope of our work,
here we nevertheless give a brief discussion of the spectrum of the
Hamiltonian $\mathbb{H}$ in the scaling limit in the regime with
$\frac{\pi}{2}(1-\frac{1}{r})<\gamma<\frac{\pi}{2}(1+\frac{1}{r})$ and $r$
odd. The analysis is mostly informed by a study of the free fermion point with
$\gamma=\frac{\pi}{2}$.

####  2.5.1 Low energy spectrum for the ${\cal Z}_{r}$ invariant spin chain

For a given solution of the Bethe Ansatz equations (1.11), with the
inhomogeneities obeying the $r$ \- site periodicity condition (1.12), the
corresponding eigenvalue of the Hamiltonian $\mathbb{H}$ reads as

${\cal E}=2{\rm
i}\sum_{\ell=1}^{r}\sum_{m=1}^{N/2-S^{z}}\,\bigg{(}\frac{1}{1+\zeta_{m}\,q^{-1}\,/\eta_{\ell}}-\frac{1}{1+\zeta_{m}\,q/\eta_{\ell}}\bigg{)}\
.$ (2.101)

The Hamiltonian commutes with the $r$ \- site translation operator
$\mathbb{K}$, whose eigenvalue ${\cal K}$ is computed from the Bethe roots as
(see, e.g., sec. 6 of ref.[25] for details)

${\cal K}=\mbox{e}^{{\rm i}r\pi{\tt
k}}\,q^{-r(\frac{N}{2}-S^{z})}\,\prod_{\ell=1}^{r}\prod_{m=1}^{N/2-S^{z}}\frac{\zeta_{m}+\eta_{\ell}\,q^{+1}}{\zeta_{m}+\eta_{\ell}\,q^{-1}}\
.$ (2.102)

The lowest energy state in the sector $S^{z}=0$ is translationally invariant,
i.e., ${\cal K}_{\rm GS}=1$. As was mentioned in the Introduction, in the
${\cal Z}_{r}$ invariant case with

$\eta_{\ell}=(-1)^{r}\,\mbox{e}^{\frac{{\rm i}\pi}{r}\,(2\ell-1)}\ \ \
\qquad\qquad(\ell=1,\ldots,r)$ (2.103)

the corresponding Bethe roots are simply related to those for the ground state
of the XXZ spin chain. Using well known facts about the latter, it is
straightforward to deduce the large - $N$ behaviour for the ground state
energy:

${\cal E}^{(0)}_{\rm GS}=Ne_{\infty}+\frac{2\pi rv_{{\rm
F}}}{N}\,\Big{(}\,\frac{1}{2}\,(n+r)\,{\tt
k}^{2}-\frac{r}{12}\Big{)}+o(N^{-1})\,,$ (2.104)

where the specific bulk energy $e_{\infty}$ and the Fermi velocity $v_{\rm F}$
are given by

$\displaystyle e_{\infty}=-\frac{2v_{{\rm F}}}{\pi}\,\int_{0}^{\infty}{\rm
d}t\
\frac{\sinh\big{(}\frac{rt}{n}\big{)}}{\sinh\big{(}\frac{(n+r)\,t}{n}\big{)}\,\cosh(t)}\,,\qquad\qquad\qquad
v_{{\rm F}}=\frac{r\,(n+r)}{n}\ .$ (2.105)

The superscript “$\scriptstyle(0)$” has been introduced in order to emphasize
that the above formulae are valid for the ${\cal Z}_{r}$ invariant spin chain.
Eq. (2.104) leads one to expect that in the scaling limit $\mathbb{H}^{(0)}$
can be written as

$\mathbb{H}^{(0)}\asymp Ne_{\infty}+\frac{2\pi rv_{{\rm F}}}{N}\,\hat{H}_{\rm
CFT}+o(N^{-1})\,.$ (2.106)

Here $\hat{H}_{\rm CFT}$ stands for the Hamiltonian of the underlying CFT,

$\hat{H}_{\rm
CFT}=\int_{0}^{2\pi}\frac{\mbox{d}x}{2\pi}\,\big{(}T+\overline{T}\,\big{)}\,,$
(2.107)

with $T$ and $\overline{T}$ being the chiral components
($\partial\overline{T}=\bar{\partial}T=0$) of the canonically normalized
energy-momentum tensor.

For the XXZ spin chain (1.1), corresponding to $r=1$, the low energy states
for any value of the anisotropy parameter $-1<\Delta<1$ can be classified in
the same way as for the free fermion point with $\Delta=0$. Namely, one should
choose the basis of eigenstates of the Hamiltonian to be the Bethe states
$\Psi_{N}$, which are labelled by the solution of the Bethe Ansatz equations.
In the scaling limit $\Psi_{N}$ takes the form of the tensor product of the
chiral states

$\operatorname*{slim}_{N\to\infty}\Psi_{N}=|\alpha\rangle\otimes|\bar{\alpha}\rangle\,,$
(2.108)

where the symbol “$\operatorname*{slim}$” is used as a reminder that the
formula applies for the class of low energy states only. Also, strictly
speaking, the existence of the limit requires that the Bethe state $\Psi_{N}$
be properly normalized (for details see ref.[16]). As was discussed in [10]
(see also [34]), in a given sector with fixed $S^{z}$ the chiral state
$|\alpha\rangle$ can be labelled by two sets of non-negative integers

$1\leq n_{1}^{\pm}<n_{2}^{\pm}<\ldots<n_{M^{\pm}}^{\pm}\,.$ (2.109)

Similarly $|\bar{\alpha}\rangle$ is labelled by

$1\leq\bar{n}_{1}^{\pm}<\bar{n}_{2}^{\pm}<\ldots<\bar{n}_{\bar{M}^{\pm}}^{\pm}$
(2.110)

and the difference $\bar{M}^{+}-\bar{M}^{-}$ should coincide with
$M^{-}-M^{+}$. The latter can be identified with the so-called winding number:

${\tt w}=M^{-}-M^{+}=\bar{M}^{+}-\bar{M}^{-}\ .$ (2.111)

Formula (2.108) implies that the eigenvalues of the CFT Hamiltonian
$\hat{H}_{\rm CFT}$ split into two terms, $I_{1}+\bar{I}_{1}$ , where $I_{1}$
($\bar{I}_{1}$) stands for the contribution of the chiral state
$|\alpha\rangle$ ($|\bar{\alpha}\rangle$). In the case of the scaling limit of
the XXZ spin chain,

$I_{1}=\frac{p^{2}}{n+1}-\frac{1}{24}+{\tt
L}\,,\qquad\qquad\bar{I}_{1}=\frac{\bar{p}^{2}}{n+1}-\frac{1}{24}+\bar{\tt
L}\qquad\qquad(r=1)\ .$ (2.112)

Here

$p=\frac{S^{z}}{2}+\frac{n+1}{2}\ ({\tt k}+{\tt
w})\,,\qquad\qquad\bar{p}=\frac{S^{z}}{2}-\frac{n+1}{2}\ ({\tt k}+{\tt w})\ ,$
(2.113)

while the non-negative integers ${\tt L}$, $\bar{\tt L}$ are the conformal
levels. In terms of the sets (2.109), (2.110) they are given by

$\displaystyle{\tt L}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{M^{-}}\big{(}n_{j}^{-}-\tfrac{1}{2}\big{)}+\sum_{j=1}^{M^{+}}\big{(}n_{j}^{+}-\tfrac{1}{2}\big{)}-\sum_{j=1}^{|{\tt
w}|}\,(j-\tfrac{1}{2})$ (2.114)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ .$
$\displaystyle\bar{\tt L}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{\bar{M}^{-}}\big{(}\bar{n}_{j}^{-}-\tfrac{1}{2}\big{)}+\sum_{j=1}^{\bar{M}^{+}}\big{(}\bar{n}_{j}^{+}-\tfrac{1}{2}\big{)}-\sum_{j=1}^{|{\tt
w}|}\,(j-\tfrac{1}{2})$

As for the eigenvalue of the lattice translation operator for the low energy
state $\Psi_{N}$, one has

${\cal K}=(-1)^{{\tt w}}\,\exp\bigg{(}\frac{2\pi{\rm i}r}{N}\
\big{(}\,I_{1}-\bar{I}_{1}\,\big{)}\bigg{)}\qquad\qquad(r=1)\,.$ (2.115)

In Appendix D we briefly discuss the spectrum of the Baxter $Q$ \- operator
$\mathbb{A}_{+}(\zeta)$ (2.27) in the scaling limit for the ${\cal Z}_{r}$
invariant spin chain at the free fermion point $\gamma=\frac{\pi}{2}$. For a
low energy Bethe state in the sector with given $S^{z}$, the chiral states
$|\alpha\rangle$ and $|\bar{\alpha}\rangle$ appearing in the r.h.s. of (2.108)
can be labelled by the positive integers

$\displaystyle 1\leq
n_{1,a}^{\pm}<n_{2,a}^{\pm}<\ldots<n_{M^{\pm}_{a},a}^{\pm}$ (2.116)

and

$\displaystyle
1\leq\bar{n}_{1,a}^{\pm}<\bar{n}_{2,a}^{\pm}<\ldots<\bar{n}_{\bar{M}^{\pm}_{a},a}^{\pm}\
,$ (2.117)

respectively, where $a=1,2,\ldots,r$. The number of elements in the sets
$\\{n_{j,a}^{\pm}\\}$ and $\\{\bar{n}_{j,a}^{\pm}\\}$, i.e., $M^{\pm}_{a}$ and
$\bar{M}^{\pm}_{a}$, are subject to the constraint

$\sum_{a=1}^{r}\big{(}M^{-}_{a}-M^{+}_{a}\big{)}+\sum_{a=1}^{r}\big{(}\bar{M}^{-}_{a}-\bar{M}^{+}_{a}\big{)}=0\,.$
(2.118)

For $r>1$, the relation (2.111) is replaced by

$\frac{1}{r}\,\sum_{a=1}^{r}\big{(}M_{a}^{-}-M_{a}^{+}\big{)}=\frac{1}{r}\,\sum_{a=1}^{r}\big{(}\bar{M}_{a}^{+}-\bar{M}_{a}^{-}\big{)}={\tt
w}+\frac{{\tt m}}{r}\ ,$ (2.119)

which, together with the winding number ${\tt w}\in\mathbb{Z}$, involves

${\tt m}\ :\ \ \ \ \ {\tt m}=0,1,2,\ldots,r-1\ .$ (2.120)

Recall that $N$ is always assumed to be divisible by the odd integer $r$. For
this reason, in taking the limit (2.108) the size of the lattice should be
increased by $2r$ in order to remain in the sector with given $S^{z}$. As
such, one can introduce an extra quantum number

$\displaystyle{\tt s}=\frac{N}{2}-S^{z}\ \ \ \ \ \ \ ({\rm mod}\ r)\,,$
(2.121)

which takes the values

${\tt s}=0,1,2,\ldots,r-1\ .$ (2.122)

It is expected that, away from the free fermion point with
$\frac{\pi}{2}(1-\frac{1}{r})<\gamma<\frac{\pi}{2}(1+\frac{1}{r})$, the chiral
states can be characterized in the same way and, in particular, are labelled
by $S^{z},{\tt w}\in\mathbb{Z}$ and ${\tt m},{\tt s}=0,1,\ldots,r-1$.

In order to describe the eigenvalues of the conformal Hamiltonian (2.106), one
should distinguish between the following four cases:

$\begin{array}[]{ll}{\rm(i)}\ \ \,\frac{{\tt s}+S^{z}}{r}-S^{z}\ {\rm even\
and}\ {\tt k}\in(-\frac{1}{2},0)\,,&{\rm(ii)}\ \ \frac{{\tt
s}+S^{z}}{r}-S^{z}\ {\rm even\ and}\ {\tt k}\in(0,\frac{1}{2})\\\\[5.69046pt]
{\rm(iii)}\ \frac{{\tt s}+S^{z}}{r}-S^{z}\ {\rm odd\ and}\ {\tt
k}\in(-\frac{1}{2},0)\,,&{\rm(iv)}\ \frac{{\tt s}+S^{z}}{r}-S^{z}\ {\rm odd\
and}\ {\tt k}\in(0,\frac{1}{2})\end{array}$ (2.123)

(note that $\frac{{\tt s}+S^{z}}{r}-S^{z}$ is always an integer, as follows
from the definition (2.121) and that $r$ is odd). Also, introduce the integers
$\tilde{\tt m}$ and $\tilde{\tt s}$ as

${\tilde{\tt m}}=\begin{cases}|r-{\tt s}-{\tt m}|&\qquad{\rm for\ cases\ (i)\
\&\ (iv)}\\\ |{\tt s}-{\tt m}|&\qquad{\rm for\ cases\ (ii)\ \&\
(iii)}\end{cases}\ ,\qquad\qquad{\tilde{\tt s}}=\begin{cases}+{\tt
s}&\qquad{\rm for\ case\ (i)}\\\ -{\tt s}&\qquad{\rm for\ case\ (ii)}\\\
r-{\tt s}&\qquad{\rm for\ case\ (iii)}\\\ {\tt s}-r&\qquad{\rm for\ case\
(iv)}\end{cases}\ .$ (2.124)

Based on an analysis of the free fermion point and numerical work we come to

Conjecture: The eigenvalues of $\hat{H}_{\rm CFT}$ (2.106) are given by
$I_{1}+\bar{I}_{1}$ with

$I_{1}=\frac{p^{2}}{n+r}-\frac{r}{24}+\frac{{\tt m}\,(r-{\tt m})}{2r}+{\tt
L}\,,\qquad\qquad\bar{I}_{1}=\frac{\bar{p}^{2}}{n+r}-\frac{r}{24}+\frac{{\tilde{\tt
m}}\,(r-{\tilde{\tt m}})}{2r}+\bar{{\tt L}}\ ,$ (2.125)

in the cases (i) and (iv) and

$I_{1}=\frac{p^{2}}{n+r}-\frac{r}{24}+\frac{{{\tilde{\tt m}}}\,(r-{\tilde{\tt
m}})}{2r}+{\tt
L}\,,\qquad\qquad\bar{I}_{1}=\frac{\bar{p}^{2}}{n+r}-\frac{r}{24}+\frac{{{\tt
m}}\,(r-{{\tt m}})}{2r}+\bar{{\tt L}}\ ,$ (2.126)

in the cases (ii) and (iii). Here $p$ and $\bar{p}$ are expressed in terms of
the quantum numbers $S^{z}$, ${\tt w}$, ${\tt m}$ and ${\tt s}$ according to

$p+\bar{p}=S^{z}\,,\qquad\qquad\frac{p-\bar{p}}{n+r}={\tt k}+{\tt
w}+\frac{{\tt m}}{r}+\frac{\tilde{\tt s}}{2r}\ ,$ (2.127)

while ${\tt L}$, $\bar{\tt L}$ are non-negative integers. The eigenvalues of
the $r$ \- site lattice translation operator for the low energy states are
given by

${\cal K}=(-1)^{\frac{r-1}{2}(\frac{N}{2}-S^{z})+{\tt m}+r{\tt w}}\
\mbox{e}^{-\frac{{\rm i}\pi}{2}{\tilde{\tt s}}}\ \,\exp\bigg{(}\frac{2\pi{\rm
i}r}{N}\ \big{(}\,I_{1}-\bar{I}_{1}\,\big{)}\bigg{)}\ .$ (2.128)

.

####  2.5.2 Some comments on the low energy spectrum away from the ${\cal
Z}_{r}$ invariant point

Let ${\cal E}_{\rm vac}$ (vacuum energy) be the lowest energy of the spin
chain in the sector with given $S^{z}$ and $N/2-S^{z}$ is divisible by $r$,
i.e., the quantum number ${\tt s}$ (2.121) is equal to zero. Note that by the
ground state energy, ${\cal E}_{\rm GS}$, what is meant is the vacuum energy
for $S^{z}=0$. We studied the large - $N$ behaviour of ${\cal E}_{\rm vac}$,
where the inhomogeneities were taken to satisfy the conditions (2.33)-(2.35)
with $\mathfrak{a}_{2j+1}$ being treated as fixed, $N$ \- independent
parameters. It was found that the extensive part does not depend on the values
of the RG invariants and, of course, $S^{z}$ so that the specific bulk energy
$e_{\infty}$ is still given by eq. (2.105). However, the difference between
${\cal E}_{\rm vac}$ and ${\cal E}_{\rm vac}^{(0)}$ (the vacuum energy for the
${\cal Z}_{r}$ invariant spin chain) goes at large $N$ as

${\cal E}_{\rm vac}-{\cal E}^{(0)}_{\rm vac}=\frac{2\pi rv_{\rm
F}}{N}\,f(\mathfrak{a}_{1},\ldots,\mathfrak{a}_{r-2})+o(N^{-1})\ ,$ (2.129)

where $f$ is a certain polynomial in $\mathfrak{a}_{2j+1}$ whose coefficients
depend on $n$. Note that $f$ is independent of the value of $S^{z}$ as well as
the twist parameter ${\tt k}$. This motivated us to separate out the factor
$2\pi rv_{\rm F}$ in (2.129), where the Fermi velocity is defined by the same
formula as for the ${\cal Z}_{r}$ invariant case, i.e., (2.105). For $r=3$,

$\displaystyle f=-C\,\mathfrak{a}_{1}^{3}$ (2.130)

with $C$ being a positive $n$ \- dependent coefficient which vanishes as
$n^{-1}$ for $n\to\infty$ and diverges $\sim n^{-2}$ as $n\to 0$. Also, at the
free fermion point $C\big{|}_{n=3}=\frac{1}{2}$.

It turns out that the form of the leading behaviour (2.129) is not universal
in that it depends on the way in which one performs the scaling limit. For
example, as was discussed, the limit can be organized differently such that
the inhomogeneities obey (2.90), (2.94) and the spin chain Hamiltonian
possesses both ${\cal CP}$ and ${\cal T}$ symmetries for any number of lattice
sites. If $r=3$, there would be a single real valued RG invariant
$\mathfrak{b}_{1}$, which is defined by eqs. (2.96) - (2.98). Then, we
observed that

$\displaystyle{\cal E}_{\rm vac}-{\cal E}^{(0)}_{\rm
vac}=-\frac{C_{1}}{N^{\frac{1}{3}}}\,\mathfrak{b}^{2}_{1}-\frac{2\pi rv_{{\rm
F}}}{N}\,C_{2}\,{\mathfrak{b}}^{3}_{1}+o(N^{-1})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ (r=3)\ ,$ (2.131)

where $C_{1}$, $C_{2}$ are some positive constants depending on $n$.

Regardless of the scheme in which the scaling limit is taken, the operator
$N(\mathbb{H}-{\cal E}_{\rm GS})$, when restricted to the class of low energy
states, is expected to possess the $N\to\infty$ limit. The low energy states
can be chosen to be the Bethe states, which are specified by solutions to the
Bethe Ansatz equations. It seems to be possible to unambiguously and
continuously deform the corresponding Bethe roots away from the ${\cal Z}_{r}$
invariant case, at least in some domain of the values of the RG invariants. A
preliminary numerical study shows that the spectrum of the operator

$\displaystyle\hat{H}\equiv\frac{1}{2}\,(n+r)\,{\tt
k}^{2}-\frac{r}{12}+\frac{1}{2\pi rv_{\rm F}}\
\operatorname*{slim}_{N\to\infty}\,N\big{(}{\mathbb{H}}-{\cal E}_{\rm
GS}\big{)}$ (2.132)

depends on the RG invariants. For $n=r$, i.e., at the free fermion point, this
can also be confirmed via a straightforward extension of the analysis
presented in Appendix D. Thus, unlike the case $A=r-1$ previously discussed in
ref.[21] (see also sec. 3.1 below), the conformal structure is affected by our
scaling procedure. Its detailed investigation is an interesting problem that
lies beyond the scope of this paper.

##  3 Basic constraints on ODEs for spin chain with anisotropy $0<\gamma<\pi$

The quantum Wronskian relation(s), obeyed by the $Q$ \- operators, first
appeared in the series of works [4, 5, 6] and has since then become a central
notion in integrability in Quantum Field Theory and Statistical Mechanics. In
the case of the inhomogeneous six-vertex model this would be the relation
(2.28) satisfied by $\mathbb{A}_{\pm}(\zeta)$ and it holds true for any values
of the anisotropy parameter $q$ and inhomogeneities $\eta_{J}$. It is expected
that, with a properly defined scaling limit, the eigenvalues of
$\mathbb{A}_{\pm}(\zeta)$ corresponding to the low energy states would become
the functions $D_{\pm}(E)$, which obey the functional relation inherited from
the quantum Wronskian relation. For the ground state of the spin chain this
would be (2.20) with $p=\frac{n+r}{2}\,{\tt k}$. In developing the ODE/IQFT
correspondence the main step is to identify $D_{\pm}(E)$ with the spectral
determinants of a certain 2nd order ODE. Unfortunately, we do not know of a
systematic way to do this at the moment. While the previous section was
devoted to the case with $r$ odd and $q=\mbox{e}^{{\rm i}\gamma}$, where
$\gamma\in\big{(}\frac{\pi}{2}(1-\frac{1}{r}),\frac{\pi}{2}(1+\frac{1}{r})\big{)}$,
it contains the essential hints for obtaining the differential equation
appearing in the scaling limit of the spin chain for the other domains

$\gamma\in\big{(}\tfrac{\pi}{r}A,\tfrac{\pi}{r}(A+1)\big{)}$ (3.1)

labelled by the integer $A=0,1,\ldots,r-1$. Among them is that the class of
possible ODEs is greatly reduced by the requirement that the differential
equation be invariant w.r.t. the transformation

$\hat{\Omega}\ :\ y\mapsto y+\frac{2\pi{\rm i}}{n+r}\,,\qquad\qquad E\mapsto
q^{-2}\,E\,,$ (3.2)

which is used to derive the quantum Wronskian relation for the spectral
determinants. Our proposal is that the scaling limit can be defined in such a
way that the ground state of the spin chain is described in terms of the
differential equation (1.23), i.e.,

$\displaystyle\Big{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-\delta U(y)\,\Big{]}\,\psi=0$ (3.3)

with

$\displaystyle\delta U(y)=\sum_{(\mu,j)\in\bm{\Xi}_{r,A}}c_{\mu,j}\,E^{\mu}\
\mbox{e}^{\big{(}(A\mu-rj)\,\frac{n+r}{r}+\mu\big{)}y}$ (3.4)

and its barred counterpart. It is easy to check that the ODE is indeed
invariant w.r.t. the transformation $\hat{\Omega}$ with
$q=\mbox{e}^{\frac{{\rm i}\pi}{r}A+\frac{{\rm i}\pi}{n+r}}$. Formally this
holds true for integer $j$ and arbitrary values of $\mu$. However, in order
for the spectral determinant to be an entire function of $E$, $\mu$ should be
a non-negative integer. To specify the ODE completely one needs to describe
the set $\bm{\Xi}_{r,A}$, which occurs in the sum in (3.4). For the case $r$
odd and $A=\frac{r-1}{2}$, this has already been carried out in the previous
section,

$\bm{\Xi}_{r,\frac{r-1}{2}}=\big{\\{}(\mu,j)\,:\ \ \mu=2j+1\ \ \&\ \
j=0,\ldots,\tfrac{r-3}{2}\big{\\}}\,\qquad\qquad(r\ {\rm odd})\,.$ (3.5)

Then the differential equation (3.3) becomes (2.18) with $c_{2j+1}\equiv
c_{2j+1,j}$

###  3.1 The case $A=r-1$

The ODE/IQFT correspondence in the regime $\pi(1-\frac{1}{r})<\gamma<\pi$ was
proposed in ref.[21]. Here we give a brief summary of the results of that work
relevant to our current discussion.

For $A=r-1$ the integers $\mu$ and $j$ are specified as

$\displaystyle\mu=j+1\ \ \ \ \ \ \ \quad{\rm with}\quad\ \ \ \ \ \
j=0,\ldots,A-1\,.$ (3.6)

Then the differential equation (3.3) takes the form

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}+(-E)^{r}\,\mbox{e}^{ry}-\sum_{\mu=1}^{r-1}c_{\mu}\,E^{\mu}\
\mbox{e}^{ry+\frac{n}{r}(r-\mu)y}\,\bigg{]}\,\psi=0\,.$ (3.7)

Performing the change of variables

$\displaystyle z=-E^{-1}\,\mbox{e}^{\frac{n}{r}y}\ ,\ \ \ \
\Psi=\mbox{e}^{\frac{n}{2r}y}\,\psi$ (3.8)

one arrives at

$\displaystyle\Big{[}-\partial^{2}_{z}+z^{-2}\,\Big{(}\big{(}\tfrac{rp}{n}\big{)}^{2}-\tfrac{1}{4}\Big{)}+\kappa^{2}z^{-2+\xi
r}\,P_{r}(z)\Big{]}\,\Psi=0\,.$ (3.9)

Here

$\displaystyle\kappa^{2}=(-1)^{r}\,\Big{(}\frac{rp}{n}\Big{)}^{2}\,E^{\frac{r(n+r)}{n}}\
,\qquad\qquad\xi=\frac{r}{n}\,,$ (3.10)

while $P_{r}(z)$ is a polynomial of degree $r$ given by

$\displaystyle
P_{r}(z)=z^{r}-\sum_{\mu=1}^{r-1}(-1)^{\mu}\,c_{\mu}\,z^{r-\mu}+1\ .$ (3.11)

If one assumes that all the coefficients $c_{\mu}$ are sufficiently small,
then the roots of this polynomial would be close to $\mbox{e}^{\frac{{\rm
i}\pi}{r}(2a+1)}$ and hence simple. This way, $P_{r}(z)$ can be written as

$\displaystyle P_{r}(z)=\prod_{a=1}^{r}(z-z_{a})\qquad\qquad{\rm
with}\qquad\qquad\prod_{a=1}^{r}(-z_{a})=1\,,$ (3.12)

where $z_{a}$ are distinct.

It was proposed in section 12 of ref.[21] that the differential equation (3.9)
(3.12) describes the scaling limit of the Bethe roots for the ground state of
the spin chain with $q=-\mbox{e}^{\frac{{\rm i}\pi}{r}(\beta^{2}-1)}$, where
$0<\beta^{2}\equiv\frac{r}{n+r}<1$. In taking the limit, the inhomogeneities
are assumed to be fixed and sufficiently close to their values for the ${\cal
Z}_{r}$ invariant case. Also, without loss of generality, the normalization
condition was imposed as

$\prod_{\ell=1}^{r}\eta_{\ell}=1\,.$ (3.13)

In our notation, this prescription for the scaling limit corresponds to
setting the exponent $d_{s}$ in the definition of the RG invariants
$\mathfrak{a}_{s}$ (1.22) to be zero, i.e.,

$\displaystyle{\mathfrak{a}}_{s}=\frac{1}{sr}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-s}\qquad\qquad(s=1,2,\ldots,r-1)\,.$ (3.14)

The relation between the set $\\{\mathfrak{a}_{s}\\}_{s=1}^{r-1}$ and the
coefficients $\\{c_{\mu}\\}_{\mu=1}^{r-1}$ of the differential equation (3.7)
is rather complicated. It is described explicitly for the case $r=2$ in
section 12.3 in [21].

For $A=r-1$ the full class of ODEs which appear in the scaling limit of the
low energy excited states of the spin chain is known [21]. It will not be
presented here. We just mention that though the integrable structure
underlying the CFT depends on the RG invariants, they have no effect on the
conformal structure. In particular, the eigenvalues and degeneracies of the
Hamiltonian $\hat{H}$ defined through (2.132) does not depend on
$\mathfrak{a}_{s}$ (3.14) provided that their absolute values are sufficiently
small.

###  3.2 General restrictions on integers $\mu$ and $j$

As was mentioned in the Introduction, our work is devoted to the study of the
scaling limit of the spin chain with “softly” broken ${\cal Z}_{r}$ symmetry.
At the level of the differential equation (3.3), this can be formulated as the
requirement that the term $\delta U$ (3.4) may be treated as a perturbation.
In view of this, we impose the following constraints:

  1. (a)

The function

$\displaystyle\delta U(y)=\sum_{(\mu,j)}c_{\mu,j}\,E^{\mu}\
\mbox{e}^{\big{(}(A\mu-rj)\,\frac{n+r}{r}+\mu\big{)}y}$ (3.15)

vanishes as $y\to-\infty$:

$\displaystyle\lim_{y\to-\infty}\delta U(y)=0\ \ \ \ \ {\rm for}\ \ \ \forall
n>0\ .$ (3.16)

Then one can still introduce the Jost solutions of the ODE, which are involved
in the derivation of the quantum Wronskian relation (2.20), through the
asymptotic conditions

$\displaystyle\psi_{\pm p}(y)\to\mbox{e}^{\pm py}\ \ \ \ \ \ \qquad{\rm as}\ \
\ \qquad y\to-\infty\qquad\qquad\big{(}\Re e(\pm p)\geq 0\big{)}\,.$ (3.17)

It is easy to see that (3.16) implies that both the integers $\mu$ and $j$ can
not be simultaneously zero as well as the inequalities

$\displaystyle\mu\geq{\rm max}\Big{(}0,\frac{rj}{A}\Big{)}\qquad{\rm
for}\qquad A\geq 1$ (3.18a) $\displaystyle\mu\geq 0\ \&\ j\leq 0\qquad\ \ \
{\rm for}\qquad A=0\,.$ (3.18b) Here we also took into account that $\mu$ must
be non-negative, which is required in order that the spectral determinant
admits a power series expansion at $E=0$.

  2. (b)

The other constraint is that $\delta U(y)$ should grow slower than the term
$\propto\mbox{e}^{(n+r)y}$ in eq. (3.3) as $y\to+\infty$, i.e.,

$\displaystyle\lim_{y\to+\infty}\mbox{e}^{-(n+r)y}\,\delta U(y)=0\ \ \ \ \
{\rm for}\ \ \ \forall n>0\ .$ (3.19)

This implies

$\displaystyle\mu<\frac{r}{A+1}\ (j+1)\ \ \ \ \ \ \ \ \ \ \ \ \ \qquad(0\leq
A<r-1)\ .$ (3.20)

It follows from (3.19) that the zeroes of the spectral determinants would
accumulate along the same rays as for the ${\cal Z}_{r}$ invariant case. For
$A=r-1$ it was discussed in the previous subsection that $\mu=j+1$ and the
inequality (3.20) is saturated. In turn, the Stokes lines are considerably
more cumbersome to describe (see, e.g., Appendix B in ref.[21]). This makes
the case $A=r-1$ somewhat special.

The conditions (3.18) and (3.20) have some immediate consequences. In
particular, one has for the integer $j$ that

$0\leq j\leq A-1\qquad\qquad(A=1,\ldots,r-2)\,.$ (3.21)

Since we exclude the case when both $\mu$ and $j$ are zero, it follows from
(3.18a) that $\mu$ is a positive integer and, furthermore,

$1\leq\mu<\frac{A}{A+1}\ r\,,\qquad\qquad\qquad(A=1,\ldots,r-2)\,.$ (3.22)

As for $A=0$, the possible values of the pair $(\mu,j)$ are described as

$\mu=1,2,\ldots,r-1\ \&\ j=0\qquad\qquad\qquad(A=0)\,.$ (3.23)

If $\delta U(y)$ is given by (3.15), where the summation is taken over the
pairs $(\mu,j)$ satisfying the above conditions, the spectral determinants for
the ODE (3.3) can be introduced along the same lines as was done in secs. 2.2
and 2.3.666Strictly speaking formula (2.23) used to specify the subdominant at
$y\to+\infty$ solution $\chi$, entering into the definition (2.24) of the
spectral determinant, is not valid in general. Nevertheless, $\chi$ may still
be unambiguously introduced via its large - $y$ asymptotic obtained within the
WKB approximation. The proof that they satisfy the quantum Wronskian relation
(2.20), which is based on the invariance of the differential equation w.r.t.
the transformation (3.2), remains unchanged. Further, $D_{\pm}(E)$ are entire
functions of $E$.

##  4 ODE with single coefficient $c_{\mu,j}$ for $A\mu>rj$ and
$A=1,\ldots,r-2$

We now come up against the problem of giving a precise description of the
scaling limit of the spin chain that would result in the differential equation
belonging to the class (3.3) and (3.4). This involves specifying the exponents
$d_{s}$ entering into the definition of the RG invariants $\mathfrak{a}_{s}$
(1.22) as well as relating their values to the coefficients $c_{\mu,j}$ of the
ODE. However, when carried out in full generality, one runs into, besides the
technical complexity, also conceptual issues. For this reason, we first focus
our attention on formulating the scaling limit such that the corresponding
differential equation is given by

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-c\,E^{\mu}\ \mbox{e}^{\big{(}(A\mu-
rj)\,\frac{n+r}{r}+\mu\big{)}y}\,\bigg{]}\,\psi=0\qquad\qquad(c\equiv
c_{\mu,j})\,,$ (4.1)

where the pair of integers $(\mu,j)$ obey the conditions

$\frac{rj}{A}<\mu<\frac{r}{A+1}\ (j+1)\ \&\ j\geq 0\qquad\qquad(1\leq
A<r-1)\,.$ (4.2)

Notice that we have excluded the case $A\mu=rj$ which, as will be discussed
below, requires a separate treatment.

###  4.1 Bohr-Sommerfeld quantization condition

The experience gained in sec. 2.3 suggests to start with the analysis of the
large - $E$ asymptotic behaviour of the spectral determinant. A computation,
based on the WKB approximation, leads to the asymptotic formula

$\displaystyle\log D_{+}$ $\displaystyle\asymp$
$\displaystyle\frac{N_{0}\mbox{e}^{\theta}}{\cos(\frac{\pi
r}{2n})}-\sum_{1\leq k\leq\frac{r}{2M}}\mbox{e}^{(r-2kM)\frac{\theta}{r}}\
\mbox{e}^{\frac{{\rm i}\pi}{r}(2a-1-A)k\mu}\
\big{(}(-1)^{(r-1)\mu}\,c\big{)}^{k}$ $\displaystyle\times$
$\displaystyle\frac{\Gamma\big{(}\frac{k\mu}{r}-\frac{n+r}{2nr}\
(r-2kM)\big{)}\Gamma\big{(}k-\frac{1}{2}-\frac{k\mu}{r}+\frac{n+r}{2nr}\
(r-2kM)\big{)}}{2n\sqrt{\pi}\,k!}-\frac{2np\,\theta}{r(n+r)}+\log
C_{p}+o(1)\,.$

Here $N_{0}$ and $C_{p}$ are given by (1.17) and (2.51), respectively, while

$\displaystyle M=(j+1)\,r-(A+1)\,\mu>0$ (4.4)

is a positive integer (which follows from (4.2)). Again, $\theta$
parameterizes $E$ in the wedge labelled by the integer $a$ as

$\displaystyle E=(-1)^{r-1}\ \mbox{e}^{\frac{{\rm i}\pi}{r}(2a-1-A)}\
\mbox{e}^{\frac{2n\theta}{r(n+r)}}$ (4.5)

and the asymptotic formula is valid for $\Re e(\theta)\to+\infty$ with
$\big{|}\Im m(\theta)\big{|}<\frac{\pi(n+r)}{2n}$.

One can also derive the Bohr-Sommerfeld quantization condition, from which the
location of the zeroes $E_{m}^{(a)}$ of the spectral determinant $D_{+}(E)$
for $m\gg 1$ may be deduced. Writing the zeroes in terms of $\theta_{m}^{(a)}$
according to (2.66), one finds that

$\displaystyle X=\mbox{e}^{-\frac{{\rm i}\pi\mu}{2rM}(2a-A)}\
\mbox{e}^{\frac{1}{r}\theta_{m}^{(a)}}$ (4.6)

with $a=1,2,\ldots,r$ and $m\gg 1$ solves the equation

$\displaystyle X^{r}+\sum_{k=1}^{[\frac{r}{2M}]}{G}_{2kM}\
X^{r-2kM}+o(1)=Y^{r}\ .$ (4.7)

Here $[\ldots]$ stands for the integer part and the notation

$\displaystyle Y=\mbox{e}^{-\frac{{\rm i}\pi\mu}{2rM}(2a-A)}\
\bigg{(}\frac{\pi}{N_{0}}\
\Big{(}m-\frac{1}{2}+\frac{p}{r}\,\Big{)}\bigg{)}^{\frac{1}{r}}\ ,\ \ \ \ \
{G}_{2kM}={\mathfrak{g}}_{2k}\ \Big{(}(-1)^{(r-1)\mu}\ c\Big{)}^{k}$ (4.8)

as well as

$\displaystyle{\mathfrak{g}}_{2k}=\frac{\Gamma(\frac{3}{2}+\frac{r}{2n})}{\Gamma(\frac{r}{2n})}\
\frac{\Gamma\big{(}\,\frac{k}{r}\,L+\frac{1}{2n}\
(r-2kM)\big{)}}{k!\,\Gamma\big{(}\frac{3}{2}-k+\frac{k}{r}\,L+\frac{1}{2n}\
(r-2kM)\big{)}}$ (4.9)

is being used. Also, together with the positive integer $M$ (4.4), we
introduce

$L=A\mu-rj>0\,,$ (4.10)

which is a positive integer due to the condition (4.2). Notice that the
argument of the $\Gamma$ \- function occurring in the numerator of the second
fraction in $\mathfrak{g}_{2k}$ is always positive since, as follows from the
upper limit in the sum in (4.7), $k\leq\frac{r}{2M}$. It vanishes when this
bound is saturated and $L=0$, which is only possible if both $A$ and $r$ are
divisible by $2K$ and

$j=\frac{A}{2K}\,(2K-1)\,,\ \ \qquad\qquad\ \
\mu=\frac{r}{2K}\,(2K-1)\,\qquad\qquad(K=1,2,\ldots)\ .$ (4.11)

Then, the last term in the sum in (4.7) would be proportional to $\log(X/{\rm
const})$. This, among other things, makes the case $L=0$ somewhat special and
is one of the reasons why it will be excluded from the discussion of this
section.

If the positive integer $M$ entering into the quantization condition (4.7)
exceeds $r/2$, the sum vanishes and hence, as $m\to+\infty$,

$\displaystyle E_{m}^{(a)}=(-1)^{r-1}\,\mbox{e}^{\frac{{\rm
i}\pi}{r}\,(2a-A)}\ \bigg{(}\frac{\pi}{N_{0}}\
\Big{(}m-\frac{1}{2}+\frac{p}{r}\,\Big{)}\bigg{)}^{\frac{2n}{r(n+r)}}\
\Big{(}1+o(m^{-1})\Big{)}\ \ \ \ \ \ \ \ \ \ (M>r/2)\ .$ (4.12)

If $0<M\leq r/2$, i.e.,

$\displaystyle\frac{r}{A+1}\ \Big{(}j+\
\frac{1}{2}\Big{)}\leq\mu<\frac{r}{A+1}\ (j+1)\ ,$ (4.13)

the sum in (4.7) is not trivial and the leading large - $m$ asymptotic formula
for $E_{m}^{(a)}$ reads as

$\displaystyle E_{m}^{(a)}=(-1)^{r-1}\,\mbox{e}^{\frac{{\rm
i}\pi}{r}\,(2a-A)}\ \bigg{(}\frac{\pi}{N_{0}}\
\Big{(}m-\frac{1}{2}+\frac{p}{r}\,\Big{)}\bigg{)}^{\frac{2n}{r(n+r)}}\
\bigg{(}1+O\Big{(}m^{-\frac{2M}{r}}\Big{)}\bigg{)}\ \ \ \ \ \ \ (0<M\leq
r/2)\,.$ (4.14)

In both cases the zeroes accumulate along the same Stokes lines as in the
${\cal Z}_{r}$ invariant spin chain. However, the subleading term in (4.14),
which evaluates to a complex number, decays considerably slower than
$o(m^{-1})$ from eq. (4.12) so that the zeroes approach to the rays less
rapidly when the value of the integer $\mu$ is restricted to the interval
(4.13). This, in turn, impacts the convergence of the series

$\displaystyle J_{s}=\frac{1}{s}\
\sum_{m=1}^{\infty}\sum_{a=1}^{r}\big{(}E_{m}^{(a)}\big{)}^{-s}\,,$ (4.15)

which would require, in general, the introduction of counterterms. Similar to
what was discussed in sec. 2.3, one can make use of this to try to deduce the
relation between the RG invariants (1.22) and the coefficient $c$ in (4.1).

###  4.2 Specification of the scaling limit

Let’s first consider the case when the pair $(\mu,j)$ satisfies the inequality

$\mu\geq\frac{r}{A+1}\ \Big{(}j+\ \frac{1}{2}\Big{)}$ (4.16)

in addition to (4.2). Numerical work shows that it is possible to organize the
scaling limit in such a way that the scaled Bethe roots for the ground state
of the spin chain are described by the ODE of the form (4.1). The exponent
$d_{s}$, which appears in the definition of the RG invariants,

$\displaystyle{\mathfrak{a}}_{s}=\frac{1}{s}\
\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{d_{s}}\ \frac{1}{r}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-s}\,,$ (4.17)

should be taken to be

$d_{s}=\frac{2M}{r}\,i_{s}\,,$ (4.18)

where $i_{s}$ denotes the smallest positive integer such that

$s=\mu\,i_{s}\ \ \ ({\rm mod}\ r)$ (4.19)

and, as before, $M=(j+1)\,r-(A+1)\,\mu$. If $s$ is not divisible by $\mu$
modulo $r$, i.e., the integer $i_{s}$ does not exist, the corresponding RG
invariant can be set to zero:

$\mathfrak{a}_{s}=0\qquad{\rm if}\qquad s\nmid\mu\ ({\rm mod}\ r)\,.$ (4.20)

In addition, we assume the normalization condition for the inhomogeneities

$\frac{1}{r}\,\sum_{\ell=1}^{r}(\eta_{\ell})^{-r}=(-1)^{r-1}\,.$ (4.21)

The inequality (4.16) guarantees that there are RG invariants
$\mathfrak{a}_{s}$ for which $0<d_{s}\leq 1$. In this case one can use the
trick discussed in sec. 2.3 to deduce the relation

$\displaystyle{\mathfrak{a}}_{s}=(-1)^{s-(j+1)i_{s}+{\cal N}_{s}}\ Z_{s}\
c^{i_{s}}\ .$ (4.22)

Here the integer ${\cal N}_{s}$ is defined by the condition

$\displaystyle s=\mu\,i_{s}+r\,{\cal N}_{s}\ .$ (4.23)

The coefficients $Z_{s}$ are built from $\mathfrak{g}_{2k}$ (4.9) as

$\displaystyle Z_{s}=\frac{1}{s}\sum_{m=1}^{i_{s}}(-1)^{m}\
\frac{\Gamma(s\nu+m)}{m!\Gamma(s\nu)}\sum_{k_{1},\ldots,k_{m}\geq 1\atop
k_{1}+k_{2}+\ldots+k_{m}={i_{s}}}{\mathfrak{R}}_{k_{1}}\ldots\mathfrak{R}_{k_{m}}$
(4.24)

with

$\mathfrak{R}_{k}=\frac{1}{r}\sum_{\alpha_{1}+2\alpha_{2}+\ldots+k\alpha_{k}=k\atop\alpha_{1},\ldots,\,\alpha_{k}\geq
0}\frac{(-1)^{\alpha_{1}+\ldots+\alpha_{k}}}{\alpha_{1}!\alpha_{2}!\ldots\alpha_{k}!}\
\frac{\Gamma(\alpha_{1}+\ldots+\alpha_{k}-\frac{2kM-1}{r})}{\Gamma(1-\frac{2kM-1}{r})}\
\
{\mathfrak{g}}_{2}^{\alpha_{1}}{\mathfrak{g}}_{4}^{\alpha_{2}}\dots\mathfrak{g}_{2k}^{\alpha_{k}}$
(4.25)

and

$\displaystyle\nu=\frac{2n}{n+r}\,.$ (4.26)

It needs to be emphasized that the argument used to obtain (4.22) can not be
literally applied to the case when the exponent $d_{s}$ for the RG invariant
$\mathfrak{a}_{s}$ exceeds one. Nevertheless, we found that under certain
conditions which will be clarified in the examples below, the scaling limit of
the Bethe roots is described by the differential equation (4.1) if the RG
invariants are specified as in (4.20), (4.22) - (4.26) for any $s$. Moreover,
numerical work shows that the values of some of the RG invariants
$\mathfrak{a}_{s}$, where $s$ is divisible by $\mu$ modulo $r$ and $d_{s}>1$
have no effect on the scaling limit and can be set to zero. In particular, for
$r\leq 10$ this was always observed if the corresponding exponent $d_{s}\geq
2$. Unfortunately, we do not have at hand a full description of the minimal
set of non-trivial RG invariants that works for generic values of $r$ and
$A=1,\ldots,r-2$.777Here and below, by “non-trivial RG invariants” we mean
those $\mathfrak{a}_{s}$, whose values affect the scaling limit of the Bethe
roots. For this reason, we list them in the tables in Appendix F, along with
the corresponding exponents $d_{s}$ for all the cases with $r=3,4,\ldots,10$,
which were deduced from a numerical study.

####  4.2.1 Example: $r=5$, $A=1$ and $(\mu,j)=(2,0)$

In this case we found that there are three non-trivial RG invariants
$\mathfrak{a}_{2}$, $\mathfrak{a}_{4}$ and $\mathfrak{a}_{1}$, whose
corresponding exponents are given by $d_{2}=\tfrac{2}{5}$,
$d_{4}=\tfrac{4}{5}$ and $d_{1}=\tfrac{6}{5}$, respectively. A specialization
of formulae (4.22) - (4.26) yields

$\displaystyle{\mathfrak{a}}_{2}$ $\displaystyle=$
$\displaystyle-\tfrac{\nu}{5}\,\mathfrak{g}_{2}\ c$
$\displaystyle{\mathfrak{a}}_{4}$ $\displaystyle=$
$\displaystyle\tfrac{\nu}{5}\,\big{(}\mathfrak{g}_{4}+\tfrac{1}{10}\,(4\nu-1)\
\mathfrak{g}_{2}^{2}\big{)}\ c^{2}$ $\displaystyle{\mathfrak{a}}_{1}$
$\displaystyle=$
$\displaystyle-\tfrac{\nu}{5}\,\big{(}\mathfrak{g}_{6}+\tfrac{1}{5}\,(\nu+1)\,\mathfrak{g}_{2}\mathfrak{g}_{4}+\tfrac{1}{150}\,(\nu+1)(\nu-4)\,\mathfrak{g}_{2}^{3}\big{)}\
c^{3}\,,$ (4.27b)

where $\nu=\frac{2n}{n+5}$ and the functions $\mathfrak{g}_{2k}$ are defined
in eq.(4.9) with $M=1$ and $L=2$. Then, according to the above discussion, the
ODE appearing in the scaling limit is expected to be

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+5)y}+E^{5}\
\mbox{e}^{5y}-c\,E^{2}\,\mbox{e}^{\frac{2}{5}(n+10)\,y}\,\bigg{]}\,\psi=0\ ,$
(4.28)

where recall that $p=\frac{1}{2}\,(n+5)\,{\tt k}$, see (1.19). One may observe
that both $\mathfrak{g}_{2}$ and $\mathfrak{g}_{4}$ are non-singular functions
of $n$ for any $n>0$. In contrast $\mathfrak{g}_{6}$, which reads explicitly
as

$\mathfrak{g}_{6}=\frac{\Gamma(\frac{3}{2}+\frac{5}{2n})}{\Gamma(\frac{5}{2n})}\,\frac{\Gamma(\frac{6}{5}-\frac{1}{2n})}{6\,\Gamma(-\frac{3}{10}-\frac{1}{2n})}\,,$
(4.29)

possesses simple poles when $n=\frac{5}{2(1+5k)}$ with $k=1,2,3,\ldots\,$,
accumulating at $n=0$. As a result, the proportionality coefficient
$Z_{1}^{(\mu=2)}=-\mathfrak{a}_{1}/c^{3}$ also contains poles at the same
values of $n$, see eq. (4.27b). Moreover, $Z_{1}^{(\mu=2)}$ has zeroes with
the largest one located at $n=\frac{5}{2}$. From the numerical work, we were
able to confirm that the scaling limit of the Bethe roots for the ground state
is governed by the ODE (4.28) for $n>\frac{5}{2}$. As for the domain
$0<n\leq\frac{5}{2}$, it turns out to be difficult to explore for reasons to
be discussed below.

####  4.2.2 Example: $r=5$, $A=1$ and $(\mu,j)=(1,0)$

So far, we have taken $\mu$ and $j$ to obey the inequality (4.16) in addition
to (4.2). For $\mu<\frac{r}{A+1}\ (j+\ \frac{1}{2})$ or, equivalently, the
integer $M$ (4.4) exceeds $\frac{r}{2}$, it follows from formula (4.18) that
all the exponents corresponding to the non-trivial RG invariants are greater
than one. The numerical study shows that the relations obtained assuming $\mu$
satisfies both $\frac{r}{A+1}\ (j+\ \frac{1}{2})\leq\mu$ and
$\frac{rj}{A}<\mu<\frac{r}{A+1}\ (j+1)$ remain valid if the former condition
is dropped. This can be illustrated on the example $r=5$, $A=1$ and
$(\mu,j)=(1,0)$. Then, there is only one non-trivial RG invariant
$\mathfrak{a}_{1}$ with corresponding exponent $d_{1}=\frac{6}{5}>1$. Formula
(4.22), specialized to the case at hand, reads as

$\mathfrak{a}_{1}=Z_{1}^{(\mu=1)}\,c\qquad\qquad{\rm with}\qquad\qquad
Z_{1}^{(\mu=1)}=\frac{\Gamma(\frac{1}{2}+\frac{5}{2n})}{5\,\Gamma(\frac{5}{2n})}\,\frac{\Gamma(\frac{1}{5}-\frac{1}{2n})}{\Gamma(\frac{7}{10}-\frac{1}{2n})}\
.$ (4.30)

Notice that the factor $Z_{1}^{(\mu=1)}$ remains finite and non-vanishing as
$n>\frac{5}{2}$. In this domain, it was numerically verified that the
corresponding differential equation is given by

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+5)y}+E^{5}\
\mbox{e}^{5y}-c\,E\,\mbox{e}^{\frac{1}{5}(n+10)\,y}\,\bigg{]}\,\psi=0\ .$
(4.31)

####  4.2.3 Example: combining the cases $\mu=1$ and $\mu=2$ for $r=5$, $A=1$,
$j=0$

We found it rather surprising that in the prescription for taking the scaling
limit, relations (4.22)-(4.26) hold true for the RG invariants
$\mathfrak{a}_{s}$ with $d_{s}>1$ at least for sufficiently large $n$. For
this reason, we performed an additional numerical study with a slightly
different setup. The exponents $d_{s}$ for all the non-trivial RG invariants
(4.17) were taken as in (4.18). Furthermore, the value of the RG invariants
$\mathfrak{a}_{s}$ with $0<d_{s}\leq 1$ were fixed according to (4.22), while
those with $d_{s}>1$ were kept as free parameters. It was observed that for
sufficiently large $n$ the scaling limit of the Bethe roots for the ground
state is described in terms of the differential equation of the form

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-c\,E^{\mu}\ \mbox{e}^{\big{(}(A\mu-
rj)\,\frac{n+r}{r}+\mu\big{)}y}-\delta\tilde{U}(y)\,\bigg{]}\,\psi=0\,,$
(4.32)

where

$\delta\tilde{U}(y)=\sum_{s\in{\Sigma}}\ \tilde{c}_{s}\,E^{s}\
\mbox{e}^{\big{(}(As-rj)\,\frac{n+r}{r}+s\big{)}y}$ (4.33)

and

$\Sigma=\Big{\\{}s:\ s=\mu\,i\ ({\rm mod}\,r)\ \ {\rm with}\
i=\big{[}\tfrac{r}{2M}\big{]}+1,\ldots,r-1\ \ \ \&\ \ \ \tfrac{rj}{A}\leq
s<\tfrac{r}{A+1}\,(j+\tfrac{1}{2})\Big{\\}}\ .$ (4.34)

The appearance of such an ODE may be expected since it belongs to the class
(3.3), (3.4), while the presence of the term $\delta\tilde{U}$ does not affect
the Bohr-Sommerfeld quantization condition (4.7), which was derived for the
case with $\delta\tilde{U}=0$. The coefficients $\tilde{c}_{s}$ depend on $c$
as well as the values of $\mathfrak{a}_{s}$ with $d_{s}>1$. They can be
obtained numerically, for instance, from the study of the first few sum rules,
along the lines of sec. 2.2.

As an example, consider again the case with $r=5$, $A=1$, $(\mu,j)=(2,0)$. The
RG invariants $\mathfrak{a}_{2}$ and $\mathfrak{a}_{4}$, with corresponding
exponents $d_{2}=\frac{2}{5}$ and $d_{4}=\frac{4}{5}$, are taken as in (4.27),
while the value of

$\displaystyle{\mathfrak{a}}_{1}=\bigg{(}\frac{N}{5N_{0}}\bigg{)}^{\frac{6}{5}}\
\frac{1}{5}\ \sum_{\ell=1}^{5}(\eta_{\ell})^{-1}$ (4.35)

is kept free now. Then, the scaling limit of the Bethe roots for the ground
state would be described by the ODE

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+5)y}+E^{5}\
\mbox{e}^{5y}-c\,E^{2}\,\mbox{e}^{\frac{2}{5}(n+10)\,y}-\tilde{c}_{1}\,E\,\mbox{e}^{\frac{1}{5}(n+10)y}\,\bigg{]}\,\psi=0\
.$ (4.36)

Among other things, the differential equation implies

$\bigg{(}\frac{N}{5N_{0}}\bigg{)}^{-\frac{2n}{5(n+5)}}\,h_{1}^{(N,{\rm
reg})}=\tilde{c}_{1}\,\frac{(n+5)^{\frac{2}{n+5}-\frac{8}{5}}}{\Gamma^{2}(\frac{4}{5}-\frac{1}{n+5})}\,f_{1}\big{(}\tfrac{{\tt
k}}{2},\tfrac{1}{5}+\tfrac{1}{n+5}\big{)}+o(1)\qquad\qquad(N\to\infty)\,,$
(4.37)

where

$\displaystyle h_{1}^{(N,{\rm
reg})}=\sum_{m=1}^{N/2}(\zeta_{m})^{-1}+\frac{N}{2r\cos(\frac{\pi}{r}+\frac{\pi}{n+r})}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-1}$ (4.38)

with $r=5$ and the function $f_{1}$ is defined in formula (A.1) in Appendix A.
The dependence of the l.h.s. on $c$, the twist parameter ${\tt k}$ and the RG
invariant $\mathfrak{a}_{1}$ can be investigated numerically. One finds that

$\tilde{c}_{1}=P_{0}\,c^{3}+P_{3}\,\mathfrak{a}_{1}\qquad\qquad{\rm
for}\qquad\qquad n>\tfrac{5}{2}$ (4.39)

with certain $n$ dependent constants $P_{0}$ and $P_{3}$. Since
$\tilde{c}_{1}$ must vanish as $\mathfrak{a}_{1}=-c^{3}Z_{1}^{(\mu=2)}$, one
obtains

$P_{0}/P_{3}=\tfrac{\nu}{5}\,\big{(}\mathfrak{g}_{6}+\tfrac{1}{5}\,(\nu+1)\,\mathfrak{g}_{2}\mathfrak{g}_{4}+\tfrac{1}{150}\,(\nu+1)(\nu-4)\,\mathfrak{g}_{2}^{3}\big{)}\
.$ (4.40)

Also notice that when $c=0$, the RG invariants
$\mathfrak{a}_{2}=\mathfrak{a}_{4}=0$ and one gets back the case with $\mu=1$
discussed above. This yields the relation $P_{3}=1/Z_{1}^{(\mu=1)}$ or,
explicitly,

$P_{3}=\frac{5\,\Gamma(\frac{5}{2n})}{\Gamma(\frac{1}{2}+\frac{5}{2n})}\,\frac{\Gamma(\frac{7}{10}-\frac{1}{2n})}{\Gamma(\frac{1}{5}-\frac{1}{2n})}\
.$ (4.41)

In fig. 2, numerical data for $P_{0}$ and $P_{3}$ obtained from the study of
the sum rule (4.37) with $\tilde{c}_{1}$ as in (4.39) is compared against the
predictions (4.40), (4.41). One observes that there is apparent deviation
between the numerical results and the analytic expressions as $n$ becomes
smaller than $\frac{5}{2}$. This is due to the slow decay of the remainder
term $o(1)$ in eq. (4.37). A simple fit of the form $b_{1}+b_{2}\,N^{-\delta}$
of the numerical data, which worked well for large $n$, yields a small value
of the exponent $\delta$ and, in fact, becomes highly unstable. As such, we
were not able to determine with confidence the limiting value of the l.h.s. of
(4.37) as $n<\frac{5}{2}$. In all likelihood, this means that its large - $N$
behaviour can not be adequately described as $b_{1}+b_{2}\,N^{-\delta}$
containing a single small exponent $\delta$.

$n$$P_{3}$$P_{0}/P_{3}$$n$ Figure 2: The crosses represent numerical data for
the ratio $P_{0}/P_{3}$ (left panel) and $P_{3}$ (right panel) obtained
through the study of the sum rule for the Bethe roots (4.37) - (4.39). The
scaling limit is taken for the spin chain with $r=5$, where the RG invariants
are specified as in eqs. (4.17) and (4.18) with $A=1$, $j=0$ and $\mu=2$. The
values of $\mathfrak{a}_{2}$ and $\mathfrak{a}_{4}$ were fixed according to
(4.27), $\mathfrak{a}_{3}$ was set to zero, while the computation was
performed for various values of $\mathfrak{a}_{1}$. The solid line is a plot
of the analytic formula given by (4.40) and (4.41) for the left and right
panels, respectively. Note that there is a loss of accuracy as
$n\lesssim\frac{5}{2}$.

####  4.2.4 Example: $r=7$, $A=1$ and $(\mu,j)=(1,0)$

To understand better the problems that occurred for the cases with $r=5$
considered above, it is instructive to turn to another example, where $r=7$,
$A=1$ and $(\mu,j)=(1,0)$. Then there is only a single non-trivial RG
invariant

$\displaystyle{\mathfrak{a}}_{1}=\bigg{(}\frac{N}{rN_{0}}\bigg{)}^{d_{1}}\
\frac{1}{r}\ \sum_{\ell=1}^{r}(\eta_{\ell})^{-1}$ (4.42)

with $r=7$ and $d_{1}=\frac{10}{7}$. The scaling limit of the Bethe roots for
the ground state is described by the ODE

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+7)y}+E^{7}\
\mbox{e}^{7y}-c\,E\,\mbox{e}^{\frac{1}{7}(n+14)\,y}\,\bigg{]}\,\psi=0\,,$
(4.43)

where $c$, therein, is related to the value of $\mathfrak{a}_{1}$ as

$\mathfrak{a}_{1}=\frac{\Gamma(\frac{1}{2}+\frac{7}{2n})}{7\,\Gamma(\frac{7}{2n})}\,\frac{\Gamma(\frac{1}{7}-\frac{3}{2n})}{\Gamma(\frac{9}{14}-\frac{3}{2n})}\,c\
.$ (4.44)

The proportionality coefficient is free of singularities and zeroes as
$n>\frac{21}{2}$. It turns out that when the position of the first pole in the
$\mathfrak{a}_{s}$ — $c$ relation occurs at sufficiently large values of $n$,
one can obtain more accurate numerical data in its vicinity. In particular, we
found with confidence that

$\bigg{(}\frac{N}{7N_{0}}\bigg{)}^{-\frac{2n}{7(n+7)}}\,h_{1}^{(N,{\rm
reg})}=c\
\frac{(n+7)^{\frac{2}{n+7}-\frac{12}{7}}}{\Gamma^{2}(\frac{6}{7}-\frac{1}{n+7})}\,f_{1}\big{(}\tfrac{{\tt
k}}{2},\tfrac{1}{7}+\tfrac{1}{n+7}\big{)}+o(1)\qquad\qquad(n>\tfrac{21}{2})\,,$
(4.45)

where

$o(1)=\begin{cases}O\big{(}c\,N^{\frac{6}{7}-\frac{4n}{49}}\big{)}&\ \ \ {\rm
for}\ \ \ \frac{21}{2}<n<35\\\\[5.69046pt] O\big{(}c\,N^{-2}\big{)}&\ \ \ {\rm
for}\ \ \ n>35\end{cases}\,.$ (4.46)

It follows from the relation (4.44) that when $n$ approaches $\frac{21}{2}$
from above with the RG invariant $\mathfrak{a}_{1}$ kept fixed, $c$ must tend
to zero. Then, the differential equation (4.43) becomes the one corresponding
to the ${\cal Z}_{7}$ invariant case at $n=\frac{21}{2}$. We also performed a
numerical study in the domain $0<n\leq\frac{21}{2}$. It was found that if
$\mathfrak{a}_{1}$ is still specified as in (4.42) with $d_{1}=\frac{10}{7}$,
the scaling limit of the Bethe roots does not depend on the value of this RG
invariant and is the same as for $\mathfrak{a}_{1}=0$. However, it impacts the
rate of convergence, in particular,

$\bigg{(}\frac{N}{7N_{0}}\bigg{)}^{-\frac{2n}{7(n+7)}}\,h_{1}^{(N,{\rm
reg})}=O\big{(}c\,N^{-\frac{6}{7}+\frac{4n}{49}}\big{)}\qquad\qquad(0<n<\tfrac{21}{2})\,.$
(4.47)

In fig. 3 the numerical data in support of eqs. (4.46) and (4.47) is
presented.

$\delta$$n$

Figure 3: The figure corresponds to the case $r=7$, $A=1$ and $(\mu,j)=(1,0)$.
In order to numerically analyze the correction terms in eqs. (4.45) and (4.47)
the quantity $(\frac{N}{7N_{0}})^{-\frac{2n}{7(n+7)}}\,h_{1}^{(N,{\rm reg})}$
was computed for increasing $N$. This was fitted with $h_{1}^{(\infty)}+{\rm
const}\times N^{-\delta}$, where $h_{1}^{(\infty)}$ stands for the known
limiting value, i.e., given in (4.45) for $n>\frac{21}{2}$ and zero for
$n<\frac{21}{2}$. The results for $\delta$ are represented by the blue dots.
The green triangles depict data for $\delta$ that was obtained by applying the
similar procedure to $(\frac{N}{7N_{0}})^{-\frac{2n}{7(n+7)}}\,(h_{1}^{(N,{\rm
reg})}|_{{\tt k}=0.02}-h_{1}^{(N,{\rm reg})}|_{{\tt k}=0.04})$ instead. The
solid red line is a plot of $\delta=|\frac{6}{7}-\frac{4n}{49}|$. Note that
near $n\approx 10.5$ the value of the exponent $\delta$ becomes small and the
simple fit with a single exponent seems to be unreliable. For $n>35$ one
observes $\delta=2$ with good accuracy. Again, in the vicinity $n\lesssim 35$
the fit with a single exponent looks inadequate, which may explain the
significant spread in the numerical data obtained via the fitting procedures.

We explored the possibility of achieving a non-trivial scaling limit by making
the exponent in (4.42) $n$ \- dependent such that $0<d_{1}<\frac{10}{7}$. It
was found that if $d_{1}$ is replaced by $\tilde{d}_{1}$ where

$\tilde{d}_{1}=\frac{4}{7}+\frac{4n}{49}\qquad\qquad(0<n<\tfrac{21}{2})$
(4.48)

then the limits

$h_{i}^{(\infty)}=\lim_{N\to\infty}\bigg{(}\frac{N}{7N_{0}}\bigg{)}^{-\frac{2in}{7(n+7)}}\,h_{i}^{(N,{\rm
reg})}\qquad\qquad\qquad(i=1,2,\ldots)$ (4.49)

exist and are non-vanishing for generic values of the RG invariant
$\mathfrak{a}_{1}$ in (4.42). This suggests that the Bethe roots for the
ground state exhibit the scaling behaviour.

####  4.2.5 The case $\frac{rj}{A}<\mu<\frac{r}{A+1}\,(j+\frac{1}{2})$

For $A=1$, and $(\mu,j)=(1,0)$ it is possible to give a uniform description
for any $r\geq 5$, which incorporates the case $r=7$ just discussed. With
regards to this, it should be pointed out that $\mathfrak{a}_{1}$ has exponent
$d_{1}=\frac{2}{r}(r-2)$, while formula (4.18) implies that for all the other
RG invariants $d_{s}>2$. As was already mentioned, their values do not affect
the scaling limit and they can be set to zero. This way, the scaling of the
Bethe roots depends only on the single RG invariant $\mathfrak{a}_{1}$ (4.42).
Similar to what was discussed for $r=7$ we confirmed that, in order to achieve
a non-trivial scaling limit, the exponent $d_{1}$ should be taken as

$d_{1}=\begin{cases}\frac{2}{r}\,(r-2)\ \ &{\rm for}\ \ \
n>\frac{1}{2}\,r\,(r-4)\\\\[5.69046pt] \frac{4}{r^{2}}\,(n+r)\ \ &{\rm for}\ \
\ 0<n<\frac{1}{2}\,r\,(r-4)\end{cases}\,.$ (4.50)

For $n>\frac{1}{2}\,r\,(r-4)$ the scaled Bethe roots are described by the
differential equation

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}+E^{r}\
\mbox{e}^{ry}-c\,E\,\mbox{e}^{\frac{1}{r}(n+2r)\,y}\,\bigg{]}\,\psi=0\,.$
(4.51)

We believe that formulae (4.50) and (4.51) are applicable for $A=1$,
$(\mu,j)=(1,0)$ and any $r\geq 5$, despite that for $r=5$ they are difficult
to confirm numerically in the domain $0<n<\frac{5}{2}$.

If one accepts that all the RG invariants with $d_{s}\geq 2$ may be ignored,
an interesting prediction is reached for the case of generic values of $r$,
$A$ and $j$ with $\mu$ restricted by the condition

$\frac{rj}{A}<\mu<\frac{r}{A+1}\ \Big{(}j+\frac{1}{2}\Big{)}\,.$ (4.52)

The latter implies that the only non-trivial RG invariant is
$\mathfrak{a}_{\mu}$ with the corresponding exponent
$d_{\mu}=\frac{2}{r}\,\big{(}(j+1)r-(A+1)\,\mu\big{)}$ such that
$1<d_{\mu}<2$. It is expected that a scaling limit, which is different than
the one for the ${\cal Z}_{r}$ invariant case, can be achieved provided

$d_{\mu}=\begin{cases}\frac{2M}{r}\ \ \ \ \ \ \ \ &{\rm for}\ \ \ n>n_{\rm
min}\\\\[5.69046pt] \frac{4Ln}{r^{2}}+\frac{2}{r}(r-M)\ \ \ \ \ \ \ &{\rm
for}\ \ \ 0<n<n_{\rm min}\end{cases}\,,$ (4.53)

where

$n_{\rm min}=\tfrac{r}{L}\,(M-\tfrac{r}{2})\qquad\qquad(L>0)$ (4.54)

and we use the integers $M$ and $L$ defined in eqs. (4.4) and (4.10),
respectively. For $n>n_{\rm min}$, the differential equation appearing in the
scaling limit would be (4.1). If in the domain $0<n<n_{\rm min}$ the exponent
$d_{\mu}$ is not assigned an $n$ dependence but kept fixed as
$d_{\mu}=\frac{2M}{r}$, the scaled Bethe roots would be the same as for the
${\cal Z}_{r}$ invariant case. When the exponent is modified according to the
second line of (4.53) we observed in many examples that the Bethe roots
exhibit non-trivial scaling behaviour. Describing the corresponding ODE for
$0<n<n_{\rm min}$ is beyond the scope of this work.

It should be emphasized that formula (4.54) does not make sense for $L=0$.
This is again a signal that the case $A\mu=rj$ requires special attention.

Let us give a short summary of this section. For sufficiently large $n>n_{\rm
min}$, the scaling limit can be organized such that the scaled Bethe roots are
described by the differential equation (4.1), (4.2). The specialization of the
RG invariants is given by eqs. (4.17), (4.20) with exponents (4.18). The
values of the non-vanishing $\mathfrak{a}_{s}$ are related to the parameters
of the differential equation as in (4.22) - (4.26). The values of some of the
RG invariants with $d_{s}>1$ turn out to have no effect on the scaling limit
and can be taken to be zero. For given $A$ and $(\mu,j)$ obeying the
inequalities (4.2), the minimal sets of non-trivial $\mathfrak{a}_{s}$, along
with their exponents, are listed in the tables contained in Appendix F for
$3\leq r\leq 10$. Therein, we also quote the corresponding values of the lower
bound $n_{\rm min}$.

##  5 The case $A\mu=rj$ with $A=1,\ldots,r-2$

A formal specification (ignoring the constraint (4.2)) of the ODE (4.1) to the
case at hand yields

$\displaystyle\Big{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-c\,E^{\mu}\ \mbox{e}^{\mu y}\,\Big{]}\,\psi=0\ .$ (5.1)

In view of the condition $A\mu=rj$, $\mu$ is not an arbitrary integer from the
segment $[1,r-1]$. Instead, if $H$ is the greatest common divisor of $r$ and
$A$,

$H={\rm gcd}(r,A)\,,$ (5.2)

then

$\displaystyle\mu=\frac{rJ}{H}\ \ \ \ \ \ \ (J=1,\ldots,H-1)\ .$ (5.3)

It turns out that the case $A\mu=rj$ with $A=1,\ldots,r-2$ falls outside of
the scope of the discussion in the previous section and the ODE (5.1) does not
actually appear in the scaling limit of the spin chain (except when both $r,A$
are even and $\mu=\frac{r}{2}$).

###  5.1 Special features of the scaling limit

Let’s first consider the case when the positive integer $\mu<\frac{r}{2}$ or,
equivalently, $J<\frac{H}{2}$. Similarly to what was discussed in sec. 4.2.5
one may expect that all of the RG invariants have no effect on the scaling
limit apart from $\mathfrak{a}_{\mu}$, whose corresponding exponent
$d_{\mu}=\frac{2}{H}\,(H-J)$ is such that $1<d_{\mu}<2$. In fact, a numerical
study shows that the scaling behaviour does not depend on the value of
$\mathfrak{a}_{\mu}$ as well, and the scaled Bethe roots coincide with those
for the ${\cal Z}_{r}$ invariant case.888A non-trivial scaling behaviour for
the Bethe roots for the ground state can be achieved by swapping
$d_{\mu}=\frac{2}{H}\,(H-J)$ for $\tilde{d}_{\mu}=\frac{2J}{H}<1$ in the
definition of the RG invariant $\mathfrak{a}_{\mu}$ (4.17). This follows from
the relation between the case $A\mu=rj>0$ and $A=0$ described in sec. 5.2
together with the results of sec. 7.

For $r$ even, $\mu=\frac{r}{2}$ is among the admissible values (5.3). Then it
is sufficient to keep only one RG invariant $\mathfrak{a}_{\frac{r}{2}}$ to be
non-vanishing with exponent $d_{\frac{r}{2}}=1$. However formulae
(4.22)-(4.26) for the coefficient $Z_{\frac{r}{2}}$, which relates the value
of $\mathfrak{a}_{\frac{r}{2}}$ with the parameter $c$ in the differential
equation, give infinity. This suggests to modify the definition of the RG
invariant as

$\displaystyle{\mathfrak{a}}_{\frac{r}{2}}\equiv\frac{2}{r^{2}}\
\frac{1}{\log\big{(}\frac{N}{rN_{0}}\big{)}}\
\bigg{(}\frac{N}{rN_{0}}\bigg{)}\
\sum_{\ell=1}^{r}(\eta_{\ell})^{-\frac{r}{2}}\ .$ (5.4)

With such a setup, the investigation of the scaling behaviour becomes
difficult to carry out numerically. Even $N\sim 5000$ turns out to be
insufficient to reliably observe the existence of a scaling limit.
Nevertheless, in Appendix E we present a numerical scheme, which dramatically
improves convergence and allows one to confirm that the differential equation
occurring in the scaling limit is given by (5.1) with $\mu=\frac{r}{2}$. The
parameter $c$ is related to the value of the RG invariant as

$\displaystyle{\mathfrak{a}}_{\frac{r}{2}}=(-1)^{\frac{1}{2}(r-A+2)}\
\frac{2}{r\sqrt{\pi}}\
\frac{\Gamma(\frac{1}{2}+\frac{r}{2n})}{\Gamma(1+\frac{r}{2n})}\ c\ .$ (5.5)

Note that the argument of the logarithm in (5.4) involves the constant
$rN_{0}$. It can be replaced by any positive number without affecting the
scaling limit.

In the case $\mu>\frac{r}{2}$ the integer $J$ in formula (5.3) is restricted
as $\frac{H}{2}<J\leq H-1$. It turns out that the minimal set of non-trivial
RG invariants consists of

$\displaystyle\mathfrak{a}_{s}:\ \ s=\tfrac{rJ}{H}\,i_{s}\ \ \ \ ({\rm mod}\
r)\,,\qquad\qquad i_{s}=1,2,\ldots,\big{[}\tfrac{H}{2(H-J)}\big{]}\,.$ (5.6)

The corresponding exponents are given by

$\displaystyle d_{s}=\frac{2}{H}\,(H-J)\,i_{s}$ (5.7)

and are all less than or equal to one. A numerical study shows that the Bethe
roots for the ground state develop a scaling behaviour, which differs from the
one that occurs in the ${\cal Z}_{r}$ invariant case. However, the scaling
limit is not described by the ODE (5.1). For instance, when $r=9$, $A=3$
($H=3$) and $\mu=6$ there is only one non-trivial RG invariant
$\mathfrak{a}_{6}$ with exponent $d_{6}=\frac{2}{3}$. The differential
equation yields the prediction that the limits (2.59) with $s=1,2,\ldots,8$
vanish except for

$h_{6}^{(\infty)}=\lim_{N\to\infty}\bigg{(}\frac{9N_{0}}{N}\bigg{)}^{\frac{4n}{3(n+9)}}\,h_{6}^{(N,{\rm
reg})}\ .$ (5.8)

However, we found that together with $h_{6}^{(\infty)}$, the limit

$h_{3}^{(\infty)}=\lim_{N\to\infty}\bigg{(}\frac{9N_{0}}{N}\bigg{)}^{\frac{2n}{3(n+9)}}\,h_{3}^{(N,{\rm
reg})}$ (5.9)

exists and is not zero.

The first idea that may come to mind for describing the scaling limit for
$\mu=\frac{rJ}{H}>\frac{r}{2}$ is to replace (5.1) by the more general
differential equation of the form

$\Big{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-c\,E^{\mu}\ \mbox{e}^{\mu y}-\sum_{i}{b}_{i}\,E^{\mu_{i}}\
\mbox{e}^{\mu_{i}y}\,\Big{]}\,\psi=0\,,$ (5.10)

where the summation runs over the positive integer
$i=\big{[}\tfrac{H}{2(H-J)}\big{]}+1,\,\ldots,\,\big{[}\tfrac{H-1}{H-J}\big{]}$
and

$\displaystyle\mu_{i}=i\,\mu-(i-1)\,r\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\big{(}\mu=\tfrac{rJ}{H}\big{)}\ .$ (5.11)

Numerical data supports this conjecture, where the values of the RG invariants
$\mathfrak{a}_{s}$ (5.6) are related to $c$ as in eqs. (4.22)-(4.26). However,
while the coefficients $b_{i}$ depend simply on $c$, they turn out to be
complicated functions of $p=\frac{1}{2}(n+r)\,{\tt k}$ and $n$:

$b_{i}=c^{i}\ {\hat{b}}_{i}(p,n)\,.$ (5.12)

For instance, the differential equation (5.10), specialized to the case $r=9$,
$A=3$ ($H=3$) and $\mu=6$, leads to the prediction

$\displaystyle h^{(\infty)}_{3}$ $\displaystyle=$ $\displaystyle
b_{1}\,\frac{(n+9)^{-\frac{2(n+6)}{n+9}}}{\Gamma^{2}\big{(}\frac{n+6}{n+9}\big{)}}\,f_{1}\big{(}\tfrac{{\tt
k}}{2},\tfrac{3}{n+9}\big{)}$ (5.13) $\displaystyle h^{(\infty)}_{6}$
$\displaystyle=$ $\displaystyle
c\,\frac{(n+9)^{-\frac{2(n+3)}{n+9}}}{\Gamma^{2}\big{(}\frac{n+3}{n+9}\big{)}}\,f_{1}\big{(}\tfrac{{\tt
k}}{2},\tfrac{6}{n+9}\big{)}+b^{2}_{1}\
\frac{(n+9)^{-\frac{4(n+6)}{n+9}}}{\Gamma^{4}\big{(}\frac{n+6}{n+9}\big{)}}\,f_{2}\big{(}\tfrac{{\tt
k}}{2},\tfrac{3}{n+9}\big{)}\,.$

Eliminating $b_{1}$, one obtains

$\displaystyle h^{(\infty)}_{6}$ $\displaystyle=$ $\displaystyle c\
\frac{(n+9)^{-\frac{2(n+3)}{n+9}}}{\Gamma^{2}\big{(}\frac{n+3}{n+9}\big{)}}\,f_{1}\big{(}\tfrac{{\tt
k}}{2},\tfrac{6}{n+9}\big{)}+\Bigg{(}\frac{h^{(\infty)}_{3}}{f_{1}\big{(}\tfrac{{\tt
k}}{2},\tfrac{3}{n+9}\big{)}}\Bigg{)}^{2}\,f_{2}\big{(}\tfrac{{\tt
k}}{2},\tfrac{3}{n+9}\big{)}\ .$ (5.14)

Here the functions $f_{1}$ and $f_{2}$ are given by eqs. (A.1)-(A) in the
Appendix A, while the $\mathfrak{a}_{s}$ — $c$ relation reads explicitly as

$\displaystyle\mathfrak{a}_{6}=-\frac{\Gamma\left(\frac{3}{2n}\right)\Gamma\left(\frac{3}{2}+\frac{9}{2n}\right)}{(n+9)\Gamma\left(1+\frac{9}{2n}\right)\Gamma\left(\frac{1}{2}+\frac{3}{2n}\right)}\
c\ .$ (5.15)

Numerical data in support of (5.14) is presented in tab. 1.

In tab. 2 we list the cases where $A\mu=rj$ and $\mu>\frac{r}{2}$ for all
possible values of $r\leq 10$. Together with the integers $s$ labeling the
non-trivial RG invariants and the corresponding exponents $d_{s}$, presented
are the values of the integers $\mu_{i}$ entering into the differential
equation (5.10). Note that the two cases $A=3$ and $A=6$ with $r=9$ lead to
the same ODE. In fact, it is easy to check that the two sets of Bethe roots
for the ground state are related as:
$\\{\zeta_{j}\\}|_{A=3}=\\{\mbox{e}^{\frac{{\rm
i}\pi}{3}}\,\zeta_{j}\\}|_{A=6}$.

| ${\tt k}=0.02$ | ${\tt k}=0.04$
---|---|---
$n$ |  $h_{6}^{(\infty)}$ from (5.8) |  $h_{6}^{(\infty)}$ from (5.14) |  $h_{6}^{(\infty)}$ from (5.8) |  $h_{6}^{(\infty)}$ from (5.14)
10 | 0.189683 | 0.189499 | 0.182812 | 0.182641
20 | 0.096347 | 0.096342 | 0.089699 | 0.089695
30 | 0.070086 | 0.070084 | 0.063443 | 0.063442
40 | 0.057235 | 0.057232 | 0.050550 | 0.050548
Table 1: Compared is numerical data obtained for $h_{6}^{(\infty)}$ in two different ways for $r=9$, $A=3$ and various values of ${\tt k}$ and $n$. In the first case, $h_{6}^{(\infty)}$ was calculated from the solution of the Bethe Ansatz equations via formulae (5.8), (2.55) and (2.36). In the second, it was found from the r.h.s. of (5.14), where $h_{3}^{(\infty)}$ was calculated numerically using eq. (5.9). The value of the single non-trivial RG invariant $\mathfrak{a}_{6}$, with exponent $d_{6}=\frac{2}{3}$, was fixed according to (5.15), where $c=1$. Computations were performed for spin chains with $N\lesssim 5000$.  |  | $\mu$ | $s$ | $d_{s}$ | $\mu_{i}$
---|---|---|---|---|---
$\ r=6$ | $\ A=3\ $ | 4 | $\\{4\\}$ | $d_{4}=\frac{2}{3}$ | $\mu_{2}=2$
$\ r=8$ | $\ A=4\ $ | 6 | $\\{6,4\\}$ | $d_{6}=\frac{1}{2},\,d_{4}=1$ | $\mu_{3}=2$
$r=9$ | $\ A=3\ $ | 6 | $\\{6\\}$ | $d_{6}=\frac{2}{3}$ | $\mu_{2}=3$
$\ A=6\ $ | 6 | $\\{6\\}$ | $d_{6}=\frac{2}{3}$ | $\mu_{2}=3$
$r=10$ | $\\!\\!A=5$ | 6 | $\\{6\\}$ | $d_{6}=\frac{4}{5}$ | $\mu_{2}=2$
8 | $\\{8,6\\}$ | $d_{8}=\frac{2}{5},\,d_{6}=\frac{4}{5}$ | $\mu_{3}=4,\,\mu_{4}=2$
Table 2: Listed are all possible cases where $A\mu=rj$ with $\mu>\frac{r}{2}$
up to $r=10$. Apart from the integers $s$ and exponents $d_{s}$, which specify
the RG invariants $\mathfrak{a}_{s}$ belonging to the set (5.6), the values of
the integers $\mu_{i}$ (5.11) entering into the differential equation (5.10)
are given in the last column.

In sec. 4.1 a subtlety was mentioned that occurs when both $r$, $A$ are even
with $\mu$ and $j$ given by (4.11). Setting $K=1$ into the formula
$\mu=\frac{r}{2K}\,(2K-1)$ one gets back the case $\mu=\frac{r}{2}$ discussed
above. If $K>1$ then $\mu>\frac{r}{2}$. For $r$ even we expect that the RG
invariant $\mathfrak{a}_{\frac{r}{2}}$ should be specified as in (5.4).
Moreover, it seems likely that (5.5) is generalized as

$\displaystyle{\mathfrak{a}}_{\frac{r}{2}}=(-1)^{\frac{1}{2}(r-A+2)}\ \
\frac{2\Gamma(K-\frac{1}{2})}{r\pi K!}\
\frac{\Gamma(\frac{1}{2}+\frac{r}{2n})}{\Gamma(1+\frac{r}{2n})}\ \ c^{K}\ .$
(5.16)

Unfortunately, we were not able to numerically verify this relation for $K>1$.
Notice that the first such case corresponds to $r=8$, $A=4$ and $\mu=6$, see
tab. 2.

###  5.2 Relation to $A=0$

It turns out that the cases $A\mu=rj$ with $A=1,2,\ldots,r-2$ and $A=0$ are
related. The link already occurs for the Bethe roots corresponding to the
ground state of the spin chain at finite lattice size $N$. It allows one to
carry over the discussion of the case $A=0$ contained in sec. 7 which, among
other things, gives a relation between the coefficients $b_{i}$ in the ODE
(5.10) with the eigenvalues of the so-called quasi-shift operators.

Recall that the integer $\mu$ is used to specify the minimal set of non-
trivial RG invariants $\mathfrak{a}_{s}$. In view of (4.20), when performing
the scaling limit, one may take

$\sum_{\ell=1}^{r}(\eta_{\ell})^{-s}=0\qquad\qquad{\rm for}\qquad\qquad
s\neq\sigma,2\sigma,\ldots,r\ ,$ (5.17)

where $\sigma>1$ stands for the greatest common divisor of $\mu$ and $r$:

$\sigma={\rm gcd}(r,\mu)\,.$ (5.18)

The condition (5.17) may be fulfilled if one sets

$\eta_{\ell+\frac{r}{\sigma}}=\mbox{e}^{\frac{2\pi{\rm
i}}{\sigma}}\,\eta_{\ell}\ .$ (5.19)

When $s$ is a multiple of $\sigma$, i.e., $s=\sigma m$ with
$m=1,2,\ldots,r/\sigma-1$, then instead of (5.17) one has

$\frac{1}{rm}\sum_{\ell=1}^{r/\sigma}(\eta_{\ell})^{-\sigma
m}=\begin{cases}\big{(}\frac{rN_{0}}{N}\big{)}^{d_{\sigma
m}}\,\mathfrak{a}_{\sigma m}\qquad\qquad&{\rm if}\qquad\sigma m\in
S_{+}\cup\bar{S}_{-},\ \sigma m\neq\frac{r}{2}\\\\[5.69046pt]
\frac{1}{\log(\frac{N}{rN_{0}})}\,\big{(}\frac{rN_{0}}{N}\big{)}\,\mathfrak{a}_{\frac{r}{2}}\qquad\qquad&{\rm
if}\qquad\sigma m=\frac{r}{2}\in S_{+}\\\\[5.69046pt] 0\qquad&{\rm
otherwise}\end{cases}\,,$ (5.20)

while the normalization condition (4.21) becomes

$\frac{\sigma}{r}\ \sum_{\ell=1}^{r/\sigma}\,(\eta_{\ell})^{-r}=(-1)^{r-1}\,.$
(5.21)

For given values of $\mathfrak{a}_{\sigma m}$ and $N$, this may be treated as
a system of $r/\sigma$ equations for determining the set
$\\{\eta_{1},\ldots,\eta_{r/\sigma}\\}$ and hence, in view of (5.19), all the
inhomogeneities.

The $N/2$ Bethe roots for the ground state ($S^{z}=0$) are split into groups
with roughly equal phases, which can be labelled by the integer
$a=1,2,\ldots,r$ such that
$\arg(\zeta_{m}^{(a)})\approx\frac{\pi}{r}\,(2a-2-A)\mod 2\pi$. If (5.19) is
obeyed and

$\zeta_{m}^{(a+\frac{r}{\sigma})}=\mbox{e}^{\frac{2\pi{\rm
i}}{\sigma}}\,\zeta_{m}^{(a)}\,,$ (5.22)

then the algebraic system of $N/2$ Bethe Ansatz equations (1.11) reduces to a
system of $\tilde{N}/2$ equations with

$\tilde{N}=N/\sigma\ .$ (5.23)

Namely,

$\Bigg{(}\prod_{\ell=1}^{\tilde{r}}\frac{\tilde{\eta}_{\ell}+\tilde{q}^{+1}\,\,\tilde{\zeta}_{j}}{\tilde{\eta}_{\ell}+\tilde{q}^{-1}\,\tilde{\zeta}_{j}}\Bigg{)}^{\frac{\tilde{N}}{\tilde{r}}}=-\omega^{2}\,\prod_{i=1}^{\tilde{N}/2}\,\frac{\tilde{\zeta}_{i}-\tilde{q}^{+2}\,\tilde{\zeta}_{j}}{\tilde{\zeta}_{i}-\tilde{q}^{-2}\,\tilde{\zeta}_{j}}\,\qquad\qquad\big{(}j=1,2,\ldots,\tfrac{\tilde{N}}{2}\big{)}\,,$
(5.24)

where

$\\{\tilde{\zeta}_{j}\\}_{j=1}^{\tilde{N}/2}=\bigcup_{a=1}^{\tilde{r}}\,\Big{\\{}(-1)^{I}\,\big{(}\zeta_{j}^{(a)}\big{)}^{\sigma}\Big{\\}}_{j=1}^{N/(2r)}$
(5.25)

and we also use

$\tilde{r}=r/\sigma\,,\qquad\qquad\quad\tilde{q}=(-1)^{I}\,q^{\sigma}\,,\qquad\qquad\quad\tilde{\eta}_{\ell}=(-1)^{\sigma-1}\,(\eta_{\ell})^{\sigma}\qquad\quad(\ell=1,2,\ldots,\tilde{r})\,.$
(5.26)

The above formulae involve the sign factor $(-1)^{I}$ with
$I=\frac{A\sigma}{r}$, which is an integer due to the condition $A\mu=rj$.

From the definition (5.26), it follows that

$\tilde{q}=\exp\Big{(}{\frac{{\rm
i}\pi}{\tilde{n}+\tilde{r}}}\Big{)}\qquad{\rm
with}\qquad\tilde{n}=n/\sigma\,,$ (5.27)

i.e., the argument of $\tilde{q}$ lies in the interval
$(0,\frac{\pi}{\tilde{r}})$. Thus the algebraic system (5.24) coincides with
the Bethe Ansatz equations for the $\tilde{r}$ \- site periodic spin chain
with the set of inhomogeneities
$\\{\tilde{\eta}_{\ell}\\}_{\ell=1}^{\tilde{r}}$ in the regime $\tilde{A}=0$.
In view of (5.20), the RG invariants $\tilde{\mathfrak{a}}_{s}$ are expressed
in terms of the original ones as

$\tilde{\mathfrak{a}}_{s}=(-1)^{(\sigma-1)s}\,\sigma^{1-d_{\sigma
s}}\,\mathfrak{a}_{\sigma s}\qquad\qquad{\rm
and}\qquad\qquad\tilde{d}_{s}=d_{\sigma s}\ .$ (5.28)

It turns out that if $\\{\zeta_{j}\\}$ is the solution to the Bethe Ansatz
equations (1.11) corresponding to the ground state, then $\tilde{\zeta}_{j}$
(5.25) are the Bethe roots for the ground state of the spin chain with
$(N,r,q,\eta_{\ell})$ swapped for
$(\tilde{N},\tilde{r},\tilde{q},\tilde{\eta}_{\ell})$ defined in eqs. (5.23)
and (5.26).

##  6 General differential equation for $A=1,2,\ldots,r-2$

Consider the equation

$\displaystyle\bigg{[}-\partial_{y}^{2}+p^{2}+\mbox{e}^{(n+r)y}-(-1)^{A}\,E^{r}\
\mbox{e}^{ry}-\sum_{(\mu,j)\in\,\bm{\Xi}_{r,A}}c_{\mu,j}\,E^{\mu}\
\mbox{e}^{\big{(}(A\mu-rj)\,\frac{n+r}{r}+\mu\big{)}y}\,\bigg{]}\,\psi=0\ ,$
(6.1)

where the sum goes over pairs of non-negative integers belonging to the set

$\displaystyle\bm{\Xi}_{r,A}$ $\displaystyle=$
$\displaystyle\big{\\{}(\mu,j)\,:\ \ \tfrac{rj}{A}<\mu<\tfrac{r}{A+1}\ (j+1)\
\ \&\ \ j\geq 0\big{\\}}\ .$ (6.2)

We have excluded $(\mu,j)$ for which the condition $A\mu=rj$ is satisfied for
the reasons discussed in the previous section. It turns out that for any
$r\geq 3$ and $A=1,2,\ldots,r-2$ the number of admissible pairs obeys the
condition

$|\bm{\Xi}_{r,A}|=|\bm{\Xi}_{r,r-1-A}|$ (6.3)

and is bounded as999We checked (6.4) up to $r=200$. If $r$ is a prime number
then $|\bm{\Xi}_{r,A}|=\frac{r-1}{2}$ for any $A=1,2,\ldots,r-2$, while as
$A=\frac{r}{2}-1$ and $A=\frac{r}{2}$ ($r$ even) one has
$|\bm{\Xi}_{r,A}|=[\tfrac{r}{4}]$.

$[\tfrac{r+1}{4}]\leq|\bm{\Xi}_{r,A}|\leq[\tfrac{r-1}{2}]\ .$ (6.4)

This way, the amount of terms in the sum in (6.1) is less than the number of