aainstitutetext: Department of Physics and Astronomy,
University of Kentucky,
Lexington, KY, USA 40506
bbinstitutetext: Department of Physics, Nagoya University, Nagoya, Aichi
464-8602, Japan
ccinstitutetext: Institute for Advanced Research, Nagoya University, Nagoya,
Aichi 464-8601, Japan

# Spectrum of quantum KdV hierarchy in the semiclassical limit

Anatoly Dymarsky a Ashish Kakkar  Kirill Pavlenko a, b, c and Sotaro Sugishita

###### Abstract

We employ semiclassical quantization to calculate spectrum of quantum KdV
charges in the limit of large central charge $c$. Classically, KdV charges
$Q_{2n-1}$ generate completely integrable dynamics on the co-adjoint orbit of
the Virasoro algebra. They can be expressed in terms of action variables
$I_{k}$, e.g. as a power series expansion. Quantum-mechanically this series
becomes the expansion in $1/c$, while action variables become integer-valued
quantum numbers $n_{i}$. Crucially, classical expression, which is homogeneous
in $I_{k}$, acquires quantum corrections that include terms of subleading
powers in $n_{k}$. At first two non-trivial orders in $1/c$ expansion these
‘‘quantum’’ terms can be fixed from the analytic form of $Q_{2n-1}$ acting on
the primary states. In this way we find explicit expression for the spectrum
of $Q_{2n-1}$ up to first three orders in $1/c$ expansion. We apply this
result to study thermal expectation values of $Q_{2n-1}$ and free energy of
the KdV Generalized Gibbs Ensemble.

##  1 Introduction

Conformal invariance in two dimensions is a very powerful tool which gives
rise to many non-pertubative relations constraining dynamics of 2d CFTs. Among
them is universality of stress-energy tensor sector BPZ , namely any
correlation function which includes only stress-energy tensor and its
descendants depends only on central charge $c$ but not on any other details of
the theory. An analytic form of all such correlators can in principle be found
in a recursive form zamolodchikov1989conformal . The stress-energy sector can
be regarded as integrable, even if the whole theory is understood to be
chaotic hartman2014universal . This can be justified formally by noting there
is an infinite number of mutually commuting quantum KdV charges
bazhanov1996integrable ; bazhanov1997integrable ; bazhanov1999integrable  –
local charges $Q_{2n-1}$ of the form

$\displaystyle Q_{2n-1}={1\over
2\pi}\int_{0}^{2\pi}T_{2n}(\varphi)\,d\varphi,$ (1)

where the densities $T_{2n}$ are appropriately regularized polynomials in
stress-energy tensor $T(\varphi)$ and its derivatives. First charge

$\displaystyle Q_{1}=L_{0}-{c\over 24}={1\over
2\pi}\int_{0}^{2\pi}T\,d\varphi$ (2)

is the CFT Hamiltonian. (Here and below we consider 2d CFT on a cylinder.
Because of standard factorization into left and right-moving sectors we
restrict the discussion to one sector only.) Interest in integrable structure
of 2d CFT stress-energy sector has been reignited recently in the context of
Eigenstate Thermalization Hypothesis (ETH) srednicki1994chaos . Following
original works lashkari2018eigenstate ; lin2016thermality ; basu2017thermality
; He2017txy ; He2017vyf ; lashkari2018universality ; Guo2018pvi ;
maloney2018generalized ; GGE  it has been conjectured and confirmed in GETH
that 2d CFTs exhibit generalized ETH with the local equilibrium being
described by qKdV Generalized Gibbs Ensemble (GGE). Schematically the role of
qKdV charges is as follows. The CFT Hamiltonian (2) is highly degenerate with
all CFT descendant states of the form

$\displaystyle|E\rangle=L_{-m_{1}}\dots
L_{-m_{k}}|\Delta\rangle,\qquad\sum_{i=1}^{k}m_{i}=m$ (3)

sharing the same energy $E=\Delta+m-c/24$. Since all $Q_{2n-1}$ commute, they
can be simultaneously diagonalized giving rise to mathematically unique
‘‘integrable’’ basis of eigenstates. Unlike the energy eigenstates of the form
(3), which fail the ETH, integrable eigenstates carry specific values of
$Q_{2k-1}$-charges and obey generalized ETH. This novel role of qKdV
symmetries motivates the question of ‘‘solving’’ integrable structure, i.e.
evaluating spectrum of qKdV charges and finding integrable eigenstates, which
would allow detailed studies of generalized ETH and qKdV GGE thermodynamics.

In certain sense the question of finding qKdV spectra can be regarded as
solved: there is not one but two distinct ways to write an algebraic Bethe-
ansatz reducing the problem of finding spectra to a bunch of algebraic
equations BazhanovEq ; LitvinovEq . In practice complexity of these equations
grows very rapidly with the level $m$ (3), making this approach useless in the
context of ETH, at least so far. The ETH holds in thermodynamic limit, it may
not and does not hold beyond that regime. Thermodynamic limit assumes the
length of the spatial circle $L$ goes to infinity, with the energy density
$E/L$ kept fixed. Using rescaling, one can always bring the circle to unit
radius, the notations we use throughout the paper. The energy $E$ then must go
to infinity as $L^{2}$ with $L\rightarrow\infty$ being an auxiliary parameter
keeping track of corrections to various ETH-related identities. For any given
primary state $|\Delta\rangle$ this essentially means the descendant level $m$
must be taken to infinity, i.e. we arrive exactly at the limit where algebraic
Bethe equations become most difficult.

A progress was achieved by taking an additional limit of large central charge.
In this case $Q_{2k-1}$-eigenstates, akin to (3), can be parametrized by a set
of natural numbers, which can be conveniently combined into a Young tableau
GGE .111Appearance of $n_{k}$ to parametrize the eigenstates can be understood
from the Virasoro algebra, which in the large $c$ limit reduces to a product
of Heisenberg algebras, with $n_{k}$ being the corresponding quantum numbers
Witten ; Brehm_2020 . It is most convenient to use representation when
$n_{k}\geq 0$ for $k=1,2,\dots$ counts the number of rows of length $k$,

$\displaystyle|n_{i}\rangle\equiv|n_{1},\dots\rangle,\qquad\sum_{k=1}^{\infty}k\,n_{k}=m.$
(4)

We emphasize (4) are eigenstates of $Q_{2n-1}$ and thus differ from (3).
Corresponding eigenvalues at leading order were conjectured in GGE2 222Since
Ref. GGE2  was working in the regime of both large central charge and
thermodynamic limit $Q_{1}\propto cL^{2}$, it only conjectured the term linear
in $n_{k}$, as the $n_{k}$-independent term is $1/L^{2}$ suppressed.

$\displaystyle Q_{2n-1}|n_{i}\rangle$ $\displaystyle=$ $\displaystyle{\sf
Q}_{2n-1}|n_{i}\rangle,$ (5) $\displaystyle{\sf Q}_{2n-1}$ $\displaystyle=$
$\displaystyle\tilde{\Delta}^{n}+\sum_{p=0}^{n-1}\xi_{n}^{p}\,\tilde{\Delta}^{n-1-p}\,\tilde{c}^{p}\left(\sum_{k=1}^{\infty}k^{2p+1}n_{k}+{\zeta(-2p-1)\over
2}\right)+O(\tilde{c}^{n-2}),$ $\displaystyle\xi_{n}^{p}$ $\displaystyle=$
$\displaystyle{(2n-1)\sqrt{\pi}\,\Gamma(n+1)\over\,2\,\Gamma(p+3/2)\Gamma(n-p)},\qquad\tilde{\Delta}=\Delta-{\tilde{c}},\qquad\tilde{c}={c-1\over
24}.$ (6)

Here we assume the scaling when $c\rightarrow\infty$ while
$\tilde{\Delta}/\tilde{c}=h$ is kept fixed. No thermodynamic limit is assumed.
This is the limit of holographic correspondence, when CFT is dual to
semiclassical gravity. The holographic picture provides an easy derivation for
the leading $1/c$ terms in (1) and provides interpretation for $n_{k}$ as the
boson occupation numbers of the boundary gravitons, see Appendix A. From the
mathematical point of view simplicity of eigenstates parametrization with help
of Young tableaux as well as relatively simple form of (1) can be readily
understood from the semiclassical quantization of the co-adjoint orbit of
Virasoro algebra. Indeed, as is explained in Witten  in the large $c$ limit
Virasoro algebra can be understood in quasi-classical terms, as quantization
of the Kirillov-Kostant-Souriau symplectic form. Because of $U(1)$ symmetry
semiclassical quantization of $Q_{1}$ is exact,

$\displaystyle Q_{1}=\tilde{\Delta}+\left(\sum_{k=1}^{\infty}k\,n_{k}-{1\over
24}\right)=\Delta+m-{c\over 24},$ (7)

but for all higher $Q_{2n-1},n>1$ it is not. It is a perturbation series in
$1/\tilde{c}$, which plays the effective role of Planck constant. In this
paper we develop a perturbative scheme to obtain the spectrum of $Q_{2n-1}$ as
a series in $1/\tilde{c}$ expansion and calculate first two non-trivial terms.
The result is summarized in (4).

In the strict $c\rightarrow\infty$ limit when the problem becomes classical,
CFT stress-tensor $T$ can be substituted by an element of the co-adjoint orbit
of Virasoro algebra ${24\over c}u$, where $u$ is a potential of an auxiliary
periodic Schr$\ddot{\rm o}$dinger equation. Then quantum KdV charges (1)
reduce to conventional KdV Hamiltonians of the periodic problem

$\displaystyle Q_{2n-1}={1\over 2\pi}\int_{0}^{2\pi}(u^{n}+\dots)\,d\varphi,$
(8)

which we denote the same as the quantum ones, as it clear from the constant
which, classical or quantum version, we had in mind. For the states with large
but finite level $m$ number of non-zero $n_{k}$ will also be finite. At the
classical level this corresponds to finite-zone potentials $u$, which form a
finite-dimensional symplectic manifold equipped with the structure of a
completely integrable system. Hamiltonians ${Q}_{2n-1}$ can be re-expressed in
terms of the action variables $I_{k}$ and the orbit invariant $h$,

$\displaystyle
Q_{2n-1}=h^{n}+\sum_{k=1}^{\infty}\sum_{j=0}^{n-1}\xi_{n}^{j}\,h^{n-1-j}k^{2j+1}\,I_{k}+O(I^{2})$
(9)

which at semiclassical level become integral quantum numbers
$I_{k}\rightarrow{n_{k}/\tilde{c}}$. It is then easy to see that (9) becomes
(1), up to an overall factor $\tilde{c}^{n}$ and certain corrections. At each
power of $1/\tilde{c}$ classical expression $Q_{2n-1}(h,I_{k})$ predicts only
leading power of $n_{k}$ while all subleading powers are ‘‘quantum
corrections’’ which must be fixed separately.

At leading $1/\tilde{c}$ order quantum correction is just $n_{k}$-independent
constant term proportional to $\zeta(-2p-1)/2$, see (1). It can be fixed
trivially by introducing Maslov index $I_{k}\rightarrow(n_{k}+1/2)/\tilde{c}$,
such that constant term can be formally rewritten as the vacuum energy of
‘‘quantum oscillators’’ with frequencies $\omega_{k}$ and occupation numbers
$n_{k}$

$\displaystyle{\sf
Q}_{2n-1}={\tilde{\Delta}}^{n}+\sum_{k=1}^{\infty}(n_{k}+1/2)\,\omega_{k}+O(\tilde{c}^{n-2}),\quad\omega_{k}=\sum_{j=0}^{n-1}\xi_{n}^{j}\,\tilde{\Delta}^{n-1-j}\,\tilde{c}^{j}\,k^{2p+1}.$
(10)

Unfortunately this simple trick fails beyond the leading order in
$1/\tilde{c}$. At $1/\tilde{c}^{2}$ order one has to fix both constant and
linear in $n_{k}$ terms, while simple $I_{k}\rightarrow(n_{k}+1/2)/\tilde{c}$
substitution leads to incorrect results.

We propose and verify up to $1/\tilde{c}^{2}$ order that the subleading
‘‘quantum correction’’ terms can be unambiguously fixed starting from the
analytic expression in terms of $1/\tilde{c}$ pertubative series of the
eigenvalues ${\sf Q}^{0}_{2n-1}$ of $Q_{2n-1}$ acting on the primary state
$|\Delta\rangle$. For leading $1/\tilde{c}$ term this statement is trivial –
taking all $n_{k}=0$ yields the constant term, which is simply leading
$1/\tilde{c}$ term in ${\sf Q}^{0}_{2n-1}$. At the $1/\tilde{c}^{2}$ order
this statement is more nuanced: naively ${\sf Q}^{0}_{2n-1}$ only fixes the
constant term with all $n_{k}=0$, but we show linear in $n_{k}$ terms can be
also fixed starting from ${\sf Q}^{0}_{2n-1}$. As a result we obtain spectrum
of $Q_{2n-1}$ up to first three orders in $1/c$ expansion, including the
leading $\Delta^{n}$ term. We then apply the obtained result to evaluate
thermal expectation values of $Q_{2n-1}$, free energy of the KdV Generalized
Gibbs Ensemble, and the asymptotic expansion of the quantum transfer matrix
acting on a primary state $|\Delta\rangle$, all at first few leading orders in
$1/c$.

The paper is organized as follows. In section 2 we discuss classic completely
integrable system associated with the finite zone potentials and evaluate
$Q_{2n-1}(h,I_{k})$ as a perturbative series in $I_{k}$. In section 3 we
discuss analytic form of $Q_{2k-1}$ acting on primary states. These two pieces
are combined in section 4 where we employ semiclassical quantization to obtain
the spectrum of qKdV charges in the first three orders of $1/\tilde{c}$
expansion. We also perform consistency checks, confirming our result. Section
5 is devoted to applications of the obtained result. In section 5.1 we
calculate thermal expectation values of $Q_{2n-1}$ and fix two leading orders
in $1/c$ of the associated differential operator ${\mathcal{D}}_{n}$

$\displaystyle{\rm
Tr}_{\Delta}(Q_{2k-1}\,q^{Q_{1}})={\mathcal{D}}_{n}\chi_{\Delta},\qquad\chi_{\Delta}\equiv{\rm
Tr}_{\Delta}(q^{Q_{1}}).$ (11)

In section 5.2 we discussed KdV Generalized Gibbs Ensemble and calculate its
free energy $-\ln Z_{\rm GGE}$,

$\displaystyle Z_{\rm GGE}={\rm Tr}\,e^{-\sum_{n}\mu_{2n-1}\,Q_{2n-1}},$ (12)

at leading order in $1/c$. In section 5.3 we use the asymptotic expansion to
calculate the quantum transfer matrix acting on a primary state at first two
orders in $1/c$ expansion. We use analytic continuation to extend the validity
beyond the asymptotic regime, but notice that certain non-pertubative terms
are missing. We conclude with a discussion in section 6. The paper also
includes a number of appendices. Appendix A provides an easy derivation of (1)
by quantizing boundary gravitons of semiclassical gravity in AdS3. Appendix B
evaluates $Q_{2n-1}(h,I_{k})$ at first two orders in $I_{k}$ by explicitly
introducing normal coordinates at the origin of the co-adjoint orbit of the
Virasoro algebra. Appendix C provides technical details concerning Novikov’s
one-zone potentials. Appendix D develops the technique of dealing with the
multi-zone potentials in the limit of the infinitesimally small zones.
Finally, appendix E provides the details of calculating the spectrum of
$Q_{2n-1}$ acting on primary states based on ODE/IM correspondence.

##  2 Calculation of $Q_{2n-1}(h,I_{k})$

In this section our goal is to find expression for $Q_{2n-1}$ in terms of the
orbit invariant $h$ and action variables $I_{k}$, by expanding pertubatively
up to cubic order in $I_{k}$,

$\displaystyle Q_{2n-1}$
$\displaystyle=h^{n}+\sum_{k}f^{(n,1)}_{k}I_{k}+f^{(n,2)}_{k}I_{k}^{2}+f^{(n,3)}_{k}I_{k}^{3}+\sum_{\begin{subarray}{c}k<\ell\end{subarray}}f^{(n,2)}_{k,\ell}I_{k}I_{\ell}+$
(13)
$\displaystyle\qquad\quad\,\,\,\sum_{k\neq\ell}f^{(n,3)}_{k,\ell}I_{k}^{2}I_{\ell}+\sum_{k<\ell<p}f^{(n,3)}_{k,\ell,p}I_{k}I_{\ell}I_{p}+\mathcal{O}(I^{4}).$

Coefficients $f$ are $h$-dependent. First three
$f^{(n,1)}_{k},f^{(n,2)}_{k},f^{(n,3)}_{k}$ will be found using one-zone
potentials in section 2.2. Using two-zone potentials we will find
$f^{(n,2)}_{k,\ell}$ and $f^{(n,3)}_{k,\ell}$ in section 2.3, while
coefficient $f^{(n,3)}_{k,\ell,p}$ will be fixed using three-zone potentials
in section 2.4. An alternative brute-force derivation of (13) up to quadratic
order in $I_{k}$ is given in the appendix B.

###  2.1 Finite zone potentials: an introduction

The starting point is the ‘‘Schr$\ddot{\rm o}$dinger’’ equation

$\displaystyle-\psi^{\prime\prime}+{u\over 4}\,\psi=\lambda\,\psi,$ (14)

with the periodic real-valued potential $u(\varphi+2\pi)=u(\varphi)$. For any
real $\lambda$ there are two linearly-independent quasi-periodic solutions

$\displaystyle\psi_{\pm}(\varphi+2\pi)=e^{\pm 2\pi
i\,p(\lambda)}\psi_{\pm}(\varphi).$ (15)

Here quasi-momentum $p(\lambda)$ could be either real or pure imaginary.
Values of $\lambda\in{\mathbb{R}}$ for which $p(\lambda)$ is imaginary are
called ‘‘forbidden zone.’’ At the end of forbidden zones $p(\lambda)$ is
integer or half-integer such that $\psi_{\pm}$ become periodic or antiperiodic
and linearly dependent. Normally, for such $\lambda$, another linearly
independent singular solution appears. Yet occasionally there are two linearly
independent regular periodic or antiperiodic solutions for the same $\lambda$.
In this case forbidden zone degenerates and disappears, with $p(\lambda)$
being real everyone in the vicinity of that point. We provide examples below.

A general potential $u$ would have an infinite number of forbidden zones, but
there are special classes when only a finite number of forbidden zones are
non-degenerate, Such $u$ are called finite zone potentials. They were
introduced in a famous work novikov1974periodic  and often refereed to as
Novikov potentials.

### Example: zero zone potential

Let us consider a constant potential $u=4\lambda_{0}=Q_{1}$ with some real
$Q_{1}$. A solution to (14) can be readily found

$\displaystyle\psi_{\pm}(\varphi)=e^{\pm ip(\lambda)\varphi},\qquad
p(\lambda)=\sqrt{\lambda-\lambda_{0}}.$ (16)

For any $\lambda>Q_{1}/4$ quasi-potential is real, i.e. there are no forbidden
zones, except for $\lambda\in(-\infty;Q_{1}/4)$. The solutions (16) are
linearly independent, including $\lambda=(Q_{1}+k^{2})/4$ for natural $k$,
when $\psi_{\pm}$ are (anti)periodic. Values $\lambda=(Q_{1}+k^{2})/4$ mark
the ends of degenerate forbidden zones.

### Example: ‘‘opening’’ a zone

Let us now consider the potential
$u=Q_{1}+\epsilon\cos(k\varphi)+O(\epsilon^{2})$ where $Q_{1}$ is a constant,
$k$ is positive integer, and $\epsilon$ is some infinitesimal parameter. Using
quantum mechanics perturbation theory we find at leading order that all
eigenvalues of periodic and anti-periodic problems remain the same and double-
degenerate, except for $\lambda_{k}$ which splits into

$\displaystyle\lambda_{k}^{\pm}={Q_{1}+k^{2}\over 4}\pm{\epsilon\over 2}.$
(17)

Hence now there are two forbidden zones, $(-\infty,Q_{1}/4)$ and
$(\lambda_{k}^{-},\lambda_{k}^{+})$.

Finite-zone potentials are characterized by the ends of non-degenerate zones
$\lambda_{i}$. For the zero-zone potential above there is only one parameter
$\lambda_{0}=Q_{1}/4$. After one zone is opened, there are three parameters:
‘‘energy’’ of the ground state $\lambda_{0}$, $\lambda_{1}=\lambda_{k}^{-}$
and $\lambda_{2}=\lambda_{k}^{+}$. In general an $m$-zone potential is
characterized by

$\displaystyle\lambda_{0}<\lambda_{1}<\dots<\lambda_{2m},$ (18)

with the forbidden zones $(-\infty,\lambda_{0})$ and
$(\lambda_{2i-1},\lambda_{2i})$, $i=\overline{1,m}$. For each set
$\\{\lambda_{i}\\}$ we can define a hyperelliptic curve

$\displaystyle y^{2}=\prod_{i=0}^{2m}(\lambda-\lambda_{i}),$ (19)

while the quasi-momentum $p$ being fixed in terms of its differential

$\displaystyle dp={\lambda^{m}+r_{m-1}\lambda^{n-1}+\dots r_{0}\over
2\,y}d\lambda,\qquad p(\lambda_{0})=0.$ (20)

The latter is defined in such a way that the integrals of $dp$ over $a$-cycles
vanish

$\displaystyle\oint\limits_{a_{i}}dp=2\int_{\lambda_{2i-1}}^{\lambda_{2i}}dp=0.$
(21)

This fixes $m$ coefficients $r_{0},\dots,r_{m-1}$. Furthermore for the
potential associated with $\\{\lambda_{i}\\}$ to be $2\pi$-periodic we must
additionally require integrals over $b$-cycles

$\displaystyle
w_{i}=\oint_{b_{i}}dp=2\int_{\lambda_{2i-2}}^{\lambda_{2i-1}}dp$ (22)

to be integer-valued

$\displaystyle w_{i}=k_{i}-k_{i-1}.$ (23)

Here natural $k_{i}$ satisfying $k_{i+1}>k_{i}$, $k_{0}\equiv 0$, label opened
zones. These are additional $m$ constrains, which reduce the total number of
independent parameters $\lambda_{i}$ to $m+1$.

A given set $\\{\lambda_{i}\\}$ which satisfies (21,23), such that only $m+1$
parameters are independent, defines periodic potential $u(\varphi)$, but in a
non-unique way. Individual potentials are labeled by points of the Jacobian of
curve (19), with all of them sharing the same spectrum. In other words
isospectral potentials form an $m$-dimensional torus, while full space of
$m$-zone potentials is therefore $2m+1$ dimensional.

At this point we would like to make a connection with the Virasoro algebra.
Consider Hill’s equation, which is ‘‘Schrodinger’’ equation (14) with
$\lambda=0$,

$\displaystyle-\psi^{\prime\prime}+{u\over 4}\,\psi=0.$ (24)

One can re-parametrize the circle going from $\varphi$ to
$\tilde{\varphi}(\varphi)$ such that
$\tilde{\varphi}(\varphi+2\pi)=\tilde{\varphi}(\varphi)+2\pi$. Then wave-
function and the potential also change accordingly

$\displaystyle\tilde{\psi}(\tilde{\varphi})$ $\displaystyle=$
$\displaystyle\psi(\varphi)\left({d\tilde{\varphi}\over
d\varphi}\right)^{-1/2},$ (25) $\displaystyle\tilde{u}(\tilde{\varphi})$
$\displaystyle=$ $\displaystyle\left({d\tilde{\varphi}\over
d\varphi}\right)^{-2}\left(u+2(S\tilde{\varphi})(\varphi)\right),$ (26)

where Schwarzian derivative

$\displaystyle(S\theta)(\varphi)\equiv{\theta^{\prime\prime\prime}\over\theta^{\prime}}-{3\over
2}\left({\theta^{\prime\prime}\over\theta^{\prime}}\right)^{2}.$ (27)

From (26) it is clear that $u$ is an element from the co-adjoint orbit of
Virasoro algebra with the Schwarzian derivative term appearing because of
central extension Witten . All potentials $u(\varphi)$ related by circle
reparametrizations, i.e. belonging to the same co-adjoint orbit share the same
invariant – quasi-momentum at zero,

$\displaystyle\psi(2\pi)/\psi(0)=e^{2\pi ip(0)},$ (28)

which is evident from (25). In other words

$\displaystyle-4p(0)^{2}=h$ (29)

is the invariant of $u$ characterizing the orbit itself. By choosing an
appropriate $\tilde{\varphi}$ the potential always333An implicit assumption
here is that $u$ belongs to the regular orbit ${\rm
diff}(\mathbb{S}^{1})/\mathbb{S}^{1}$, which upon quantization, becomes Verma
module.  can be brought to a constant form, in which case

$\displaystyle\tilde{u}(\tilde{\varphi})=h.$ (30)

The co-adjoint orbit is a symplectic space equipped with the Kirillov-Kostant-
Souriau bracket

$\displaystyle{c\over
24}\\{u(\varphi),u(\varphi^{\prime})\\}=-2\pi{\mathcal{D}}\delta(\varphi-\varphi^{\prime}),\qquad{\mathcal{D}}=\partial
u+2u\partial-2\partial^{3}.$ (31)

Here, using linearity of symplectic form we introduce a formal parameter $c$,
which later will be identified with the CFT central charge. Any Hamiltonian
flow defined by (31) leaves $h$ invariant.

There is an infinite tower of the so-called KdV Hamiltonians $Q_{2k-1}$, which
can be defined recursively with help fo Gelfand-Dikii polynomials $R_{n}$,

$\displaystyle Q_{2n-1}={1\over
2\pi}\int_{0}^{2\pi}R_{n}d\varphi\,\qquad\qquad\partial R_{n+1}={n+1\over
2n+1}{\mathcal{D}}R_{n},\quad$ (32) $\displaystyle R_{0}=1,\quad R_{1}=u,\quad
R_{2}=u^{2}-{4\over 3}\partial^{2}u,\quad
R_{3}=u^{3}-4u\partial^{2}u-2(\partial u)^{2}+{8\over 5}\partial^{4}u,\dots$

Their Hamiltonian flows generate isospectral deformations of $u$

$\displaystyle\delta u={c\over 24}\\{Q_{2n-1},u\\}=(2n-1)\partial R_{n},$ (33)

while they all remain in involution $\\{Q_{2n-1},Q_{2\ell-1}\\}=0$.

We now consider a space of all $m$-zone potentials sharing the same $h$. This
is a $2m$-dimensional subspace within the orbit parametrized by $h$, which we
will denote as ${\mathcal{F}}_{m}(h)$. The pullback of the symplectic form on
this space is non-degenerate, hence it is also a symplectic manifold equipped
with the Poisson bracket. Isospectral flows leave this manifold invariant.
Upon restricting to ${\mathcal{F}}_{m}(h)$, only first $n$ KdV Hamiltonians
remain algebraically independent. The flows they generate move $u$ along the
Jacobian of (19), which is the Liouvillian torus of a completely integrable
dynamical system defined by $Q_{2n-1}$, $n\leq m$. In other words the geometry
of ${\mathcal{F}}_{m}(h)$ is a $m$-dimensional torus parametrized by angle
variables fibered above a base parametrized by $m$ variables $Q_{2n-1}$.
Alternatively, one can introduce $m$ action variables $I_{k}$ parameterizing
the base and forming canonical conjugate pairs with angle variables.

In terms of $dp$ (20) values of KdV charges are given by an expansion at
infinity

$\displaystyle Q_{2n-1}={2\Gamma(n+1)\Gamma(1/2)\over\Gamma(n+1/2)}{4^{n}\over
2\pi i}\oint\limits_{\infty}dp\,\lambda^{n-1/2},$ (34)

while the action variables are

$\displaystyle
I_{k}={i\over\pi}\oint\limits_{a_{k}}p{d\lambda\over\lambda}={1\over
i\pi}\oint\limits_{a_{k}}dp\ln\lambda.$ (35)

Functional dependence of $Q_{2n-1}$ for $n>m$ on the first $m$ ones readily
follows from (34) and the form of $dp$ (20).

Our task is conceptually trivial: we want to learn an explicit change of
variables on the base of ${\mathcal{F}}_{m}(h)$ from $Q_{2n-1}$ to $I_{k}$.
The expressions for $Q_{2n-1}(h,I_{k})$ is not available in the closed form,
we therefore will find first few orders by expanding it in powers of $I_{k}$.
There is one notable exception, using Riemann bilinear relation with two one-
forms $dp$ and $pd\lambda/\lambda$ one can show in full generality

$\displaystyle Q_{1}=h+\sum_{k}k\,I_{k}.$ (36)

Our main approach will be based on parameterizing both $Q_{2n-1}$ and $I_{k}$
in terms of the spectral curve $i=\overline{0,m}$, with the infinitesimal
$\lambda_{2i}-\lambda_{2i-1}$, and then re-expressing $Q_{2n-1}$ in terms of
$I_{k}$. There is an alternative straightforward approach, to parametrize the
potential $u(\varphi)$ in terms of its Fourier modes $u_{\ell}$, and then
express both $Q_{2n-1}$ and $I_{k}$ in terms of $u_{\ell}$. We develop this
method in the appendix B and confirm the expansion (13) up to second order in
$I_{k}$.

###  2.2 One-zone potentials

Before we consider one-zone potential in detail, we revisit the zero-zone
potential $u=Q_{1}\equiv 4\lambda_{0}$ and readily find differential

$\displaystyle dp={d\lambda\over 2\sqrt{\lambda-\lambda_{0}}}$ (37)

to be defined on a Riemann sphere. This is the simplest possible case. In this
case $p=\sqrt{\lambda-\lambda_{0}}$, $u(\varphi)=h=4\lambda_{0}$ and the whole
symplectic space ${\mathcal{F}}_{0}(h)$ shrinks to a point. All KdV
Hamitonians are fixed by $h$, $Q_{2n-1}=h^{n}$ with all action variables
identically equal to zero.

Next, we consider the differential

$\displaystyle dp={(\lambda-r)d\lambda\over
2\sqrt{(\lambda-\lambda_{0})(\lambda-\lambda_{1})(\lambda-\lambda_{2})}}$ (38)

parameterized by $\lambda_{i},r_{0}$. It is defined on a torus – a Riemann
curve of genus one. We assume that $(\lambda_{2},\lambda_{1})$ correspond to
$k$-th zone. After satisfying (21) and (23), which requires evaluating
elliptic integrals, we find one-parametric family

$\displaystyle\lambda_{2}=\lambda_{0}+{k^{2}\over
4}\theta_{3}(\tau)^{4},\quad\lambda_{1}=\lambda_{0}+{k^{2}\over
4}\theta_{4}(\tau)^{4},\quad r=\lambda_{0}+{k^{2}\over
4}\theta_{4}(\tau)^{4}\left(1+2{\partial\ln\theta_{3}^{2}(\tau)\over\partial\ln
m}\right),$ (39)

where $m=\theta_{2}^{4}(\tau)/\theta_{3}^{4}(\tau)$ and $\tau=i\tau_{2}$ with
positive $\tau_{2}$. In what follows we use444Our definition of $q$ is aligned
with Wolfram Mathematica. In this section $q$ denotes modular parameter of the
genus one elliptic curve $y(\lambda)$. In section 5.1 we use $q$ to denote
modular parameter of the CFT spacetime torus. $q=e^{i\pi\tau}$ such that
$\theta_{2}=\sum_{n}q^{(n+1/2)^{2}}$, $\theta_{3}=\sum_{n}q^{n^{2}}$,
$\theta_{4}=\sum_{n}(-1)^{n}q^{n^{2}}$.

To impose the orbit constraint $-4p(0)^{2}=h$ it is more convenient to use the
following trick. First we evaluate

$\displaystyle Q_{1}=4(\lambda_{0}+\lambda_{1}+\lambda_{2})-8r,$ (40)

which expresses $\lambda_{0}$ in terms of $Q_{1}$ and $q$ expansion,

$\displaystyle
4\lambda_{0}=Q_{1}-k^{2}\left(\theta_{2}^{4}-4\theta_{4}^{4}{\partial\ln\theta_{3}^{2}(\tau)\over\partial\ln
m}\right)=Q_{1}-32k^{2}q^{2}\left(1+2q^{2}+4q^{4}+4q^{6}+\dots\right),\qquad$
(41)

and then use (35) to evaluate action variable perturbatively in $q$,

$\displaystyle I_{k}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int^{\lambda_{2}}_{\lambda_{1}}\frac{d\lambda(\lambda-r)\log\lambda}{\sqrt{(\lambda-\lambda_{0})(\lambda-\lambda_{1})(\lambda_{2}-\lambda)}}$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{2(-1)^{n}(\lambda_{1}-\lambda_{0})^{n+1}}{n\sqrt{\lambda_{2}-\lambda_{0}}\lambda_{0}^{n}}\left[\frac{F\left(\frac{3}{2},\frac{1}{2},1;m\right)}{F\left(\frac{1}{2},\frac{1}{2},1;m\right)}F\left(n+\frac{1}{2},\frac{1}{2},1;m\right)-F\left(n+\frac{3}{2},\frac{1}{2},1;m\right)\right].$

Here $F\equiv{}_{2}F_{1}$ is the hypergeometric function such that
$F\left(\frac{3}{2},\frac{1}{2},1;m\right)=\theta_{3}^{2}$.

An infinite sum over $n$ above has to be evaluated individually for each term
in $q$ expansion. This gives $I_{k}$ as a function of $\lambda_{0}$ and $q$,
$I_{k}={32k^{3}q^{2}\over k^{2}+4\lambda_{0}}+O(q^{4})$, which with help of
(41) can be expressed as a function of $Q_{1}$ and $q$,

$\displaystyle I_{k}$ $\displaystyle=$
$\displaystyle\frac{32k^{3}}{k^{2}+Q_{1}}q^{2}+\frac{64k^{3}\left(17k^{4}+12k^{2}Q_{1}+3Q_{1}^{2}\right)}{\left(k^{2}+Q_{1}\right)^{3}}q^{4}$
$\displaystyle+\frac{128k^{3}\left(5k^{2}+Q_{1}\right)\left(77k^{6}+69k^{4}Q_{1}+27k^{2}Q_{1}^{2}+3Q_{1}^{3}\right)}{\left(k^{2}+Q_{1}\right)^{5}}q^{6}+\mathcal{O}(q^{8}).$

At this point we use (36), which is exact, $Q_{1}=h+kI_{k}.$. Using $I_{k}$
given as a $q$-series expansion with $Q_{1}$-dependent coefficients (2.2),
with help of (2.2) we express $Q_{1}$ as a series in $q$ with $h$-dependent
coefficients by iteratively substituting $Q_{1}$ written as an $h$-dependent
series in $q$. Once we find $Q_{1}=Q_{1}(h,q)$, $I_{k}$ can be deduced from
(2.2),

$\displaystyle I_{k}$ $\displaystyle=$
$\displaystyle\frac{32k^{3}}{h+k^{2}}q^{2}+\frac{64\left(3h^{2}k^{3}+12hk^{5}+k^{7}\right)}{\left(h+k^{2}\right)^{3}}q^{4}$
(45)
$\displaystyle+\frac{128k^{3}\left(3h^{4}+42h^{3}k^{2}+108h^{2}k^{4}-58hk^{6}+k^{8}\right)}{\left(h+k^{2}\right)^{5}}q^{6}+\mathcal{O}(q^{8}).$

At this point it is straightforward to re-express $q$ as a $h$-dependent power
series in $I_{k}$, $q^{2}={h+k^{2}\over 32k^{3}}I_{k}+O(I_{k}^{2})$.

To obtain coefficients $f^{(n,i)}$ (13) we act as follows. From the definition
(34) we can find $Q_{2n-1}$ as a polynomial in $\lambda_{i}$ and $r$. Using
expressions for $\lambda_{i},r$ (39) and (41), where $Q_{1}$ is understood as
a function of $h,q$ we write $Q_{2n-1}$ as an $h$-dependent power series in
$q$. After that it is straightforward to use $q^{2}=q^{2}(h,I_{k})$ to re-
express $Q_{2n-1}$ as an $h$-dependent power series in $I_{k}$,

$\displaystyle Q_{2n-1}$
$\displaystyle=h^{n}+f^{(n,1)}_{k}I_{k}+f^{(n,2)}_{k}I_{k}^{2}+f^{(n,3)}_{k}I_{k}^{2}+\mathcal{O}(I_{k}^{3}),$
(46)

thus fixing $f^{(n,i)}$,

$\displaystyle f^{(n,1)}_{k}$
$\displaystyle=\sum_{j=0}^{n-1}\frac{\sqrt{\pi}(2n-1)\Gamma(n+1)}{2\Gamma(j+\frac{3}{2})\Gamma(n-j)}h^{n-1-j}k^{2j+1}=\sum_{j=0}^{n-1}\xi_{n}^{j}\,h^{n-1-j}k^{2j+1},$
(47) $\displaystyle f^{(n,2)}_{k}$
$\displaystyle=\sum_{j=0}^{n-1}\frac{\sqrt{\pi}(2n-1)\Gamma(n+1)(j(2n+1)-2n+2)}{16\,\Gamma\left(j+\frac{3}{2}\right)\Gamma(n-j)}h^{n-1-j}k^{2j},$
(48) $\displaystyle f^{(n,3)}_{k}$ $\displaystyle={(2n-1)n(n-1)\over
64k^{3}}h^{n}+$ (49)
$\displaystyle\quad\sum_{j=0}^{n-1}\frac{\sqrt{\pi}(2n-1)\Gamma(n+1)\,{\sf
p}}{1536\,\Gamma\left(j+\frac{5}{2}\right)\Gamma(n-j)}h^{n-1-j}k^{2j-1},$
$\displaystyle{\sf p}$
$\displaystyle=4j^{3}(2n+1)(2n+3)-2j^{2}(2n+1)(10n-21)-3j(2n+3)(10n-7)+36(n-1)(2n-1).$

More technical details concering one-zone calculation can be found in appendix
C.

###  2.3 Two-zone potentials

In case of two zones the differential

$\displaystyle dp={(\lambda-r_{1})(\lambda-r_{2})d\lambda\over
2\sqrt{\prod_{i=0}^{4}(\lambda-\lambda_{i})}}$ (50)

depends on seven parameters subject to 4 constraints (21) and (23).
Corresponding integrals can not be evaluated analytically. We therefore
proceed by expanding perturbatively, assuming both zones, and hence
corresponding action variables, are small. We introduce two infinitesimal
variables $\epsilon_{1},\epsilon_{2}$ of the same order, such that
$\lambda_{2}-\lambda_{1}$ is of order $\epsilon_{1}$ and
$\lambda_{4}-\lambda_{3}$ is of order $\epsilon_{2}$. Action variables are
quadratic in $\epsilon_{i}$,
$I_{k}\sim\epsilon_{1}^{2},I_{\ell}\sim\epsilon_{2}^{2}$, where we assumed
$(\lambda_{1},\lambda_{2})$ and $(\lambda_{3},\lambda_{4})$ correspond to
$k$-th and $\ell$-zones respectively. Our goal is to find $Q_{2n-1}$ up to
third order in the pertubative expansion in $I_{k},I_{\ell}$. Hence in what
follows we must expand all quantities in $\epsilon_{i}$ up to sixth order. The
details of this calculation can be found in in Appendix D.

After satisfying (21) and (23) we find $\lambda_{i}$ for $i\geq 1$ and $r_{i}$
in terms of $\lambda_{0}$ and $\epsilon_{1},\epsilon_{2}$, as a perturbative
expansion in $\epsilon_{i}$. Then, we evaluate $I_{k},h$ and $Q_{2n-1}$ also
as function of $\lambda_{0}$ and $\epsilon_{1},\epsilon_{2}$, similarly
expanding in $\epsilon_{i}$ up to and including sixth order. By matching both
sides of (13) we find coefficients $f^{(m,n)}_{k,\ell}$, yielding

$\displaystyle f^{(n,2)}_{k,\ell}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n-1}\frac{\sqrt{\pi}(2n-1)^{2}\Gamma(n+1)}{4\Gamma(n-j)\Gamma\left(j+\frac{3}{2}\right)}h^{n-1-j}\sum_{s=0}^{j-1}k^{2(j-s)-1}\ell^{2s+1},$
(51) $\displaystyle f^{(n,3)}_{k,\ell}$ $\displaystyle=$
$\displaystyle{\ell\over(k^{2}-\ell^{2})^{2}}\left(-{(2n-1)n(n-1)\over
4}h^{n}+\sum_{j=0}^{n-1}{\sqrt{\pi}(2n-1)^{2}\Gamma(n+1)\over
64\Gamma(n-j)\Gamma\left(j+\frac{5}{2}\right)}h^{n-1-j}\,{\sf q}\right),$ (52)
$\displaystyle{\sf q}$ $\displaystyle=$
$\displaystyle-4(2n+1){k^{2j+4}-\ell^{2j+4}\over
k^{2}-\ell^{2}}+k^{2j+2}(3+2j)(j(2n+1)-4n+5)+\ell^{2j+2}2(3j+n+5)+$
$\displaystyle k^{2j}\ell^{2}(3+2j)(j(2n+1)-2n+2)+k^{2}\ell^{2j}(3+2j)(4n-1).$

###  2.4 Three-zone potentials

Extending calculations of the previous section using the technique of appendix
D to the three-zone case we can fix

$\displaystyle
f^{(n,3)}_{k,\ell,p}=\sum_{j=0}^{n-3}{\sqrt{\pi}(2n-1)^{3}\Gamma(n+1)(n-2-j)\over
8\Gamma(n-1-j)\Gamma\left(j+{7\over
2}\right)}h^{n-3-j}\sum_{s_{1}=0}^{j}\sum_{s_{2}=0}^{j-s_{1}}k^{2j+1-2(s_{1}+s_{2})}\ell^{2s_{1}+1}p^{2s_{2}+1}.$
(53)

###  2.5 Consistency check

In case of an $m$-zone potential we can parametrize the differential $dp$ with
help of $\lambda_{0}$ and $\epsilon_{i}$, $1\leq i\leq m$, cf. (170-175),

$\displaystyle\lambda_{i}$ $\displaystyle=$
$\displaystyle\lambda_{0}+\dots,\qquad 1\leq i\leq 2m,$ (54) $\displaystyle
r_{i}$ $\displaystyle=$ $\displaystyle\lambda_{0}+\dots,\qquad 1\leq i\leq m,$
(55)

where dots stand for $\epsilon_{i}$ but not $\lambda_{0}$-dependent terms.
Similarly action variables $I_{k}$, charges $Q_{2n-1}$ and the orbit parameter
$h=-4p(0)^{2}$ will be some functions of $\lambda_{0}$ and $\epsilon_{i}$.
While dependence of $I_{k}$ and $h$ on $\lambda_{0}$ is non-trivial, since
$Q_{2n-1}$ are the coefficients of $1/\lambda$ expansion of $p(\lambda)$ at
infinity and $\lambda_{0}$ is simply the shift of the argument of
$p(\lambda)$, we find

$\displaystyle
Q_{2n-1}=\sum_{k=0}^{n}\frac{\Gamma(n+1)}{\Gamma(k+1)\Gamma(n-k+1)}(4\lambda_{0})^{n-k}Q_{2k-1}^{0}.$
(56)

Here $Q_{2k-1}^{0}$ are the charges evaluated with help of (34) taking
$\lambda_{0}=0$ in (54,55). Assuming we know $Q_{2n-1}(h,I_{k})$ where
$h=h(\lambda_{0},\epsilon_{i})$ and $I_{k}=I_{k}(\lambda_{0},\epsilon_{i})$,
one can introduce $I_{k}^{0}=I_{k}(0,\epsilon_{i})$ such that
$Q_{2k-1}^{0}=Q_{2k-1}(0,I_{k}^{0})$. Here first argument is zero simply
because $h(0,\epsilon_{i})=0$. Then both sides of equation (56) become
functions of $\lambda_{0}$ and $\epsilon_{i}$, providing a non-trivial check.

There is an alternative way to use (56) to check the consistency of the
perturbative expansion(13) with the coefficients found in the text. We can
invert $h=h(\lambda_{0},\epsilon_{i})$ and
$I_{k}=I_{k}(\lambda_{0},\epsilon_{i})$ to express both $\lambda_{0}$ and
$I_{k}^{0}$ via $h$ and $I_{k}$,

$\displaystyle\lambda_{0}$ $\displaystyle=$ $\displaystyle
h+\sum_{k}-{hI_{k}\over
k}+\frac{h\left(h+5k^{2}\right)I_{k}^{2}}{8k^{4}}-\frac{h\left(5h^{2}+30hk^{2}+41k^{4}\right)I_{k}^{3}}{128k^{7}}+\sum_{k<\ell}{hI_{k}I_{\ell}\over
k\ell}+$
$\displaystyle\sum_{k\neq\ell}\frac{hI_{k}^{2}I_{\ell}\left(h^{2}\ell^{2}-h\left(k^{4}-4k^{2}\ell^{2}+\ell^{4}\right)-5k^{6}+11k^{4}l^{2}-5k^{2}\ell^{4}\right)}{8k^{4}\ell(k-\ell)^{2}(k+\ell)^{2}}-\sum_{k<\ell<p}\frac{hI_{k}I_{\ell}I_{p}}{k\ell
p}+{\mathcal{O}}(I^{4}),$ $\displaystyle I_{k}^{0}$ $\displaystyle=$
$\displaystyle I_{k}+{hI_{k}\over
k^{2}}-\frac{hI_{k}^{2}\left(h+5k^{2}\right)}{8k^{5}}+\frac{hI_{k}^{3}\left(5h^{2}+30hk^{2}+41k^{4}\right)}{128k^{8}}-\sum_{\ell\neq
k}\frac{hI_{k}I_{\ell}\left(h+k^{2}\right)}{k^{2}\ell\left(k^{2}-\ell^{2}\right)}+$
$\displaystyle\sum_{\ell\neq
k}\frac{hI_{k}^{2}I_{\ell}\left(2h^{2}\left(-k^{4}+2k^{2}\ell^{2}+\ell^{4}\right)+hk^{2}\left(7k^{4}-14k^{2}\ell^{2}+15\ell^{4}\right)+k^{4}\left(5k^{4}-10k^{2}\ell^{2}+9\ell^{4}\right)\right)}{8k^{5}l\left(k^{2}-\ell^{2}\right)^{3}}+$
$\displaystyle\sum_{\ell\neq
k}\frac{hI_{k}I_{\ell}^{2}\left(h^{2}\left(k^{4}+5k^{2}\ell^{2}-2\ell^{4}\right)+h\left(k^{6}+10k^{4}\ell^{2}-9k^{2}\ell^{4}+6l^{6}\right)+k^{2}\ell^{2}\left(5k^{4}-7k^{2}\ell^{2}+6l^{4}\right)\right)}{8k^{2}\ell^{4}\left(k^{2}-\ell^{2}\right)^{3}}+$
$\displaystyle\sum_{p\neq\ell\neq
k}\frac{hI_{k}I_{\ell}I_{p}\left(2h^{2}+3hk^{2}+k^{4}\right)}{k^{2}\ell
p\left(k^{2}-\ell^{2}\right)\left(k^{2}-p^{2}\right)}+{\mathcal{O}}(I^{4}).$

Now $Q_{2n-1}^{0}(0,I_{k}^{0}(h,I_{k}))$ is a function of $h,I_{k}$ and (56)
provides a non-trivial check for the coefficients in (13).

This check also ensures that $Q_{2n-1}(h,I_{k})$ satisfy another identity

$\displaystyle{1\over n+1}{\partial Q_{2n+1}\over\partial u_{0}}=Q_{2n-1},$
(57)

which follows from the properties of Gelfand-Dikii polynomials (32). Here
$Q_{2n-1}[u(\varphi)]$ are understood as functionals of $u(\varphi)$ and the
derivative is with respect the zero Fourier mode of $u(\varphi)$, while all
other Fourier modes are kept fixed. The shift of $u_{0}$ with all other modes
intact is equivalent to a shift of the spectrum by a constant, hence

$\displaystyle\left({\partial\over\partial
u_{0}}\right)_{u_{\ell}}=4\left({\partial\over\partial\lambda_{0}}\right)_{\epsilon_{i}}.$
(58)

Then (57) follows immediately from the right-hand-side of (56).

For an $m$-zone potential, all higher KdV charges $Q_{2n-1}$ are some
functions of first $m+1$ charges. Thus for one-zone potentials
$Q_{5},Q_{7},\dots$ are functions of $Q_{1},Q_{3}$, see e.g. section 2.4 of
Dymarsky_2020  for details. For the three-zone potentials higher $Q_{2n-1}$
would depend on $Q_{1},Q_{3},Q_{5},Q_{7}$ In principle this provides
additional consistency check for (13). In practice the dependence is so
complicated, it doesn’t provide a useful check even for the one-zone case.

##  3 ‘‘Energies’’ of primary states via ODE/IM correspondence

In the previous section we found classical expression for $Q_{2n-1}$ in term
of action variables $I_{k}$ and the orbit invariant $h$. Following the
standard rules of semiclassical quantization $I_{k}$ should be promoted to an
integer quantum number, while $h$ will become the dimension of the highest
weight (primary) state $\Delta$, marking representation of the Virasoro
algebra. It is easy to see, this naive receipt fails already for the values of
$Q_{2n-1}$ on a primary state $|\Delta\rangle$. Indeed, taking all $I_{k}$ to
zero, we readily find $Q_{2n-1}=h^{k}$, which upon the naive quantization
yields ${\sf Q}^{0}_{2n-1}=\Delta^{n}$ where

$\displaystyle Q_{2n-1}|\Delta\rangle={\sf Q}^{0}_{2n-1}|\Delta\rangle.$ (59)

This answer is missing $c$-dependent terms. Explicit values of ${\sf
Q}^{0}_{2n-1}$ for $n\leq 8$ were calculated in bazhanov97zero  via brute-
force approach, using explicit expressions for $Q_{2n-1}$ in terms of free
field representation. The pattern is clear, while $\Delta^{n}$ is indeed the
leading term, full expression is a polynomial in both $\Delta$ and $c$ of
order $n$.

There is no known receipt to obtain exact ${\sf Q}^{0}_{2n-1}$ from the
semiclassical quantization, hence our strategy will be the following. We will
combine exact expression for ${\sf Q}^{0}_{2n-1}$ in the large $c$ limit,
which will be obtained in this section by a different method, with the
classical result of section 2, to find spectrum of excited states in the large
$c$ limit in next section.

To find ${\sf Q}^{0}_{2n-1}$ we use ODE/IM correspondence, initiated in
bazhanov2001spectral ; bazhanov2003higher  and more recently developed in
dorey2020geometric  (also see Conti:2021xzr ), which relates qKdV spectrum to
solutions of an auxiliary Schr$\ddot{\rm o}$dinger equation

$\displaystyle\partial_{x}^{2}\Psi(x)+\left(E-x^{2\alpha}-\frac{l(l+1)}{x^{2}}\right)\Psi(x)=0,$
(60)

where

$\displaystyle(l+1/2)^{2}=4(\alpha+1)\tilde{\Delta},\quad\tilde{c}=-\frac{\alpha^{2}}{4(\alpha+1)}.$
(61)

Equation (60) can be solved using WKB approximation by systematically
expanding in a small parameter. This leads to a quadratic ODE which can be
solved iteratively. We delegate all details to Appendix E and only write down
iterative relation which defines coefficients $c^{(n)}_{k}$ for $n\geq 1$,
$n\geq k\geq 0$,

$\displaystyle\sum_{j=0}^{n}\sum_{p=0}^{j}\sum_{q=0}^{n-j}\delta_{p+q,k}c_{p}^{(j)}c_{q}^{(n-j)}-2\left[n-k-u-\frac{n-2}{2\alpha}\right]c_{k-1}^{(n-1)}+\left(2k-3n+4\right)c_{k}^{(n-1)}=0,$
(62)

and we formally assumed $c^{(n)}_{-1}=c^{(n)}_{n+1}=0$,
$u^{2}=-\tilde{\Delta}/\tilde{c}$, and the starting values are

$\displaystyle c_{0}^{(0)}=-\frac{1}{\alpha},\qquad
c_{0}^{(1)}=-\frac{1}{2},\quad c_{1}^{(1)}=\frac{1}{2\alpha}-u.$ (63)

Coefficients $c_{k}^{(n)}$ determine values of $Q_{2n-1}$ acting on primaries
dorey2020geometric ,

$\displaystyle{\sf
Q}^{0}_{2n-1}=\frac{(2n-1)\Gamma(n+1)}{\sqrt{\pi}\Gamma(1-\frac{2n-1}{2\alpha})4^{n}(\alpha+1)^{n}}\sum_{k=0}^{2n}c_{k}^{(2n)}\Gamma\left(k+\frac{3}{2}-3n\right)\Gamma\left(2n-k-\frac{2n-1}{2\alpha}\right).$
(64)

Although this is not obvious, ${\sf Q}^{0}_{2n-1}$ given by (64) is a
polynomial in terms of $\tilde{\Delta}$ and $\tilde{c}$. After some algebra we
find leading order expansion

$\displaystyle{\sf Q}^{0}_{2n-1}$
$\displaystyle=\tilde{\Delta}^{n}+\sum_{j=0}^{n-1}\tilde{R}^{(1)}_{n,j}\tilde{\Delta}^{n-j-1}\,\tilde{c}^{j}+\sum_{j=0}^{n-2}\tilde{R}^{(2)}_{n,j}\tilde{\Delta}^{n-j-2}\,{\tilde{c}}^{j}+\sum_{j=0}^{n-3}\tilde{R}^{(3)}_{n,j}\tilde{\Delta}^{n-j-3}\,{\tilde{c}}^{j}+\mathcal{O}({\tilde{c}}^{n-3}).$
(65)

where

$\displaystyle\tilde{R}^{(1)}_{n,j}$
$\displaystyle=\frac{(2n-1)\sqrt{\pi}\Gamma(n+1)}{4\Gamma(j+\frac{3}{2})\Gamma(n-j)}\zeta(-2j-1)=\xi_{n}^{j}{\zeta(-2j-1)\over
2},$ (66) $\displaystyle\tilde{R}^{(2)}_{n,j}$
$\displaystyle=\frac{(2n-1)\sqrt{\pi}\Gamma(n+1)}{24\times
4\Gamma(j+\frac{5}{2})\Gamma(n-j-1)}\times$ (67)
$\displaystyle\Biggl{\\{}-6\zeta(-2j-3)\left(2j+3-(2n-1)y_{1}(j+1)\right)+3(2n-1)\zeta_{2}(j)\,\Biggr{\\}},$
$\displaystyle\tilde{R}^{(3)}_{n,j}$
$\displaystyle=\frac{(2n-1)\sqrt{\pi}\Gamma(n+1)}{24^{2}\times
4\Gamma(j+\frac{7}{2})\Gamma(n-j-2)}\Biggl{\\{}6^{2}\zeta(-2j-5)(2j^{2}+7j+5)-(2n-1)r_{n,j}\Biggr{\\}},$
$\displaystyle{r}_{n,j}$
$\displaystyle=12\zeta_{3}(j)+36\zeta_{2}(j+1)(y_{1}(j+2)+j+2)+3\left(4j^{2}+18j+23\right)\zeta(-2j-3)+$
$\displaystyle
36\zeta(-2j-5)\left(y_{1}^{2}(j+2)+2(j+2)y_{1}(j+2)+y_{2}(j+2)\right)+(2n+1)p_{j}.$
(68)

Functions $\zeta_{2},\zeta_{3},y_{1},y_{3}$ are defined in the Appendix E,
where we also give values of $p_{j}$ for $0\leq j\leq 17$.

##  4 Spectrum of quantum $Q_{2k-1}$

At this point we are ready to combine classical pertubative expression for
$Q_{2n-1}(h,I_{k})$ (13) with the ‘‘energies’’ of primary state (65) to obtain
${\sf Q}_{2n-1}$ up to first two non-trivial orders in $1/\tilde{c}$
expansion.

The naive semi-classical quantization would map the co-adjoint orbit invariant
$h$ and the actions variables $I_{k}$ on the classical side to dimension of
the primary state $\Delta$ and the excited state quantum numbers $n_{k}$
correspondingly,

$\displaystyle h\rightarrow{24\Delta\over c},\qquad
I_{k}\rightarrow{24n_{k}\over c}.$ (69)

Also classical charge $Q_{2n-1}$ should be rescaled by $(c/24)^{n}$. Starting
from (13) this correctly reproduces full quantum spectrum of $Q_{1}$ and the
leading $\Delta^{n}$ term in $Q_{2n-1}$. But it falls short of reproducing
sub-leading terms even for the primary state (65). The relation between
classical and quantum quantities (69) is only correct at the leading $c$
order. In GGE2  we observed that using $c-1$ as an expansion parameter leads
to more elegant expressions. This is confirmed by (65), which looks most
naturally if written in terms of $\tilde{\Delta}$ and $\tilde{c}$. We
therefore propose the following quantization map, which agrees with the naive
one at leading order,

$\displaystyle h\rightarrow{\tilde{\Delta}\over\tilde{c}},\qquad
I_{k}\rightarrow{n_{k}\over\tilde{c}},\qquad\tilde{\Delta}=\Delta-\tilde{c},\qquad\tilde{c}={c-1\over
24}.$ (70)

This does not solve the problem of reproducing subleadig terms in ${\sf
Q}_{2n-1}^{0}$, but this can be fixed, at least at first subleading order, by
introducing the Maslov index, $n_{k}\rightarrow\tilde{n}_{k}=n_{k}+1/2$. We
thus arrive at the following map,

$\displaystyle Q_{2n-1}(h,I_{k})\rightarrow{\sf
Q}_{2n-1}=\tilde{c}^{n}\,Q_{2n-1}(\tilde{\Delta}/\tilde{c},(n_{k}+1/2)/\tilde{c}).$
(71)

Infinite sums due to Maslov index contributing to ‘‘vacuum energy’’ should be
regularized using zeta-function regularization. It is now straightforward to
see that we immediately reproduce the leading $1/\tilde{c}$ term (66),

$\displaystyle
Q_{2n-1}=h^{n}+\sum_{k}f_{k}^{(n,1)}(h)\,I_{k}+O(I^{2})\rightarrow{\sf
Q}_{2n-1}=\tilde{\Delta}^{n}+\tilde{c}^{n-1}\sum_{k}f_{k}^{(n,1)}(\tilde{\Delta}/\tilde{c})\,\tilde{n}_{k}+O(\tilde{c}^{n-2})=$
$\displaystyle\tilde{\Delta}^{n}+\sum_{k}\sum_{j=0}^{n-1}\xi_{n}^{j}\,\tilde{\Delta}^{n-1-j}\tilde{c}^{j}k^{2j+1}(n_{k}+1/2)+O(\tilde{c}^{n-2})=$
$\displaystyle\tilde{\Delta}^{n}+\sum_{j=0}^{n-1}\xi_{n}^{j}\,\tilde{\Delta}^{n-1-j}\tilde{c}^{j}\left(\sum_{k}k^{2j+1}n_{k}+{\zeta(-2j-1)\over
2}\right)+O(\tilde{c}^{n-2})$ (72)

In other words, at first sub-leading order $\tilde{c}^{n-1}$ the quantization
prescription (71) leads to (1) which passes all available tests: matches the
spectrum of $Q_{1},Q_{3},Q_{5},Q_{7}$ (see section 4.1 below) and thermal
expectation values for $Q_{9},\dots,Q_{13}$ (see section 5.1 below) at the
order $\tilde{c}^{n-1}$.

There is another way to write (72). We can express ${\sf Q}_{2n-1}$ as ${\sf
Q}_{2n-1}^{0}$ plus the terms from the classical $Q_{2n-1}$ (13) which non-
trivially depend on $I_{k}$ using the substitution (70), i.e. without the
Maslov index,

$\displaystyle{\sf Q}_{2n-1}={\sf
Q}_{2n-1}^{0}+\tilde{c}^{n-1}\sum_{k}f_{k}^{(n,1)}(\tilde{\Delta}/\tilde{c})\,n_{k}+O(\tilde{c}^{n-2}).$
(73)

At $\tilde{c}^{n-1}$ order it is the same as (72).

To obtain the quantum spectrum at next order $\tilde{c}^{n-2}$, we could try
the prescription (71), apply the zeta-function regularization and notice that
many but not all terms from (67) are reproduced. Thus, we see that the
quantization (71) is exact only at leading $1/\tilde{c}$ order, at higher
orders the expression obtained from the classical $Q_{2n-1}$ has to be
modified as well. Indeed, starting from the classical (13) and using
substitution (70) we would find that terms contributing at the order
$\tilde{c}^{n-p}$ are homogeneous polynomials in $n_{k}$ of order $p$. This is
very restrictive and obviously incorrect. We already saw that even at the
first sub-leading order $\tilde{c}^{n-1}$ the homogeneous (linear) in $n_{k}$
terms have to be amended by a constant, i.e. $(n_{k})^{0}$ term. This suggest
the following ‘‘quantization rules’’: to obtain the quantum spectrum ${\sf
Q}_{2n-1}$ in $1/\tilde{c}$ expansion, one starts with the classical
perturbation expression (13) and make the substitution (70), together with the
overall rescaling by $\tilde{c}^{n}$. As the order $\tilde{c}^{n-p}$ this
fixes leading, homogeneous in $n_{k}$ terms of order $p$. These terms should
be amended by the sub-leading terms of order $p-1$, $p-2$, …, $0$ in $n_{k}$.
These terms should be regarded as quantum corrections and should be determined
separately, they do not follow from the classical answer in any simple way.
More explicitly,

$\displaystyle{\sf
Q}_{2n-1}=\tilde{\Delta}^{n}+\tilde{c}^{n-1}\left(\sum_{k}g^{(1)}_{k}n_{k}+g^{(1)}\right)+\tilde{c}^{n-2}\left(\sum_{k_{1},k_{2}}g^{(2)}_{k_{1},k_{2}}n_{k_{1}}n_{k_{2}}+\sum_{k}g^{(2)}_{k}n_{k}+g^{(2)}\right)+$
$\displaystyle\tilde{c}^{n-3}\left(\sum_{k_{1},k_{2},k_{3}}g^{(3)}_{k_{1},k_{2},k_{3}}n_{k_{1}}n_{k_{2}}n_{k_{3}}+\sum_{k_{1},k_{2}}g^{(3)}_{k_{1},k_{2}}n_{k_{1}}n_{k_{2}}+\sum_{k}g^{(3)}_{k}n_{k}+g^{(3)}\right)+\dots$
(74)

Here $g^{(p)}$ with different number of indexes denote different quantities.
The leading terms $g^{(p)}_{k_{1},\dots,k_{p}}$ are given by classical
expressions (13) upon the substitution (70)

$\displaystyle g^{(1)}_{k}$ $\displaystyle=$ $\displaystyle
f^{(n,1)}_{k}(\tilde{\Delta}/\tilde{c}),$ (75) $\displaystyle g^{(2)}_{k\ell}$
$\displaystyle=$ $\displaystyle{1\over 2}f^{(n,2)}_{k\ell},\quad
g^{(2)}_{kk}=f^{(n,2)}_{k},$ (76) $\displaystyle g^{(3)}_{k\ell m}$
$\displaystyle=$ $\displaystyle{1\over 6}f^{(n,3)}_{k\ell m},\quad
g^{(3)}_{kk\ell}={1\over 3}f^{(n,3)}_{k\ell},\quad
g^{(3)}_{kkk}=f^{(n,3)}_{k},$ (77)

for $k\neq\ell\neq m$ and $g^{(p)}$ are given by (66,67,68). This is
essentially the generalization of (73) to higher orders in $1/\tilde{c}$.
Coefficients $g^{(2)}_{k}$, $g^{(3)}_{k\ell},g^{(3)}_{k}$, etc. are quantum
corrections and a priory not known.

To fix $g^{(2)}_{k}$ we employ the following strategy, we will try to
‘‘salvage’’ the Maslov index quantization (71) by adding minimal possible
terms subleading in powers of $n_{k}$,

$\displaystyle{\sf
Q}_{2n-1}=\tilde{\Delta}^{n}+\tilde{c}^{n-1}\sum_{k}g^{(1)}_{k}\tilde{n}_{k}+\tilde{c}^{n-2}\left(\sum_{k_{1},k_{2}}g^{(2)}_{k_{1},k_{2}}\tilde{n}_{k_{1}}\tilde{n}_{k_{2}}+\sum_{k}\tilde{g}^{(2)}_{k}\tilde{n}_{k}+\tilde{g}^{(2)}\right)+\dots$
(78)

This expression is understood in terms of the zeta-function regularization and
$\tilde{g}^{(2)}_{k},\tilde{g}^{(2)}$ are different from
$g^{(2)}_{k},g^{(2)}$. Our goal is to reproduce ‘‘vacuum energy’’ ${\sf
Q}_{2n-1}^{0}$. There is infinitely many ways to do that, for example by
taking $g^{(2)}_{k}=0$, $\tilde{g}^{(2)}=g^{(2)}$, but we will additionally
require that the zeta-functions from (67) will become the sums of the form
$\sum_{k}k^{p}$ in (78). This leads to

$\displaystyle\tilde{g}^{(2)}_{k}=\sum_{j=0}^{n-1}{1\over
4}\xi_{n}^{j}\left((2n-1)y_{1}(j)-2j-1\right)\tilde{\Delta}^{n-1-j}\tilde{c}^{j}k^{2j+1},$
(79)

and very simple

$\displaystyle\tilde{g}^{(2)}=-\frac{n(n-1)(2n-1)\tilde{\Delta}^{n-1}}{96\tilde{c}}.$
(80)

This term is necessary to subtract $n_{k}$-independent
$\tilde{\Delta}^{n-1}\tilde{c}^{-1}$ term coming from
$\sum_{k}\tilde{g}^{(2)}_{k}\tilde{n}_{k}$ to match ${\sf Q}_{2n-1}^{0}$ (65)
which has no terms with the negative powers of $c$.

For convenience we give the full expression (78) explicitly

$\displaystyle{\sf Q}_{2n-1}$
$\displaystyle=\tilde{\Delta}^{n}+\sum_{k}\sum_{j=0}^{n-1}\frac{(2n-1)\sqrt{\pi}\,\Gamma(n+1)}{2\Gamma(j+\frac{3}{2})\Gamma(n-j)}\tilde{\Delta}^{n-1-j}\tilde{c}^{j}k^{2j+1}\tilde{n}_{k}$
(81) $\displaystyle\quad-\frac{n(n-1)(2n-1)\tilde{\Delta}^{n-1}}{96\tilde{c}}$
$\displaystyle\quad+\sum_{k}\sum_{j=0}^{n-1}\frac{(2n-1)\sqrt{\pi}\,\Gamma(n+1)}{8\Gamma(j+\frac{3}{2})\Gamma(n-j)}\left((2n-1)y_{1}(j)-2j-1\right)\tilde{\Delta}^{n-1-j}\tilde{c}^{j-1}k^{2j+1}\tilde{n}_{k}$
$\displaystyle\quad-\sum_{k}\sum_{j=0}^{n-1}\frac{(2n-1)\sqrt{\pi}\,\Gamma(n+1)(2nj+2n-3j-2)}{16\,\Gamma\left(j+\frac{3}{2}\right)\Gamma(n-j)}\tilde{\Delta}^{n-j-1}\tilde{c}^{j-1}k^{2j}\tilde{n}_{k}^{2}$
$\displaystyle\quad+\frac{1}{2}\sum_{k,\ell}\sum_{j=1}^{n-1}\frac{(2n-1)^{2}\sqrt{\pi}\,\Gamma(n+1)}{4\Gamma\left(j+\frac{3}{2}\right)\Gamma(n-j)}\tilde{\Delta}^{n-j-1}\tilde{c}^{j-1}\sum_{s=0}^{j-1}k^{2(j-s)-1}\ell^{2s+1}\tilde{n}_{k}\tilde{n}_{\ell}+\mathcal{O}(c^{n-3}).$

We conjecture this is the full quantum spectrum of $Q_{2n-1}$ up to
$\tilde{c}^{n-2}$ order and verify that it passes all available checks.

From here it is now straightforward to find ${\sf Q}_{2n-1}$ in the
representation (74). Coefficient

$\displaystyle
g^{(2)}_{k}=\sum_{j=0}^{n-1}\frac{(2n-1)\sqrt{\pi}\,\Gamma(n+1)}{8\Gamma\left(j+\frac{3}{2}\right)\Gamma(n-j)}v(n,j,k)\tilde{\Delta}^{n-j-1}\tilde{c}^{j-1},$
(82) $\displaystyle v(n,j,k)=$
$\displaystyle(2n-1)\sum_{s=0}^{j-1}\zeta(2(s-j)+1)k^{2s+1}+((2n-1)y_{1}(j)-2j-1)k^{2j+1}-{1\over
2}(2nj+2n-3j-2)k^{2j}$

is significantly more bulky than (79), while the full expression is

$\displaystyle{\sf Q}_{2n-1}$ $\displaystyle=$ $\displaystyle{\sf
Q}_{2n-1}^{0}+\sum_{k}\sum_{j=0}^{n-1}\xi_{n}^{j}\tilde{\Delta}^{n-j-1}\tilde{c}^{j}k^{2j+1}n_{k}$
$\displaystyle+$
$\displaystyle\sum_{k,\ell}\sum_{j=1}^{n-1}\xi_{n}^{j}{(2n-1)\over
4}\tilde{\Delta}^{n-j-1}\tilde{c}^{j-1}\sum_{s=0}^{j-1}k^{2(j-s)-1}\ell^{2s+1}n_{k}n_{\ell}$
$\displaystyle-$
$\displaystyle\sum_{k}\sum_{j=0}^{n-1}\xi_{n}^{j}{(2nj+2n-3j-2)\over
8}\tilde{\Delta}^{n-j-1}\tilde{c}^{j-1}k^{2j}n_{k}^{2}$ $\displaystyle+$
$\displaystyle\sum_{k}\sum_{j=0}^{n-1}\xi_{n}^{j}{v(n,j,k)\over
4}\tilde{\Delta}^{n-j-1}\tilde{c}^{j-1}n_{k}+\mathcal{O}(c^{n-3}).$

To summarize, we have found the (conjectured) spectrum of all qKdV charges at
first two sub-leading orders in $1/c$ expansion (81,4) and observed certain
patterns which may help fix the spectrum at higher orders. Let us spell the
step to find the next $1/\tilde{c}^{3}$ order, i.e. fix the terms of order
$\tilde{c}^{n-3}$ in (74). The classical result for $Q_{2n-1}$ in terms of
action variables $I_{k}$ was calculated up to cubic order in (49,52,53).
‘‘Energies’’ of primary states ${\sf Q}_{2n-1}^{0}$ were also calculated to
this order, see eq. (68). Thus $g^{(3)}_{k_{1}k_{2}k_{3}}$ and $g^{(3)}$ are
known, and to find the spectrum one would only need to fix
$g^{(3)}_{k_{1}k_{2}}$ and $g^{(3)}_{k}$. To do that one would need to find
$\tilde{g}^{(3)}_{k_{1}k_{2}}$ and $\tilde{g}^{(3)}_{k}$ from the expansion
(78) to reproduce (68) via zeta-function regularization and minimal possible
$\tilde{g}^{(3)}$, which presumably will only include terms with negative
powers of $\tilde{c}$. ‘‘Restoring’’ $\tilde{g}^{(3)}_{k_{1}k_{2}}$ and
$\tilde{g}^{(3)}_{k}$ from $\tilde{R}_{n,j}^{(3)}$ is not a mathematically
well-posed problem. We expect that all zeta-functions $\zeta(-2j-1)$ in
$\tilde{R}_{n,j}^{(3)}$ to lead to the sums $\sum_{k}k^{2j+1}\tilde{n}_{k}$ –
the rule which successfully worked at second $1/\tilde{c}$ order. At third
order this rule should be amended by others, as suggested by a non-polynomial
dependence on $k$ in (52). In practice, restoring $\tilde{g}^{(3)}_{k}$ from
$\tilde{R}_{n,j}^{(3)}$ may require establishing the analytic form of
coefficients $p_{j}$ in (68) and then reverse-engineering corresponding
$k_{1},k_{2},k_{3}$-dependent sums. Once hypothetical
$\tilde{g}^{(3)}_{k_{1}k_{2}}$ and $\tilde{g}^{(3)}_{k}$, and accordingly
${g}^{(3)}_{k_{1}k_{2}}$ and ${g}^{(3)}_{k}$ are fixed, a non-trivial set of
checks is provided by the spectrum of $Q_{3},Q_{5},Q_{7}$ generated by
computer algebra, as well as the requirement that thermal expectation values
$\langle Q_{2n-1}\rangle_{q}$ discussed in section 5.1 must have certain
modular properties.

###  4.1 Computer algebra check

For $n=1$ the expansion (4) reduces to (7) which is a simple check. A more
sophisticated check is provided by $Q_{3},Q_{5}$ and $Q_{7}$ which are known
explicitly in terms of the Virasoro algebra generators bazhanov1996integrable

$\displaystyle Q_{3}$ $\displaystyle=$ $\displaystyle\left(L_{0}^{2}-{c+2\over
12}L_{0}+{c(5c+22)\over 2990}\right)+\tilde{Q}_{3},$ (84)
$\displaystyle\tilde{Q}_{3}$ $\displaystyle=$ $\displaystyle
2\sum_{k=1}^{\infty}L_{-k}L_{k},$ $\displaystyle Q_{5}$ $\displaystyle=$
$\displaystyle\left(L_{0}^{3}-\frac{c+4}{8}L_{0}^{2}+\frac{(c+2)(3c+20)}{576}L_{0}-\frac{c(3c+14)(7c+68)}{290304}\right)+\tilde{Q}_{5},$
(85) $\displaystyle\tilde{Q}_{5}$ $\displaystyle=$
$\displaystyle\sum_{k,l=0}^{\infty}L_{-k-l}L_{k}L_{l}+2\sum_{k=1,l=0}^{\infty}L_{-k}L_{k-l}L_{l}+\sum_{k,l=1}^{\infty}L_{-k}L_{-l}L_{k+l}\,+$
(86) $\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{c+2}{6}m^{2}-\frac{c}{4}-1\right)L_{-n}L_{n}-L_{0}^{3},$

and dymarsky2019zero

$\displaystyle Q_{7}$ $\displaystyle=$
$\displaystyle\sum_{k,l,m=1}^{\infty}L_{-k}L_{-l}L_{-m}L_{k+l+m}+\sum_{k,l,m=0}^{\infty}L_{-k-l-m}L_{k}L_{l}L_{m}+$
(87) $\displaystyle 3\sum_{\begin{subarray}{c}k,l=1\\\
m=0\end{subarray}}^{\infty}L_{-k}L_{-l}L_{k+l-m}L_{m}+3\sum_{\begin{subarray}{c}k=1\\\
l,m=0\end{subarray}}^{\infty}L_{-k}L_{k-l-m}L_{l}L_{m}+$
$\displaystyle\frac{8+c}{3}\left[\sum_{k,l=1}^{\infty}(k+l)lL_{-k}L_{-l}L_{k+l}+\sum_{\begin{subarray}{c}k=1\\\
l=0\end{subarray}}^{\infty}(k-l)kL_{-k}L_{k-l}L_{l}\right]+$
$\displaystyle\frac{8+c}{3}\left[\sum_{k,l=0}^{\infty}(k+l)kL_{-k-l}L_{k}L_{l}+\sum_{\begin{subarray}{c}k=0\\\
l=1\end{subarray}}^{\infty}(k-l)kL_{-l}L_{l-k}L_{k}\right]+$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{c^{2}-c-141}{90}n^{4}-\frac{7c+59}{18}n^{2}\right)L_{-n}L_{n}-\left(\frac{1}{48}c^{2}+\frac{53}{360}c+\frac{19}{90}\right)\tilde{Q}_{3}-$
$\displaystyle\left(\frac{1}{6}c+1\right)\tilde{Q}_{5}-\frac{c+6}{6}L^{3}_{0}+\frac{15c^{2}+194c+568}{1440}L_{0}^{2}-$
$\displaystyle\frac{(c+2)(c+10)(3c+28)}{10368}L_{0}+\frac{c(3c+46)(25c^{2}+426c+1400)}{24883200}.$

Using computer algebra spectrum of $Q_{3},Q_{5},Q_{7}$ for all descendants at
a small levels $m$ can be evaluated explicitly, as an expansion in powers of
$1/c$. The resulting expressions can be compared with the spectrum following
from (4), which we will write in terms of quantum numbers $n_{k}$ packaged as
follows

$\displaystyle\qquad m_{p,r}\equiv\sum_{k}k^{p}n_{k}^{r},\qquad m_{p}\equiv
m_{p,1},\qquad m\equiv m_{1},\qquad h=\tilde{\Delta}/\tilde{c},$ (88)
$\displaystyle{\sf Q}_{3}$ $\displaystyle=$
$\displaystyle\tilde{\Delta}^{2}+\tilde{\Delta}\left(6\,m_{1}-\frac{1}{4}\right)+\tilde{c}\left(4\,m_{3}+\frac{1}{60}\right)+$
$\displaystyle\left(m_{3}-{3\over 2}m_{2}-{1\over
4}m_{1}\right)-\frac{3}{2}m_{2,2}+3m_{1}^{2}+\frac{3}{2}h(2m_{1}-m_{0}-m_{0,2})+{3\over
320}+{\mathcal{O}}(1/\tilde{c}),$

and similarly

$\displaystyle{\sf Q}_{5}$ $\displaystyle=$
$\displaystyle\tilde{\Delta}^{3}+\left(15m_{1}-\frac{5}{8}\right)\tilde{\Delta}^{2}+\tilde{\Delta}\,\tilde{c}\left(20m_{3}+\frac{1}{12}\right)+\tilde{c}^{2}\left(8m_{5}-\frac{1}{63}\right)+\,\,\,\quad\qquad\qquad\qquad$
$\displaystyle\tilde{\Delta}\left({5\over
12}(-5m_{1}-42m_{2}+44m_{3})-{35\over 2}m_{2,2}+25m_{1}^{2}+{15\over
2}h(2m_{1}-m_{0}-m_{0,2})+{23\over 192}\right)+$
$\displaystyle\tilde{c}\left({1\over
12}(m_{1}-10m_{3}-120m_{4}+64m_{5})-10m_{4,2}+20m_{1}m_{3}-{85\over
6048}\right)+\mathcal{O}(\tilde{c}^{0}),$

and

$\displaystyle{\sf Q}_{7}$ $\displaystyle=$
$\displaystyle\tilde{\Delta}^{4}+\tilde{\Delta}^{3}\left(28m_{1}-\frac{7}{6}\right)+\tilde{\Delta}^{2}\,\tilde{c}\left(56m_{3}+\frac{7}{30}\right)+\tilde{\Delta}\,\tilde{c}^{2}\left(\frac{224}{5}m_{5}-\frac{4}{45}\right)+\tilde{c}^{3}\left(\frac{64}{5}m_{7}+\frac{2}{75}\right)+$
(91) $\displaystyle\tilde{\Delta}^{2}\left({7\over
6}\left(-7m_{1}-66m_{2}+76m_{3}\right)-77m_{2,2}+98m_{1}^{2}+21h(2m_{1}-m_{0}-m_{0,2})+{259\over
480}\right)+$ $\displaystyle\tilde{\Delta}\,\tilde{c}\left({7\over
75}\left(7m_{1}-70m_{3}-960m_{4}+688m_{5}\right)-{448\over 5}m_{4,2}+{784\over
5}m_{1}m_{3}-\frac{167}{1080}\right)+$
$\displaystyle\tilde{c}^{2}\left({2\over
225}(-10m_{1}+21m_{3}-210m_{5}-3780m_{6}+1704m_{7})-{168\over
5}m_{6,2}+{112\over 5}(2m_{1}m_{5}+m_{3}^{2})+\frac{77}{2160}\right)$
$\displaystyle+$ $\displaystyle\mathcal{O}(\tilde{c}^{1}).$

We checked, these expressions are in agreement with the computer algebra
generated spectrum for $m\leq 12$, which serves as a non-trivial consistency
check of (4).

##  5 Miscellaneous results

Explicit expression for the spectrum of quantum $Q_{2n-1}$ in large $c$ limit
opens the opportunity to make progress in a number of adjacent directions. In
this section we discuss several applications of our results.

###  5.1 Thermal expectation values of $Q_{2n-1}$

Our first application is toward thermal exaction value of $Q_{2n-1}$, i.e.
averaged over the CFT Gibbs ensemble $\langle Q_{2n-1}\rangle_{q}\equiv{\rm
Tr}(q^{L_{0}-c/24}Q_{2n-1})$. This question appears naturally, though in a
more complicated form, to calculate the averaged value of $Q_{2n-1}$ over the
KdV Generalized Gibbs Ensemble (see section 5.2 below), if one wants to match
the GGE chemical potentials to describe equilibration endpoint of some initial
state. The expectation value $\langle Q_{2n-1}\rangle_{q}$, which is
essentially the one-point function of $T_{2n}$ (1) on the torus, exhibits
modular properties and can be represented as a covariant differential operator
acting on the CFT torus partition function maloney2018thermal . In fact, one
can average $Q_{2n-1}$ over a particular Verma module, $\langle
Q_{2n-1}\rangle_{\Delta}\equiv{\rm Tr}_{\Delta}(q^{L_{0}-c/24}Q_{2n-1})$,
where sum goes over all Virasoro descendants of the primary state
$|\Delta\rangle$. This sum too is a modular object and can be evaluated with
help of the same differential operator

$\displaystyle\langle Q_{2n-1}\rangle_{\Delta}$ $\displaystyle=$
$\displaystyle{\mathcal{D}}_{n}\chi_{\Delta},\qquad\chi_{\Delta}\equiv{\rm
Tr}_{\Delta}(q^{L_{0}-c/24})={q^{\tilde{\Delta}-\frac{1}{24}}/\eta},$ (92)
$\displaystyle{\mathcal{D}}_{n}$ $\displaystyle=$ $\displaystyle
D^{n}+\sum_{j=1}^{n-1}P_{n}^{j}(c,q)D^{n-j-1},\qquad D^{n}=D_{2(n-1)}\dots
D_{2}D_{0},$ (93)

and $D_{r}=q\partial_{q}-{r\over 12}E_{2}$ is Serre derivative. Each
$P_{n}^{j}$ is a degree $j$ polynomial in $c$ with each coefficient being a
modular form of weight $2j+2$,

$\displaystyle
P_{n}^{j}(c,q)=\sum_{k=1}^{j+1}P^{(k)}_{n,j}\tilde{c}^{j-k+1}E_{2j+2}^{(n,k)}(q).$
(94)

Here $P_{n,j}^{(k)}$ are numerical coefficients and $E_{2j+2}^{(n,k)}$ is some
modular form, which is a linear combination of $E_{4}^{a}E_{6}^{b}$ with
$4a+6b=2j+2$ for non-negative integer $a,b$, normalized such that
$E_{2j+2}^{(n,j)}=1+O(q)$. For $j=1,2,3,4,6$ there is a unique modular form of
the weight $2(j+1)$ and therefore for these $j$, independently of $n$ and $k$,
$E_{2j+2}^{(n,k)}=E_{2j+2}$ where

$\displaystyle
E_{2n}=1+{2\over\zeta(1-2n)}\sigma_{2n-1},\qquad\sigma_{p}=\sum_{k=1}^{\infty}{k^{p}q^{k}\over
1-q^{k}}.$ (95)

For instance, in the simplest case of $Q_{3}$ the operator $\mathcal{D}_{2}$
is given by

$\displaystyle\langle Q_{3}\rangle_{\Delta}$ $\displaystyle=$
$\displaystyle{\mathcal{D}}_{2}\chi_{\Delta}=\left[D^{2}+\frac{c}{1440}E_{4}\right]\chi_{\Delta}.$
(96)

In this case $P_{2,1}^{(1)}=1/60$ and $P_{2,1}^{(2)}=1/1440$. Explicit
expressions for $\mathcal{D}_{n}$ for $n\leq 7$ were found in
maloney2018thermal . For higher $n$ the modular form $E_{2i+2}^{(n,j)}$ and
coefficients $P_{n,j}^{(k)}$ are not known.

Strictly speaking (92,93) is an unproven ansatz proposed in maloney2018thermal
. We find it to be consistent with the large $c$ spectrum of $Q_{2n-1}$ (4)
and fix two leading in $c$ terms in $P_{n}^{j}$. To compare with (92), we need
to calculate $\langle Q_{2n-1}\rangle_{\Delta}$ starting from (4). Here the
following straightforward identities will be helpful

$\displaystyle\langle\sum_{k=1}^{\infty}n_{k}k^{p}\rangle_{\Delta}=\sigma_{p}\chi_{\Delta},\quad\langle\sum_{k=1}^{\infty}n^{2}_{k}k^{p}\rangle_{\Delta}=\left(2q\partial_{q}\sigma_{p-1}-\sigma_{p}\right)\chi_{\Delta},\quad$
(97)
$\displaystyle\langle\sum_{k=1}^{\infty}n_{k}k^{p}\sum_{\ell=1}^{\infty}n_{\ell}\ell^{p^{\prime}}\rangle_{\Delta}=\left(q\partial_{q}\sigma_{p+p^{\prime}-1}+\sigma_{p}\sigma_{p^{\prime}}\right)\chi_{\Delta},$
(98)

where by $n_{k}$ we mean the quantum numbers (4). Then (1) immediately yields

$\displaystyle\langle
Q_{2n-1}\rangle_{\Delta}=\tilde{\Delta}^{n}\chi_{\Delta}+\sum_{j=0}^{n-1}\tilde{\Delta}^{n-p-1}\tilde{c}^{p}\xi_{n}^{p}\left(\sigma_{2p+1}+{\zeta(-2p-1)\over
2}\right)\chi_{\Delta}+O(\tilde{c}^{n-2}),$ (99)

where we assumed the usual limit, $h=\tilde{\Delta}/\tilde{c}$ is kept fixed
while $\tilde{c}\rightarrow\infty$. Comparing this with (92), we immediately
see that the leading $\tilde{\Delta}^{n}$ term is coming from (we drop
$\chi_{\Delta}$ for simplicity)

$\displaystyle
D^{n}\rightarrow(q\partial_{q})^{n}\rightarrow\tilde{\Delta}^{n}.$ (100)

Similarly we can trace origin of all $\tilde{c}^{n-1}$ terms,

$\displaystyle D^{n}\rightarrow(q\partial_{q})^{n}-{n(n-1)\over
12}E_{2}(q\partial_{q})^{n-1}\rightarrow\tilde{\Delta}^{n-1}n\left(\sigma_{1}-{1\over
24}\right)-\tilde{\Delta}^{n-1}{n(n-1)\over 12}E_{2}=-{n(2n-1)\over 24}E_{2},$

which agrees with (99), and

$\displaystyle P^{(1)}_{n,j}\tilde{c}^{j}E^{(n,1)}_{2j+2}D^{n-j-1}\rightarrow
P^{(1)}_{n,j}\tilde{c}^{j}E^{(n,1)}_{2j+2}(q\partial_{q})^{n-j-1}\rightarrow
P^{(1)}_{n,j}\tilde{\Delta}^{n-j-1}\tilde{c}^{j}E^{(n,1)}_{2j+2},$ (101)

for $n-1\geq j>0$. From here immediately follows

$\displaystyle P^{(1)}_{n,j}=\tilde{R}^{(1)}_{n,j},\qquad
E_{2j+2}^{(n,1)}=E_{2j+2},\qquad n-1\geq j\geq 1.$ (102)

To fix $P^{(2)}_{n,j}$ it is convenient to take $q\rightarrow 0$ limit and
compare $\langle Q_{2n-1}\rangle_{\Delta}$ with (65), yielding

$\displaystyle
P^{(2)}_{n,1}=\tilde{R}^{(2)}_{n,0}-\frac{n(n-1)(12n^{2}-16n-1)}{3456}={n(n-1)(12n^{2}-38n+31)\over
8640},$ $\displaystyle
P^{(2)}_{n,j}=\tilde{R}^{(2)}_{n,j-1}+{(n-j)(2(n-j)-1)\over
24}P^{(1)}_{n,j-1},\quad n-1\geq j\geq 2.$

Evaluation of $E^{(n,2)}_{2j+2}$ is a more challanging task and requires first
using (97,98) and then combining pieces into modular forms to match (92,93).
We note, there are terms in (4) proportional to
$\tilde{\Delta}^{n-1}\tilde{c}^{-1}$, but (92) has no negative powers of $c$.
Hence these terms must vanish after averaging, which follows from the identity
$q\partial_{q}\sigma_{-1}-\sigma_{1}=0$ and serves as a consistency check. The
final expression reads

$\displaystyle P^{(2)}_{n,j}E^{(n,2)}_{2j+2}={(2n-1)\sqrt{\pi}\Gamma(n+1)\over
8\Gamma(j+3/2)\Gamma(n-j)}\left(((2n-1)y_{1}(j)-2j-1){\zeta(-2j-1)\over
2}E_{2j+2}\,\,-\right.$ (103)
$\displaystyle(n-1-j)\zeta(-2j+1)D_{2j}E_{2j}+\left.{(2n-1)\over
4}\sum_{s=1}^{j-2}\zeta(-2s-1)\zeta(-2(j-s)+1)E_{2s+2}E_{2(j-s)}\right).$

It is valid for $n-1\geq j\geq 2$. For $j=1$, there is a unique modular form
$E^{(n,2)}_{2j+2}=E_{4}$. Also, as was mentioned above
$E^{(n,2)}_{2j+2}=E_{2j+2}$ for $j=2,3,4,6$, which can be checked
straightforwardly. Because of the identities between modular forms there are
other ways to write (103).

Explicit form of ${\sf Q}_{2n-1}^{0}$ up to $\tilde{c}^{n-3}$ order allows us,
in principle, to calculate $P^{(3)}_{n,j}$, although calculation of
$E^{(n,3)}_{2j+2}$ would require first extending (4) to the next $1/c$ order.
Given involved form of $P^{(2)}_{n,j}$ and $E^{(2)}_{n,j}$ we do not expect
the answer to be simple.

###  5.2 Generalized Gibbs Ensemble

Spectrum of $Q_{2n-1}$ can help understand the qKdV generalized Gibbs ensemble
(GGE)

$\displaystyle\rho_{\rm GGE}=e^{-\sum_{n}\mu_{2n-1}\,Q_{2n-1}},\qquad Z_{\rm
GGE}=\operatorname{Tr}\rho_{\rm GGE},$ (104)

and corresponding (generalized) partition function and free energy. Earlier
attempts to evaluate KdV generalized free energy include de2016remarks ;
maloney2018generalized ; GGE ; GGE2 . The GGE describes local equilibrium in a
state carrying specific values of qKdV charges. It is expected on general
grounds that most initial states, upon equilibration, can be locally described
by the GEE with the appropriate values of chemical potentials $\mu_{2n-1}$
rigol2007relaxation . From the mathematical point of view, it is of great
interest to investigate modular properties of $Z_{\rm GGE}$, generalizing
modular invariance of the conventional torus partition function
$\mu_{2n-1}=0$, for $n>1$.

The explicit spectrum of $Q_{2n-1}$ in the large $c$ limit allows in principle
to calculate the generalized sum over a particular Verma module

$\displaystyle\operatorname{Tr}_{\Delta}e^{-\sum_{n}\mu_{2n-1}Q_{2n-1}}$ (105)

in the ‘‘holographic limit’’: $h=\tilde{\Delta}/\tilde{c}$,
$t_{n}:=\mu_{2n-1}/\tilde{c}^{n-1}$ fixed, $\tilde{c}\rightarrow\infty$, by
expanding the answer in powers of $1/c$. In practice sums of exponents of
quadratic or higher order expressions in $n_{k}$

$\displaystyle\sum_{n_{k}}e^{O(n^{2})}$ (106)

can not be evaluated, and we restrict our analysis to first non-trivial $1/c$
order,

$\displaystyle\operatorname{Tr}_{\Delta}e^{-\sum_{n}\mu_{2n-1}Q_{2n-1}}=e^{-\tilde{c}\sum_{n}t_{n}h^{n}}{e^{-\sum_{n}t_{n}\sum_{p=0}^{n-1}h^{n-1-p}\xi_{n}^{p}\zeta(-2p-1)/2}\over\prod_{k=1}^{\infty}\left(1-e^{-k\sum_{n}(2n-1)n\,t_{n}\,h^{n-1}{}_{2}F_{1}(1,1-n,3/2,-k^{2}/h)}\right)}.$

From here generalized partition function can be evaluated using Cardy formula
(we are only writing explicitly the chiral part),

$\displaystyle Z_{\rm GGE}$ $\displaystyle=$ $\displaystyle
e^{\tilde{c}f_{0}+f_{1}+O(1/\tilde{c})},$ (107) $\displaystyle f_{0}$
$\displaystyle=$ $\displaystyle\sum_{n=1}^{\infty}(2n-1)\,t_{n}\,h^{n},$ (108)
$\displaystyle h^{1/2}$ $\displaystyle=$ $\displaystyle{1\over
2\pi}\sum_{n=1}^{\infty}t_{n}\,n\,h^{n},$ (109) $\displaystyle f_{1}$
$\displaystyle=$
$\displaystyle-\sum_{k=1}^{\infty}\ln\left(1-e^{-\gamma}\right)-\sum_{n=2}^{\infty}t_{n}\,h^{n-1}\left(\sum_{p=0}^{n-1}\xi_{n}^{p}\,h^{-p}{\zeta(-2p-1)\over
2}-{n\over 24}\right),$ (110) $\displaystyle\gamma(k)$ $\displaystyle=$
$\displaystyle
k\sum_{n=1}^{\infty}(2n-1)n\,t_{n}\,h^{n-1}{}_{2}F_{1}(1,1-n,3/2,-k^{2}/h).$
(111)

Here $Z_{\rm GGE}$ is understood to be a function of
$t_{n}\equiv\mu_{2n-1}\tilde{c}^{1-n}$, while $h$ is a function of $t_{n}$
satisfying (109). For (107)-(111) to be valid, resulting
$\tilde{\Delta}=\tilde{c}\,h$ should be in the regime of validity of Cardy
formula. There are at least two limits when this assumption is controllable.
First, (107) is valid for any large $c$ theory in the thermodynamic limit. We
introduce the spatial circle radius $L$ (we kept $L=1$ in the paper so far)
and inverse temperature $\beta$, $\mu_{1}=t_{1}=\beta/L$. By taking
$L\rightarrow\infty$, while all other chemical potentials scale as
$\mu_{2n-1}\propto t_{n}\sim L^{1-2n}$ to ensure that values of all
$Q_{2n-1}\propto L$ are extensive, we find the saddle point value $h\sim
L^{2}$ and $f_{0},f_{1}\sim\L$. (The scaling of $f_{1}$ follows by
substituting the sum over $k$ in (110) by an integral over $\kappa=k^{2}/h$.)
In this limit second term in (110), the sum over $n$, is sub-extensive and can
be neglected. We therefore arrive at the leading (extensive) contribution to
$f_{0}$ and $f_{1}$ found in GGE2 .

Second case when (107)-(111) can be trusted is in holographic theories, i.e.
large $c$ theories satisfying HKS sparseness condition hartman2014universal .
From the holographic point of view $f_{0}$ is the free energy of BTZ black
hole in the Euclidean classical theory of gravity with the deformed boundary
conditions such that the dual CFT Hamiltonian is
$H=\sum_{n}\mu_{2n-1}Q_{2n-1}$ perez2016boundary ; GETH ; Dymarsky_2020 . The
leading correction $f_{1}$ can be interpreted as the one-loop contribution
coming from the boundary gravitons. Different solutions of (109) means
Euclidean path integral could have numerous BTZ saddles and the condition
$h(t_{n})>1/12$ necessary for the validity of Cardy formula would come
automatically as the requirement of smoothness of bulk geometry.

It is possible to fine-tune chemical potentials $t_{n}$ such that $\gamma(k)$
(111) for some $k$ will vanish. That will render $f_{1}$ divergent, indicating
higher order $1/c$ corrects are necessary to make free energy finite.
Schematically, the spectrum ${\sf Q}_{2n-1}$ is an expansion in
$n_{k}/\tilde{c}$. For the higher order corrections to contribute at the
leading order, the quantum numbers $n_{k}$ should be of order $\tilde{c}$. In
terms of the classical problem of section (2), action variables $I_{k}$ should
be of order one rather than infinitesimal. In other words, leading
contribution would come from a non-trivial saddle when classical $u(\varphi)$
is not a constant but some solitonic solution. Such saddles, describing black
holes, which are geometrically different from the BTZ configurations, were
constructed in Dymarsky_2020  and it was shown that for certain parameters
$\mu_{2n-1}$ they give leading contribution to generalized free energy. We
dubbed these configurations ‘‘KdV-charged’’ black holes to emphasize that
higher KdV charges $Q_{2n-1}$, even at leading order in $c$, are different
from $Q_{1}^{n}$, unlike for BTZ configurations for which $u(\varphi)=u_{0}$
is a constant and $Q_{2n-1}\sim u_{0}^{n}$.

Theoretical control over generalized free energy in the large $c$ limit can be
used to probe modular properties of $Z_{\rm GGE}$. The currents $T_{2n}$ (1)
have no anomalous dimension and therefore naively $Z_{\rm GGE}$ should be
invariant under modular transformation $t_{1}\rightarrow
t_{1}^{\prime}=(2\pi)^{2}/t_{1}$ accompanied by

$\displaystyle t_{n}\rightarrow(-1)^{n}\left({2\pi\over
t_{1}}\right)^{2n}t_{n},\quad n>1.$ (112)

This only holds to linear order in $t_{n},n>1$, i.e. at the level of thermal
expectation values $\langle Q_{2n-1}\rangle_{q}$ discussed in section 5.1. At
higher orders invariance is broken due to colliding $T_{2n}$
maloney2018thermal . To restore invariance of $Z_{\rm GGE}$, while working in
the $c\rightarrow\infty$ limit one may require $f_{0}$ given by (108),(109) to
be invariant under the hypothetical transformation $t_{n}\rightarrow
t_{n}^{\prime}(t_{n},t_{1})$, $n>1$. More accurately, in addition to BTZ black
holes described by (108),(109) we should include vacuum (thermal AdS3) and
KdV-charged black holes to the list of possible saddles. Given a non-trivial
diagram of the Hawking-Page phase transitions, to match leading saddles, the
hypothetical transformation $t_{n}\rightarrow t_{n}^{\prime}(t_{n},t_{1})$
should be very complicated, with numerous branches of continuity. This may
indicate that in the presence of higher KdV charges modular invariance of
$Z_{\rm GGE}$ is not mathematically natural. Similar conclusion is recently
reached in Downing:2021mfw , which evaluated $Z_{\rm GGE}$ explicitly in the
case of $c=1/2$ free fermion model. They found that to reproduce $Z_{\rm GGE}$
in the dual channel, one needs to sum over not one but three fermion Hilbert
spaces, schematically $Z_{\rm GGE}(t_{1},t_{3})\propto
Z_{1}(t_{1}^{\prime})Z_{2}(t_{1}^{\prime})Z_{3}(t_{1}^{\prime})$, a
mathematical observation (conjecture), which so far has no physical
interpretation. To summarize, failure to establish invariance of $Z_{\rm
GGE}(t)$ under modular transformation supplemented by an appropriate map
$t_{n}\rightarrow t_{n}^{\prime}$ in both infinite $c$ limit and for $c=1/2$
model may suggest that it is not mathematically natural and instead covariance
of $Z_{\rm GGE}(t)$ under (112) should be investigated.

###  5.3 Transfer Matrix

In the classical case, as follows from (34), charges $Q_{2n-1}$ encode
asymptotic expansion of the quasi-momentum $p(\lambda)$. The quasi-momentum
controls the eigenvalues $e^{\pm 2\pi ip(\lambda)}$ of the monodromy matrix of
the differential equation (14). Instead of $p(\lambda)$ one can consider the
trace of monodromy matrix

$\displaystyle T(\lambda)=2\cos(2\pi p(\lambda)).$ (113)

In case of the constant potential $u(\varphi)=h$ this becomes
$T(\lambda)=2\cos(2\pi\sqrt{\lambda-h/4})$.

In quantum case $T(\lambda)$ becomes the transfer matrix, which is related to
qKdV charges via an asymptotic expansion bazhanov1996integrable

$\displaystyle\ln T$ $\displaystyle=$ $\displaystyle\kappa\,\mu^{1/2}\,{\rm
1}-\sum_{n}C_{n}\,\mu^{1/2-n}\,Q_{2n-1},\quad\mu\rightarrow\infty,$ (114)

where

$\displaystyle\kappa$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{\pi}\Gamma\left(\frac{1}{2}-\frac{\xi}{2}\right)}{\Gamma\left(1-\frac{\xi}{2}\right)}\left(\Gamma\left(1-\beta^{2}\right)\right)^{1+\xi},$
(115) $\displaystyle C_{n}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\pi}(1+\xi)\beta^{2n}\Gamma\left((n-\frac{1}{2})(1+\xi)\right)}{\Gamma(n+1)\Gamma\left(1+(n-\frac{1}{2})\xi\right)}\left(\Gamma\left(1-\beta^{2}\right)\right)^{-(2n-1)(1+\xi)},$
(116) $\displaystyle\beta$ $\displaystyle\equiv$
$\displaystyle\sqrt{\frac{1-c}{24}}-\sqrt{\frac{25-c}{24}},\qquad\xi\equiv\frac{\beta^{2}}{1-\beta^{2}}.$
(117)

Variable $\mu$ will become spectral parameter $-\lambda$ in the classical
limit. The original paper bazhanov1996integrable  introduces another variable
$\lambda$ , defined as $\mu\equiv\lambda^{2(1+\xi)}$. We use this definition
in the reminder of this section.

We are interested in the limit $c\rightarrow\infty$, or $\beta\rightarrow 0$.
Following bazhanov1996integrable  we introduce $p^{2}=\beta^{2}\tilde{\Delta}$
which remains finite in this limit finite,
$p^{2}\rightarrow-\tilde{\Delta}/4\tilde{c}=-h/4$. (This is, obviously, a
different quantity from the quasi-momentum $p(\lambda)$ mentioned above.) We
would like to find $T$ by summing the asymptotic expansion (114) while
expanding it in powers of $\beta^{2}$ which corresponds to $1/c$ expansion. In
principle we can use the spectrum (4) to calculate $\ln T$ acting on an
excited state, but resort to a simpler calculation for $\ln T$ acting on a
primary state. In this case $Q_{2n-1}$ in (114) should be substituted by ${\sf
Q}_{2n-1}^{0}$, which we expand in powers of $\beta^{2}\propto 1/c$ (65). The
calculation is tedious and we only give the final expression

$\displaystyle(\ln T)_{\rm asympt}|\Delta\rangle$ $\displaystyle=$
$\displaystyle 2\pi i\sqrt{p^{2}-\lambda^{2}}\Phi(\lambda,p)|\Delta\rangle,$
(118) $\displaystyle\Phi(\lambda,p)$ $\displaystyle=$ $\displaystyle
1+\beta^{2}\Psi-\frac{\beta^{4}}{48\lambda^{2}(p^{2}-\lambda^{2})^{3}}\left[3\lambda^{4}p^{2}+2\pi^{2}\lambda^{4}(p^{2}-\lambda^{2})(4p^{2}-3\lambda^{2})\right]$
(119)
$\displaystyle\quad+\frac{\beta^{4}}{2\lambda^{2}}\left[(2p^{2}+\lambda^{2})\Psi^{2}-2(p^{2}-\lambda^{2})\lambda\Psi\frac{d\Psi}{d\lambda}\right]+\mathcal{O}(\beta^{6}),$
$\displaystyle\Psi(\lambda,p)$ $\displaystyle=$
$\displaystyle\frac{\lambda^{2}}{\lambda^{2}-p^{2}}\left[\gamma+\frac{1}{2}\psi\left(2\sqrt{p^{2}-\lambda^{2}}\right)+\frac{1}{2}\psi\left(-2\sqrt{p^{2}-\lambda^{2}}\right)\right].$
(120)

Here $\psi$ is the polylog function.

Given analytic form of (118) it is tempting to extend its validity from the
asymptotic regime $\lambda\rightarrow\infty$ to the vicinity of $\lambda=0$.
This is clearly wrong as even in the strict classical limit $\beta\rightarrow
0$ we do not recover correct classical expression for the trace of monodromy
matrix simply from $e^{(\ln T)_{\rm asympt}}$. Yet in the limit
$\beta\rightarrow 0$ the correct answer is reproduced by the following simple
conjectural expression (we implicitly assume this is an eigenvalue of $T$
acting on $|\Delta\rangle$),

$\displaystyle T_{\rm guess}(\lambda,\beta,p)=e^{(\ln T)_{\rm
asympt}}+e^{-(\ln T)_{\rm asympt}},$ (121)

and we would like to check if it could be valid beyond the strict $\beta=0$
limit. To that end we expand (121) in powers of $\beta$ (amended by an
expansion in $\lambda$), to find

$\displaystyle T_{\rm
guess}=2\cos(2\pi\sqrt{p^{2}-\lambda^{2}})+\beta^{2}\left[\frac{2\pi\sin(2\pi
p)}{p}\left(2\gamma+\psi(2p)+\psi(-2p)\right)\lambda^{2}+\right.$
$\displaystyle\left.\left(-\frac{2\pi^{2}\cos(2\pi
p)}{p^{2}}\left[2\gamma+\psi(2p)+\psi(-2p)\right)\right.\right.$
$\displaystyle\left.\left.+\frac{\pi\sin(2\pi
p)}{p^{3}}\left(2\gamma+\psi(2p)+\psi(-2p)-2p\psi^{(1)}(2p)+2p\psi^{(1)}(-2p)\right)\right)\lambda^{4}+O(\lambda^{6})\right]+$
$\displaystyle\beta^{4}\left[\frac{\pi\sin(2\pi
p)}{12p^{3}}\left[3+8\pi^{2}p^{2}+12p^{2}\left(2\gamma+\psi(2p)+\psi(-2p)\right)^{2}\right]\lambda^{2}+O(\lambda^{4})\right]+O(\beta^{6}).$

Small $\lambda$ expansion of the actual $T$ is given in bazhanov1996integrable
in the explicit form in terms of the integrals of free field correlators. A
comparison with $T_{\rm guess}$ reveals that, besides the classical
$\beta^{0}$ term, which matches the classical expression (113) for a constant
potential $u=h$, only $\beta^{2}\lambda^{2}$ term coincides with, while
$\beta^{2}\lambda^{4}$ and $\beta^{4}\lambda^{2}$ terms do not match the
correct result. We thus conclude that the conjectural expression (121) is
missing non-perturbative terms, which are not captured by the asymptotic
expansion (114).

##  6 Discussion

In this paper we obtained spectrum of quantum KdV charges $Q_{2n-1}$ in first
two non-trivial orders in $1/c$ expansion. Our result (81) and (4) is valid in
the semiclassical limit of large central charge $c\rightarrow\infty$ with the
ratio of $\Delta/c$ kept fixed. This limit is inspired by holographic
correspondence, when CFT is dual to weakly coupled gravity. Accordingly,
dynamics of stress-energy sector becomes semiclassical, with the leading
(classical) contribution governed by integrable dynamics on the co-adjoint
orbit of the Virasoro algebra. Under semiclassical quantization classical
action variables $I_{k}$ are promoted to integer quantum numbers $n_{k}$, and
the spectrum of $Q_{2n-1}$ looks most elegant in terms of variables
$\tilde{\Delta}$ and $\tilde{c}$ (70). At each order in $1/\tilde{c}$ the
quantum answer is a polynomial in $n_{k}$. Classical calculation fixes the
leading term with the highest power of $n_{k}$, while all other terms should
be regarded as ‘‘quantum corrections.’’ We have seen that semiclassical
quantization, combined with the values of qKdV charges $Q_{2n-1}$ acting on
primary states, is sufficient to completely fix these quantum corrections and
obtain the spectrum of excited states at least in first two orders in $1/c$.
We conjecture this quantization scheme can be extended to higher orders in
$1/c$. We laid the groundwork for the next order $1/c^{3}$ by calculating
classical $Q_{2n-1}(h,I_{k})$ as well as ‘‘energies’’ on primary states ${\sf
Q}^{0}_{2n-1}$, albeit in the latter case not all terms are known
analytically. To complete the job one would need to find analytic expressions
for ${\sf Q}^{0}_{2n-1}$ and develop a dictionary that maps each term to an
infinite sum, yielding this term back via zeta-function regularization.

It is tempting to interpret quantization of $Q_{2n-1}$ holographically, as a
semiclassical quantization of boundary gravitons in AdS3. We develop this
picture at first $1/c$ order in the appendix A, but holographic picture does
not provide any immediate insight into ‘‘quantum corrections’’ appearing at
higher orders in $1/c$.

The obtained spectrum has several immediate applications. First, in section
5.1 we calculated two leading terms in large $c$ expansion of the ‘‘thermal
expectation values’’ $\langle Q_{2n-1}\rangle_{\Delta}\equiv{\rm
Tr}_{\Delta}(q^{L_{0}-c/24}Q_{2n-1})$, where sum goes over a particular Verma
module, and compared them with the predictions of maloney2018thermal .
Covariance under modular transformation of $\langle Q_{2n-1}\rangle_{\Delta}$
in each order in $1/c$ serves as a non-trivial check of our main result (4).
We also fixed two leading terms in the differential operator
${\mathcal{D}}_{n}$ yielding thermal expectation values via $\langle
Q_{2n-1}\rangle_{\Delta}={\mathcal{D}}_{n}{\rm Tr}_{\Delta}(q^{L_{0}-c/24})$,
see (102) and (103). Second, in section 5.2 we calculated first $1/c$
correction to generalized free energy of the qKdV Generalized Gibbs Ensemble

$\displaystyle Z_{\rm GGE}=\operatorname{Tr}e^{-\sum_{n}\mu_{2n-1}Q_{2n-1}}.$
(122)

The latter describes local equilibrium of a 2d CFT in a state carrying
specific values of qKdV charges. It is of great interest to further
investigate mathematical properties of $Z_{\rm GGE}$, in particular covariance
under modular transformation. Third, in section 5.3 using asymptotic expansion
we calculated quantum transfer matrix acting on a primary state in first two
non-trivial orders in $1/c$ expansion. Unfortunately the obtained expression
is lacking terms non-perturbative in spectral parameter, which can not be
fixed from the knowledge of spectrum of $Q_{2n-1}$ alone.

There are several potential applications of our results, which we hope to
address in the future. The obtained spectrum of $Q_{2n-1}$ will be helpful to
study generalized Eigenstate Thermalization Hypothesis of 2d CFTs GETH  at the
subleading order in $1/c$. We also expect the semiclassical quantization
approach developed in this paper could be helpful in the context of
Intermediate Long Wave hierachry, which is closely related to qKdV problem.
More generally, it would be interesting to bridge the gap between the
semiclassical approach of this work with the Bethe anzatz approach of
LitvinovEq  by taking ‘‘holographic’’ limit $c\rightarrow\infty$ with fixed
$h=\tilde{\Delta}/\tilde{c}$ of the appropriate Bethe anzatz equations.

###### Acknowledgements.

We thank A. Gorsky, I. Krichiver, A. Litvinov, N. Nekrasov, A. Okounkov, V.
Pestun and A. Zamolodchikov for helpful discussions. AD, AK, and SS were
supported by the National Science Foundation under Grants No. PHY-1720374 and
PHY-2013812. AD is grateful to IHES and KITP for hospitality, where this work
was partially done. The research at KITP was supported in part by the National
Science Foundation under Grant No. PHY-1748958. The research at IHES was
funded from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (QUASIFT grant agreement 677368).
KP was supported by the RFBR grant 19-32-90173.

##  Appendix A Spectrum of linear perturbations from AdS3

In the classical (infinite central charge) limit gravity in AdS3 can be
described in terms of two functions $u(t,\varphi)$ and $\bar{u}(t,\varphi)$
living at the boundary and satisfying EOM, $\dot{u}=\partial_{\varphi}u$ and
$\dot{\bar{u}}=-\partial_{\varphi}\bar{u}$. That is in the conventional case,
when the dual CFT’s Hamitlonian is
$H=Q_{1}+\bar{Q}_{1}=L_{0}+\bar{L}_{0}-{c/12}$. Should the Hamiltonian be
chosen to be one of the higher qKdV charges, $H=Q_{2n-1}+\bar{Q}_{2n-1}$,
functions $u,\bar{u}$ will be satisfying higher KdV equations

$\displaystyle\dot{u}={c\over 24}\\{Q_{2n-1},\delta u\\}=(2n-1)\partial
R_{n},$ (123)

and similarly for $\bar{u}$ perez2016boundary ; Ojeda:2019xih ; Dymarsky_2020
. In this case the spectrum of $Q_{2n-1}$ is the spectum of small fluctuations
of $u$ above the constant backgroun $u=u_{0}$, that corresponds to unpertubed
metric in AdS3. In other words, to quantize $Q_{2n-1}$ we consider linearized
EOM for small fluctuations $u=u_{0}+\delta u$, where $\delta u\propto
e^{i\varepsilon t+ik\varphi}$ is a flat wave. We want to find energy
$\varepsilon$ of the flat wave which satisfies the equation of motion (123)

$i\varepsilon_{n}\delta u=(2n-1)\delta\partial R_{n}.$ (124)

For example in the case $n=1$ we have $\varepsilon_{1}=k$, in case $n=2$ we
have $\varepsilon_{2}=2u_{0}n+{4\over 3}k^{3}$ and so on. In general we can
get from (32)

$\displaystyle\delta\partial R_{n+1}={n+1\over 2n+1}\delta(\partial
u+2u\partial-2\partial^{3})R_{n}=i{n+1\over
2n+1}(kR_{n}(u_{0})+2(u_{0}+k^{2})\delta\partial R_{n}).\quad$ (125)

Hence we find the following iterative relation for $\varepsilon_{n}$

$\varepsilon_{n+1}=(2n-1){n+1\over
2n+1}\left[2(u_{0}+k^{2})\varepsilon_{n}+ku_{0}^{n}\right],$ (126)

where we have used that $R_{n}(u_{0})=u_{0}^{n}$. Each $\varepsilon_{n}$ is a
polynomial of the form

$\varepsilon_{n}=\sum_{p=0}^{n-1}\xi^{p}_{n}\,k^{2p+1}u_{0}^{n-1-p},$ (127)

where $\zeta_{n}^{p}$ satisfy

$\xi_{n+1}^{p}=(2n-1){2(n+1)\over
2n+1}\left(\xi_{n}^{p}+\xi_{n}^{p-1}\right),$ (128)

and we defined $\xi_{n}^{-1}\equiv 1/2$. The solution is easy to find, cf.
(6),

$\xi_{n}^{p}={(2n-1)\Gamma(n+1)\Gamma(1/2)\over 2\Gamma(p+3/2)\Gamma(n-p)}.$
(129)

To match the spectrum of individual bosons $\varepsilon_{n}(k,u_{0})$ with the
spectrum of quantum $Q_{2n-1}$ we need to restore powers of $\tilde{c}$ and
make the following identification

$Q_{2n-1}=\tilde{\Delta}^{n}+\tilde{c}^{n-1}\sum_{k}\left(n_{k}+{1\over
2}\right)\varepsilon_{n}(k,u_{0})+\dots$ (130)

where $n_{k}$ are boson occupation numbers of boundary gravitons and
$u_{0}=\tilde{\Delta}/\tilde{c}$. This reproduces the spectrum of $Q_{2n-1}$
at two first leading orders in $1/c$ and provides physical interpretation of
$n_{k}$. Unfortunately the holographic picture provides no clear path to
compute higher $1/c$ corrections to (130).

##  Appendix B Brute-force pertubative calculation

A straightforward but a laborious approach to evaluate $Q_{2n-1}$ in terms of
action variables $I_{k}$ would be to use Fourier modes $u_{k}$ of $u$,

$u(\varphi)=\sum_{k}u_{k}e^{ik\varphi},$ (131)

to parametrize the co-adjoint orbit of Virasoro algebra, i.e. the space of
potentials $u$ sharing the same orbit invariant $h$ (29). To that end $u_{0}$
should be understood as a function of $u_{k}$ Witten . Then $Q_{2n-1}$ and
$I_{k}$ can be expressed in terms of $u_{k}$, and consequently in terms of
each other.

In terms of the Fourier modes the Poisson bracket is

$i\frac{c}{24}\\{u_{k},u_{\ell}\\}=(k-\ell)u_{k+\ell}+2k^{3}\delta_{k+\ell}.$
(132)

This coincides with the Virasoro algebra upon $u_{0}$ is shifted by a
constant. At this point we introduce the orbit invariant $h(u_{k})$ and
express it in terms of $u_{k}$ by expanding in power series

$h=u_{0}+\sum_{n=2}^{\infty}U_{n},\qquad
U_{n}=\frac{1}{n!}\sum_{\begin{subarray}{c}p_{1},\cdots,p_{n}\\\
o_{1}+\dots+p_{n}=0,\,p_{k}\neq
0\end{subarray}}h_{p_{1},\cdots,p_{n}}u_{p_{1}}\cdots u_{p_{n}}.$ (133)

After imposing ${c\over 24}\\{h(u),u_{k}\\}=0$ for any $k$ we find

$\displaystyle h_{p_{1},p_{2}}$ $\displaystyle=$
$\displaystyle.-\frac{p_{1}^{2}+p_{2}^{2}+2h}{4(p_{1}^{2}+h)(p_{2}^{2}+h)},$
$\displaystyle h_{p_{1},p_{2},p_{3}}$ $\displaystyle=$
$\displaystyle\frac{p_{1}^{2}+p_{2}^{2}+p_{3}^{2}+6h}{8(p_{1}^{2}+h)(p_{2}^{2}+h)(p_{3}^{2}+h)},$
$\displaystyle h_{p_{1},p_{2},p_{3},p_{4}}$ $\displaystyle=$
$\displaystyle-\frac{15h^{4}-25h^{3}q_{2}+h^{2}(13q_{2}^{2}-9q_{4})+h(-3(q_{2}^{3}+q_{3}^{2})+8q_{2}q_{4})+q_{2}^{2}q_{4}-q_{2}q_{3}^{2}-4q_{4}^{2}}{den},\quad$

where

$\displaystyle q_{2}\equiv
p_{1}p_{2}+p_{1}p_{3}+p_{2}p_{3}+p_{1}p_{4}+p_{2}p_{4}+p_{3}p_{4}=-\frac{1}{2}\sum
p_{k}^{2},$ (134) $\displaystyle q_{3}\equiv
p_{1}p_{2}p_{3}+p_{1}p_{2}p_{4}+p_{1}p_{3}p_{4}+p_{2}p_{3}p_{4}=\frac{1}{3}\sum
p_{k}^{3},$ (135) $\displaystyle q_{4}\equiv
p_{1}p_{2}p_{3}p_{4}=-\frac{1}{4}\sum p_{k}^{4}+\frac{1}{2}q_{2}^{2},$ (136)
$\displaystyle
den=[(p_{1}+p_{2})^{2}+(p_{3}+p_{4})^{2}+2h][(p_{1}+p_{3})^{2}+(p_{2}+p_{4})^{2}+2h][(p_{1}+p_{4})^{2}+(p_{2}+p_{3})^{2}+2h]\times$
$\displaystyle\qquad\prod_{k=1}^{4}((p_{k}^{2}+h).$ (137)

Now we can get rid of $u_{0}=h-\sum_{n=2}^{\infty}U_{n}$ and express the
Poisson brackets in terms of $u_{k}$, $k\neq 0$,

$i\frac{c}{24}\\{u_{k},u_{\ell}\\}=\delta_{k+\ell}(2k)(k^{2}+h-\sum_{n=2}^{\infty}U_{n})+(1-\delta_{k+\ell})(k-\ell)u_{k+\ell}.$
(138)

Our next goal is to find symplectic form associated with the Poisson brackets

$\displaystyle\omega=\frac{c}{24}\times\frac{i}{2}\sum_{k\neq 0,\ell\neq
0}\omega_{k,\ell}\,du_{k}\wedge
du_{\ell},\qquad\omega_{k,\ell}=\sum_{n=0}^{\infty}\omega^{(n)}_{k,\ell},$
(139)

where $\omega^{(n)}_{k,\ell}$ is an order $n$ homogeneous polynomial in
$u_{k}$. We find, order by order in $u_{k}$,

$\displaystyle\omega^{(0)}_{k,\ell}$ $\displaystyle=$
$\displaystyle-\frac{1}{2k(k^{2}+h)}\delta_{k+\ell},$ (140)
$\displaystyle\omega^{(1)}_{k,\ell}$ $\displaystyle=$
$\displaystyle-(1-\delta_{k+\ell})\frac{k-\ell}{4k\ell(k^{2}+h)(\ell^{2}+h)}u_{-k-\ell},$
(141) $\displaystyle\omega^{(2)}_{k,\ell}$ $\displaystyle=$
$\displaystyle-\frac{1}{8k\ell(k^{2}+h)(\ell^{2}+h)}\sum_{m\neq
0}\left[\frac{k\delta_{k+\ell}}{m^{2}+h}u_{m}u_{-m}+\right.$ (143)
$\displaystyle\quad\left.(1-\delta_{k-m})(1-\delta_{\ell+m})\frac{(k+m)(\ell-m)}{m(m^{2}+h)}u_{-k+m}u_{-\ell-m}\right]$

We now would like to introduce (rescaled) normal coordinates $z_{k}$ near the
origin $u_{k}=0$ (which corresponds to constant $u(\varphi)=h$), such that

$\displaystyle\frac{24}{c}\omega$ $\displaystyle=\frac{i}{2}\sum_{k\neq
0}\frac{-1}{2k(k^{2}+h)}d{z}_{k}\wedge d{z}_{-k},\qquad
i\frac{c}{24}\\{{z}_{k},{z}_{\ell}\\}=\delta_{k+\ell}(2k)(k^{2}+h).$ (144)

We find

$\displaystyle{z}_{k}$
$\displaystyle=u_{k}+\frac{1}{4}\sum_{\begin{subarray}{c}p_{1}+p_{2}=k\\\
p_{i}\neq
0\end{subarray}}\frac{1}{p_{1}p_{2}}u_{p_{1}}u_{p_{2}}+\frac{1}{24}\sum_{\begin{subarray}{c}p_{1}+p_{2}+p_{3}=k\\\
p_{i}\neq
0,k\end{subarray}}\frac{p_{1}p_{2}+p_{2}p_{3}+p_{3}p_{1}}{p_{1}p_{2}p_{3}(p_{1}+p_{2})(p_{2}+p_{3})(p_{3}+p_{1})}u_{p_{1}}u_{p_{2}}u_{p_{3}}$
$\displaystyle\quad-\frac{1}{8}\sum_{\ell\neq 0,\pm
k}\frac{2\ell^{2}+h}{\ell^{2}(k^{2}-\ell^{2})(\ell^{2}+h)}u_{k}u_{\ell}u_{-\ell}-\frac{2k^{4}-k^{2}h+h^{2}}{32k^{4}(k^{2}+h)^{2}}u_{k}^{2}u_{-k}+\mathcal{O}(u^{4}).$
(145)

This expression can be inverted

$\displaystyle u_{k}$
$\displaystyle={z}_{k}-\frac{1}{4}\sum_{\begin{subarray}{c}p_{1}+p_{2}=k\\\
p_{i}\neq
0\end{subarray}}\frac{1}{p_{1}p_{2}}{z}_{p_{1}}{z}_{p_{2}}+\frac{1}{24}\sum_{\begin{subarray}{c}p_{1}+p_{2}+p_{3}=k\\\
p_{i}\neq
0,k\end{subarray}}\frac{k^{2}}{p_{1}p_{2}p_{3}(k-p_{1})(k-p_{2})(k-p_{3})}{z}_{p_{1}}{z}_{p_{2}}{z}_{p_{3}}$
$\displaystyle\quad+\frac{1}{2}\sum_{\begin{subarray}{c}p_{1}+p_{2}=0\\\
p_{i}\neq
0,k\end{subarray}}\frac{h}{p_{1}p_{2}(k-p_{1})(k-p_{2})[(p_{1}-p_{2})^{2}+4h]}{z}_{p_{1}}{z}_{p_{2}}{z}_{k}-\frac{h(5k^{2}+h)}{32k^{4}(k^{2}+h)^{2}}{z}_{k}^{2}{z}_{-k}+\mathcal{O}({z}^{4}).$
(146)

We are now ready to introduce action and angles variables $I_{k},\theta_{k}$
such that ${24\over c}\omega=\sum_{k}I_{k}\,dI_{k}\wedge d\theta_{k}$,

$\displaystyle\frac{{z}_{k}}{\sqrt{2k(k^{2}+h)}}=\sqrt{I_{k}}e^{-i\theta_{k}}.$
(147)

This leads to

$\displaystyle I_{k}=\frac{1}{2k(k^{2}+h)}{z}_{k}{z}_{-k}$ (148)
$\displaystyle=\frac{1}{2k(k^{2}+h)}u_{k}u_{-k}+\frac{1}{8k(k^{2}+h)}\sum_{\begin{subarray}{c}p_{1}+p_{2}=k\\\
p_{i}\neq
0\end{subarray}}\frac{1}{p_{1}p_{2}}u_{p_{1}}u_{p_{2}}u_{-k}+\frac{1}{8k(k^{2}+h)}\sum_{\begin{subarray}{c}p_{1}+p_{2}=-k\\\
p_{i}\neq
0\end{subarray}}\frac{1}{p_{1}p_{2}}u_{p_{1}}u_{p_{2}}u_{k}+\mathcal{O}(u^{4}).$

At this point we can go back to $u_{0}=h-\sum_{n=2}^{\infty}U_{n}$ are
represent it in terms of action variables (by expressing both sides as a
series in $z_{k}$),

$\displaystyle Q_{1}\equiv u_{0}=h+\sum_{k=1}kI_{k}+\mathcal{O}(z^{5}).$ (149)

This matches the exact relation (36) up to the fifth order in $z_{k}$,
reflecting the expansion order in (146).

To find $Q_{2n-1}$ in terms of $I_{k}$ we first write an iterative relation
for the Fourier modes of Gelfand-Dikii polynomials, which satisfy (32),

$\displaystyle R_{n,k}$ $\displaystyle\equiv$ $\displaystyle\frac{1}{2\pi}\int
d\varphi\,e^{-ik\varphi}R_{n},$ (150) $\displaystyle R_{n+1,k}$
$\displaystyle=$
$\displaystyle\frac{n+1}{2n+1}\left[2(k^{2}+u_{0})R_{n,k}+Q_{2n-1}u_{k}+\frac{1}{k}\sum_{\ell\neq
0,k}(2k-\ell)u_{\ell}R_{n,k-\ell}\right],\quad(k\neq 0).$

Then, using the relation between $Q_{2n-1}$ and $R_{n}$

$\displaystyle R_{n,k}=\frac{1}{ik(2n-1)}\frac{c}{24}\\{Q_{2n-1},u_{k}\\}$
(151)

we find

$\displaystyle\frac{c}{24}\\{Q_{2n+1},u_{k}\\}=ik(n+1)Q_{2n-1}u_{k}+\frac{2(n+1)(k^{2}+u_{0})}{2n-1}\frac{c}{24}\\{Q_{2n-1},u_{k}\\}+$
$\displaystyle\frac{n+1}{2n-1}\sum_{\ell\neq
0,k}\frac{2k-\ell}{k-\ell}u_{\ell}\frac{c}{24}\\{Q_{2n-1},u_{k-\ell}\\}.$
(152)

We use the following ansatz for $Q_{2n-1}$ in terms of $u_{k}$,

$\displaystyle Q_{2n-1}=$ $\displaystyle
h^{n}+\frac{1}{2!}\sum_{\begin{subarray}{c}p_{1}+p_{2}=0\\\ p_{i}\neq
0\end{subarray}}q^{(n)}_{p_{1},p_{2}}u_{p_{1}}u_{p_{2}}+\frac{1}{3!}\sum_{\begin{subarray}{c}p_{1}+p_{2}+p_{3}=0\\\
p_{i}\neq
0\end{subarray}}q^{(n)}_{p_{1},p_{2},p_{3}}u_{p_{1}}u_{p_{2}}u_{p_{3}}$ (153)
$\displaystyle+\frac{1}{4!}\sum_{\begin{subarray}{c}p_{1}+p_{2}+p_{3}+p_{4}=0\\\
p_{i}\neq
0\end{subarray}}q^{(n)}_{p_{1},p_{2},p_{3},p_{4}}u_{p_{1}}u_{p_{2}}u_{p_{3}}u_{p_{4}}+\mathcal{O}(u^{5}),$

and the iterative relation (152) becomes the iterative relation for
$q^{(n)}_{p_{1},\dots,p_{i}}$ for $i=2,3,4$,

$\displaystyle
q^{(n+1)}_{k,-k}=\frac{2(n+1)(k^{2}+h)}{2n-1}q^{(n)}_{k,-k}+\frac{(n+1)h^{n}}{2(k^{2}+h)},$
(154) $\displaystyle q^{(n+1)}_{p_{1},p_{2},p_{3}}$
$\displaystyle=\frac{2(n+1)(p_{1}^{2}+p_{2}^{2}+p_{3}^{2}+3h)}{3(2n-1)}q^{(n)}_{p_{1},p_{2},p_{3}}-\frac{(n+1)h^{n}(p_{1}^{2}+p_{2}^{2}+p_{3}^{2}+6h)}{8(p_{1}^{2}+h)(p_{2}^{2}+h)(p_{3}^{2}+h)}$
(155)
$\displaystyle\quad-\frac{n+1}{3(2n-1)}\left[\left(\frac{p_{1}-p_{2}}{p_{2}}+\frac{p_{1}-p_{3}}{p_{3}}\right)q^{(n)}_{p_{1},-p_{1}}+\text{symmetric
w.r.t. $p_{1},p_{2},p_{3}$}.\right]$ (156)

These can be solved as follows

$\displaystyle
q^{(n)}_{k,-k}=\frac{(h+k^{2})^{n-2}(2n)!!}{4(2n-3)!!}\sum_{m=0}^{n-1}\frac{(2m-1)!!}{(2m)!!}\left(\frac{h}{h+k^{2}}\right)^{m},$
(157) $\displaystyle q^{(n)}_{p_{1},p_{2},p_{3}}=$ (158)
$\displaystyle\frac{(2n)!!}{8(2n-3)!!p_{1}p_{2}p_{3}}\sum_{m=0}^{n-1}\frac{(2m-1)!!}{(2m)!!}h^{m}\left[p_{1}(p_{1}^{2}+h)^{n-m-2}+p_{2}(p_{2}^{2}+h)^{n-m-2}+p_{3}(p_{3}^{2}+h)^{n-m-2}\right].$

Our goal would be to match (153)with the expansion

$\displaystyle Q_{2n-1}$
$\displaystyle=h^{n}+\sum_{k=1}(f^{(n,1)}_{k}I_{k}+f^{(n,2)}_{k}I_{k}^{2})+\frac{1}{2}\sum_{\begin{subarray}{c}k,\ell=1\\\
k\neq\ell\end{subarray}}f^{(n)}_{k,\ell}I_{k}I_{\ell}+\mathcal{O}(I^{3}),$
(159)

by expressing $I_{k}$ in terms of $u_{k}$ using (148). This leads to

$\displaystyle f^{(n,1)}_{k}=2k(k^{2}+h)q^{(n)}_{k,-k}$ (160)

and the relations for $f^{(n,2)}_{k},f^{(n,2)}_{k\ell}$ in terms of
$q^{(n)}_{p_{1},\dots,p_{i}}$. To fix $f^{(2)}_{k},f^{(2)}_{k\ell}$, we would
not need $q^{(n)}_{p_{1},p_{2},p_{3},p_{4}}$ with arbitrary
$p_{1},\dots,p_{4}$, but only $q^{(n)}_{k,-k,\ell,-\ell}$, including the case
of $k=\ell$. The iterative relation for $q^{(n)}_{k,-k,\ell,-\ell}$ is
cumbersome. Instead, it is more convenient to work directly with the iterative
relation in terms of $f^{(2)}_{k}$ and $f^{(2)}_{k\ell}$. Once everything
combined together we find

$\displaystyle f^{(n,1)}_{k}$
$\displaystyle=\frac{(2n)!!k(h+k^{2})^{n-1}}{2(2n-3)!!}\sum_{m=0}^{n-1}\frac{(2m-1)!!}{(2m)!!}\left(\frac{h}{h+k^{2}}\right)^{m},$
(161) $\displaystyle f_{k}^{(n,2)}$
$\displaystyle=-\frac{(2n)!!(h+k^{2})^{n-2}}{16(2n-3)!!}\sum_{m=0}^{n-1}(3h+k^{2}-4k^{2}m)\sum_{j=0}^{m-1}\frac{(2j-1)!!}{(2j)!!}\left(\frac{h}{h+k^{2}}\right)^{j},$
(162)

and

$\displaystyle
f^{(n)}_{k,\ell}=\frac{(2n)!!k\ell}{4(2n-3)!!(k^{2}-\ell^{2})}\sum_{j=0}^{n-1}(n-1-j)\frac{(2j-1)!!h^{j}}{(2j)!!}\left[(h+k^{2})^{n-j-1}-(h+\ell^{2})^{n-j-1}\right]$
$\displaystyle+$
$\displaystyle\frac{(2n)!!k\ell}{4(2n-3)!!}\sum_{m=0}^{n-1}\sum_{j=0}^{m-1}\frac{m(2j-1)!!h^{j}}{(2j)!!}\left[(h+k^{2})^{m-j-1}(h+\ell^{2})^{n-m-1}+(h+\ell^{2})^{m-j-1}(h+k^{2})^{n-m-1}\right].$

Although written in a different form, this result is in agreement with (47),
(48), and (51).

##  Appendix C One-zone potentials: details

One-zone potentials $u$ can be found from the condition $\\{Q_{3}+\alpha
Q_{1},u\\}=0$ for some constant $\alpha$. From here we immediately find, see
section 2.4 of Dymarsky_2020 ,

$\displaystyle\lambda_{0}$
$\displaystyle=-\frac{\alpha}{24}-\frac{k^{2}}{12}(\theta_{3}(\tau)^{4}+\theta_{4}(\tau)^{4}),$
(164) $\displaystyle\lambda_{1}$
$\displaystyle=-\frac{\alpha}{24}-\frac{k^{2}}{12}(\theta_{2}(\tau)^{4}-\theta_{4}(\tau)^{4}),$
(165) $\displaystyle\lambda_{2}$
$\displaystyle=-\frac{\alpha}{24}+\frac{k^{2}}{12}(\theta_{2}(\tau)^{4}+\theta_{3}(\tau)^{4}).$
(166)

Pertubatively, i.e. in the limit of small $q=e^{i\pi\tau}$, corresponding
potential is

$\displaystyle
u=h+\frac{32k^{4}}{k^{2}+h}q^{2}-16k^{2}q\cos(k\varphi)-32k^{2}q^{2}\cos(2k\varphi)+\mathcal{O}(q^{3}).$
(167)

There are useful relations involving Jacobi elliptic functions and
hypergeometric function,

$\displaystyle m:=\theta_{2}^{4}(\tau)/\theta_{3}^{4}(\tau),\qquad
F\left(\frac{1}{2},\frac{1}{2},1;m\right)=\theta_{3}(\tau)^{2},\quad\frac{F\left(\frac{1}{2},\frac{1}{2},1;1-m\right)}{F\left(\frac{1}{2},\frac{1}{2},1;m\right)}=-\frac{1}{\pi}\log
q,$
$\displaystyle\frac{F\left(\frac{3}{2},\frac{1}{2},1;m\right)}{F\left(\frac{1}{2},\frac{1}{2},1;m\right)}=1+2{\partial\ln\theta_{3}^{2}(\tau)\over\partial\ln
m},\quad-16\sum_{n=0}^{\infty}\frac{q^{2n+1}}{(1-q^{2n+1})^{2}}+2\theta_{3}(\tau)^{4}-2\theta_{4}(\tau)^{4}\frac{F\left(\frac{3}{2},\frac{1}{2},1;m\right)}{F\left(\frac{1}{2},\frac{1}{2},1;m\right)}=0.$

We also give here more terms in the $q$-expansion of $I_{k}$,

$\displaystyle I_{k}=$
$\displaystyle\frac{32k^{3}q^{2}}{h+k^{2}}+\frac{64q^{4}\left(3h^{2}k^{3}+12hk^{5}+k^{7}\right)}{\left(h+k^{2}\right)^{3}}+\frac{128k^{3}q^{6}\left(3h^{4}+42h^{3}k^{2}+108h^{2}k^{4}-58hk^{6}+k^{8}\right)}{\left(h+k^{2}\right)^{5}}$
$\displaystyle+\frac{128k^{3}q^{8}\left(7h^{6}+156h^{5}k^{2}+1083h^{4}k^{4}+1232h^{3}k^{6}-4035h^{2}k^{8}+788hk^{10}+k^{12}\right)}{\left(h+k^{2}\right)^{7}}+\mathcal{O}(q^{10}),$

which with help of $Q_{1}=h+kI_{k}$ immediately yields

$\displaystyle Q_{1}$
$\displaystyle=h+\frac{32k^{4}}{k^{2}+h}q^{2}+\frac{64k^{4}(3h^{2}+12hk^{2}+k^{4})}{(k^{2}+h)^{3}}q^{4}$
(168)
$\displaystyle\quad+\frac{128k^{4}(3h^{4}+42h^{3}k^{2}+108h^{2}k^{4}-58hk^{6}+k^{8})}{(k^{2}+h)^{5}}q^{6}+\mathcal{O}(q^{7}).$
(169)

The relation for $I_{k}$ in terms of $q$ can be solved for $q$ in terms of
$I_{k}$ iteratively, which was used in section 2.2.

$\displaystyle q^{2}=$
$\displaystyle\frac{\left(h+k^{2}\right)}{32k^{3}}I_{k}-\frac{\left(3h^{2}+12hk^{2}+k^{4}\right)}{512k^{6}}I_{k}^{2}+\frac{\left(15h^{3}+87h^{2}k^{2}+105hk^{4}+k^{6}\right)}{8192k^{9}}I_{k}^{3}$
$\displaystyle-\frac{\left(187h^{4}+1402h^{3}k^{2}+3012h^{2}k^{4}+1606hk^{6}+k^{8}\right)}{262144k^{12}}I_{k}^{4}+\mathcal{O}(I_{k}^{5}).$

##  Appendix D Perturbative calculation for finite-zone potentials

We start with the two-zone case and parametrize corresponding differential
$dp$ with help of two infinitesimal parameters $\epsilon_{1},\epsilon_{2}$ and
$\lambda_{0}$,

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\lambda_{0}+{k^{2}\over
4}+\epsilon_{1}+a_{1}\epsilon_{1}^{2}+b_{1}\epsilon_{1}\epsilon_{2}+c1\epsilon_{2}^{2}+\dots,$
(170) $\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\lambda_{0}+{k^{2}\over
4}-a\epsilon_{1}+a_{2}\epsilon_{1}^{2}+b_{2}\epsilon_{1}\epsilon_{2}+c2\epsilon_{2}^{2}+\dots,$
(171) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\lambda_{0}+{\ell^{2}\over
4}+\epsilon_{2}+a_{3}\epsilon_{1}^{2}+b_{3}\epsilon_{1}\epsilon_{2}+c3\epsilon_{2}^{2}+\dots,$
(172) $\displaystyle\lambda_{4}$ $\displaystyle=$
$\displaystyle\lambda_{0}+{\ell^{2}\over
4}-b\epsilon_{2}+a_{4}\epsilon_{1}^{2}+b_{4}\epsilon_{1}\epsilon_{2}+c4\epsilon_{2}^{2}+\dots,$
(173) $\displaystyle r_{1}$ $\displaystyle=$
$\displaystyle\lambda_{0}+{k^{2}\over
4}+d_{1}\epsilon_{1}^{2}+e_{1}\epsilon_{1}\epsilon_{2}+f_{1}\epsilon_{2}^{2}+\dots,$
(174) $\displaystyle r_{2}$ $\displaystyle=$
$\displaystyle\lambda_{0}+{\ell^{2}\over
4}+d_{2}\epsilon_{1}^{2}+e_{2}\epsilon_{1}\epsilon_{2}+f_{2}\epsilon_{2}^{2}+\dots$
(175)

The parametrization is redundant, with different choices related by
redefinitions of $\epsilon_{1},\epsilon_{2}$. We assume
$\epsilon_{1}\sim\epsilon_{2}$ are of the same order and in what follows we
refer to expansion in $\epsilon_{1},\epsilon_{2}$ simply as $\epsilon$
expansion. While keeping two-zone case in mind for concreteness, most of the
discussion below applies to $m$-zone case with arbitrary $m$.

###  D.1 $a$-cycles

To impose $a_{1}$-cycle constraint (21), we need to integrate from
$\lambda_{1}$ to $\lambda_{2}$. By introducing $x$ via

$\displaystyle\lambda={\lambda_{2}+\lambda_{1}\over
2}+x{\lambda_{2}-\lambda_{1}\over 2},$ (176)

and then expanding in powers of $\epsilon$ we reduce the integral to standard
integrals of the form

$\displaystyle\int_{-1}^{1}{dx\,x^{2n}\over\sqrt{1-x^{2}}}={\sqrt{\pi}\Gamma(n+1/2)\over\Gamma(n+1)}.$
(177)

Provided we want to find $Q_{2n-1}$ in terms of $I_{k}$ by expanding up to
$p$-th power, we would need to keep $2p$ terms in $\epsilon$-expansion, up to
and including $\epsilon^{2p}$. This method works for any $a$-cycle integral
and any number of zones.

###  D.2 $b$-cycles

We start with the $b_{1}$-cycle, which goes from $\lambda_{0}$ to
$\lambda_{1}$, and introduce another variable $x$

$\displaystyle\lambda=\lambda_{1}-x(\lambda_{1}-\lambda_{0}).$ (178)

We can use the proximity of $\lambda_{4}$ to $\lambda_{3}$ to expand
$\sqrt{(\lambda-\lambda_{3})(\lambda-\lambda_{4})}$ in $\epsilon$. Now the
integral of interest reduced to a sum of integrals of the form

$\displaystyle\int_{0}^{1}{dxP(x)\over\sqrt{x(1-x)(16w+x)}(x-c)^{r}}$ (179)

where $16w$ is a small parameter of order $\epsilon$,

$\displaystyle 16w={\lambda_{2}-\lambda_{1}\over\lambda_{1}-\lambda_{0}},$
(180)

$P(x)$ is some polynomial and $c=1-\ell^{2}/k^{2}$ (we assumed $\ell>k$). The
integral (179) can be related to

$\displaystyle
J_{n}(c):=\int_{0}^{1}{dx\,x^{n}\over\sqrt{x(1-x)(16w+x)}(x-c)}$ (181)

by differentiating over $c$. To evaluate it, it is helpful to first introduce
the integral

$\displaystyle
I_{n}:=\int_{0}^{1}{dx\,x^{n}\over\sqrt{x(1-x)(16w+x)}}=\sum_{m=0}^{\infty}a_{m}(n)w^{m}+\sum_{m=n}^{\infty}b_{m}(n)w^{m}\ln
w,$ (182)

which can be expressed as formal series in $w$. Coefficients $a_{m}(n)$ for
$n>m$ and $b_{m}(n)$ for any $n,m$ can be found analytically

$\displaystyle
a_{m}(n)=\frac{(-16)^{m}\Gamma\left(m+\frac{1}{2}\right)\Gamma(n-m)}{\Gamma(m+1)\Gamma\left(-m+n+\frac{1}{2}\right)},\quad
n>m,$ (183) $\displaystyle
b_{m}(n)=\frac{16^{m}(-1)^{n+1}\Gamma\left(m+\frac{1}{2}\right)}{\Gamma(m+1)\Gamma(m-n+1)\Gamma\left(-m+n+\frac{1}{2}\right)}.$
(184)

To find $a_{m}(n)$ for $m\geq n$ we can use the iterative relation

$\displaystyle I_{n}=(1-2n)(I_{n}-I_{n-1})-{1\over
8}\partial_{w}(I_{n+1}-I_{n}),$ (185)

which follows from the integration by parts, and $a_{m}(0)$ which can be found
directly from (182) since the corresponding integral can be evaluated
analytically. For example we find the following iterative relation for
$a_{n}(n)$,

$\displaystyle
a_{m+1}(m+1)=\frac{(-1)^{m}2^{2m+3}\Gamma(2m+1)}{\Gamma(m+2)^{2}}-\frac{8(2m+1)a_{m}(m)}{m+1},\quad
a_{0}(0)=0.$ (186)

So far we are interested only in first $2p$ powers of $w$, we only need to
worry about $a_{m}(n)$ with $m\leq 2p$. In our case $p=3$ and we simply
tabulate values of $a_{m}(n)$ for $0\leq m\leq 6$ and $m\geq n$ for
convenience

$\displaystyle a_{m}(n)=\left(\begin{array}[]{ccccccc}0&&&&&&\\\ 8&8&&&&&\\\
-84&-104&-112&&&&\\\ \frac{2960}{3}&1152&\frac{4288}{3}&\frac{4736}{3}&&&\\\
-\frac{37310}{3}&-\frac{42040}{3}&-16368&-\frac{60992}{3}&-\frac{68224}{3}&&\\\
\frac{820008}{5}&180656&203584&\frac{1189248}{5}&\frac{1478144}{5}&\frac{1666048}{5}&\\\
-\frac{11153912}{5}&-\frac{12097344}{5}&-\frac{7995904}{3}&-\frac{9011456}{3}&-\frac{17549824}{5}&-\frac{65468416}{15}&-\frac{74166272}{15}\\\
\end{array}\right).$ (194)

Going back to (181), we can expand $(x-c)$ in the denominator into power
series in $x$, thus reducing the integral to a sum of (182). Provided $n>2p$
and so far we are only interested in terms of order $w^{r}$ and $w^{r}\ln w$
with $r\leq 2p$, only relevant contributions would come from $a_{m}(n)w^{m}$
term in (182) with $m<n$. Corresponding coefficients are known analytically,
(183), and can be re-summed yielding,

$\displaystyle
J_{n}(c)=-\sum_{m=0}^{2p}\frac{(-16)^{m}\omega^{m}\Gamma\left(m+\frac{1}{2}\right)\Gamma(l-m)\,_{2}F_{1}\left(1,l-m;l-m+\frac{1}{2};\frac{1}{c}\right)}{c\,\Gamma(m+1)}+{\mathcal{O}}(w^{2p+1}).$

Here $2F_{1}$ is regularized hypergeometric function and this expression is
only valid for $n>2p$. To extend it to smaller $n$ we use the iterative
relation, which follows from the integration by parts,

$\displaystyle J_{n}={J_{n+1}-I_{n}\over c}.$ (195)

This completes technical preliminaries as now integral over $b_{1}$ cycle can
be reduced to a number of integrals $J_{n}$ and their derivatives, so far we
are only interested in terms of order $w^{r}$ with $r\leq 2p$. Clearly, the
approach above can be used to evaluate integrals over $b_{1}$ when there are
more than two zones. In this case one would need to evaluate integrals

$\displaystyle\int_{0}^{1}{dx\,x^{n}\over\sqrt{x(1-x)(16w+x)}\prod_{i=1}^{m-1}(x-c_{i})},$
(196)

where $m$ is the number of zones. This can be reduced to (181) by noting

$\displaystyle\prod_{i=1}^{m-1}{1\over(x-c_{i})}=\sum_{i=1}^{m-1}{\alpha_{i}\over
x-c_{i}},$ (197)

with the appropriate coefficients $\alpha_{i}$.

To evaluate the integral over $b_{2}$-cycle from $\lambda_{2}$ to
$\lambda_{3}$ is more challenging because in the $\epsilon\rightarrow 0$ limit
there are singularities appearing at both boundaries. There is a
straightforward but complicated way. By appropriately changing variables and
expanding in $\epsilon$ all terms except for
$\sqrt{(\lambda-\lambda_{1})(\lambda-\lambda_{2})(\lambda-\lambda_{3})(\lambda-\lambda_{4})}$
we reduce the calculation to the integral

$\displaystyle\int_{0}^{1}dx{x^{n}\over\sqrt{x(1-x)(16w+x)(1+16u-x)}}$ (198)

for positive small $w,u$. The indefinite integral of this kind can be
evaluated analytically. Then the definite integral above can be integrated by
expanding it powers of $w,u$ (which both are of order $\epsilon$), and keeping
terms up to order $2p$. This is an involved exercise and instead one can use
one of the following shortcuts.

In the particular case of two-zone potential, instead of evaluating integral
over $b_{2}$, one can combine the integral over $b_{1}$ and $b_{2}$ such that
the contour would enclose $\lambda_{0},\dots,\lambda_{3}$. Now one can deform
the contour to go from $\lambda_{4}$ to infinity, if necessary accompanied by
a circle at infinity. At this point integrand can be expanded in $\epsilon$
such that brunch-cut from $\lambda_{1}$ to $\lambda_{2}$ disappears, yielding
pole singularities at $\lambda=\lambda_{0}+{k^{2}/4}$. At this point
corresponding integral can be rewritten as

$\displaystyle\oint_{-\infty}^{-16w}dx{P(x)\over\sqrt{x(1-x)(x+16w)}(x-c)^{r}},$
(199)

where $P(x)$ is some polynomial and $0\leq c\leq 1$. We also emphasize that to
render this integral finite, one may need to close the contour at infinity.
This integral can be decomposed into a sum of integrals of the form

$\displaystyle\oint_{-\infty}^{-16w}dx{x^{n}\over\sqrt{x(1-x)(x+16w)}},$ (200)

and

$\displaystyle\oint_{-\infty}^{-16w}dx{1\over\sqrt{x(1-x)(x+16w)}(x-c)},$
(201)

and its derivatives. First integral can be reduced to (182) by deforming the
contour to go from $0$ to $1$. Last integral can be reduced to $J_{0}$ and
$J_{1}$ with help of modular transformation mapping $\infty$ to $1$, $-16w$ to
$0$, and $0$ to $-16w$.

$\displaystyle x\rightarrow{x+16w\over x-1}.$ (202)

This shortcut works for two-zone case, but with more zones present it is not
applicable. Neveftheless there is a very simple trick which make evaluation of
$b_{2}$ and other $b$-cycles unnesessary. Indeed, to satisfy (21) and (22) for
all cycles, it is sufficient to satisfy (21) for all cycles and (22) for
$b_{1}$ and also impose that the expansion (170-175), and its generalizations
for the case of more than two zones, is invariant under permutation of indexes
and $k_{i}$ defined in (23). Say, for two zones we find

$\displaystyle\lambda_{1}$ $\displaystyle=\lambda_{0}+{k^{2}\over
4}-\epsilon_{1}+{3\epsilon_{1}^{2}\over
k^{2}}+{4\epsilon_{2}^{2}k^{2}\over\ell^{2}(k^{2}-\ell^{2})}+{\mathcal{O}}(\epsilon^{3}),$
(203) $\displaystyle\lambda_{2}$ $\displaystyle=\lambda_{0}+{k^{2}\over
4}+\epsilon_{1},$ (204) $\displaystyle r_{1}$
$\displaystyle=\lambda_{0}+{k^{2}\over 4}+{\epsilon_{1}^{2}\over
2k^{2}}+{2\epsilon_{2}^{2}k^{2}\over\ell^{2}(k^{2}-\ell^{2})}+{\mathcal{O}}(\epsilon^{3}).$
(205)

and $\lambda_{3,4},r_{2}$ related to $\lambda_{1,2},r_{1}$ by the exchange
$\epsilon_{1}\leftrightarrow\epsilon_{2}$ and $k\leftrightarrow\ell$. The same
logic with the permutation symmetry works for any number of zones.

Above we only explicitly wrote terms up to $\epsilon^{2}$, while evaluating
all terms up to $\epsilon^{6}$. The simple form of $\lambda_{2}$ above is a
parametrization choice. With this choice taking $\epsilon_{2}=0$ does not
close the second zone. One can check that taking

$\displaystyle\epsilon_{2}=-\frac{2\ell^{2}\epsilon_{1}^{2}}{k^{2}\left(k^{2}-\ell^{2}\right)}+\dots$
(206)

such that $\lambda_{4}=\lambda_{3}$ would make $I_{\ell}$ dicussed below
vanish. Alternatively one could chooe $\epsilon_{i}$ to control the size of
$\lambda_{2i}-\lambda_{2i-1}$, but with this choice both all $\lambda_{i}$
would depend on all $\epsilon_{i}$.

###  D.3 Evaluation of $I_{k}$, $h$ and $Q_{2n-1}$.

Evaluation of action variables $I_{k}$ as a pertubative series in
$\epsilon_{i}$ is straightforward. It is an integral over $a$-cycle and
therefore can be evaluated along the lines discussed above. The only
difference, in comparision with the discussion in subsection D.1, is the term
$\ln\lambda$, which needs to be expanded in powers of $\epsilon$ yielding
polynomials in $x$ in the numerator of (177),

$\displaystyle I_{k}={2\epsilon_{1}^{2}\over
k(\lambda_{0}+k^{2}/4)}+{\mathcal{O}}(\epsilon^{3}).$ (207)

Again, we only keep terms up to $\epsilon^{2}$ for simplicity.

Evaluation of $h$ is also straightforward. To that end one needs to calculate
$p(0)$, given by an integral from $0$ to $\lambda_{0}$. After expanding the
integrand in powers of $\epsilon$ it becomes the integral which can be
evaluated in a closed form, yielding

$\displaystyle h/4=\lambda_{0}+\lambda_{0}\left({2\epsilon_{1}^{2}\over
k^{2}(\lambda_{0}+k^{2}/4)}+{2\epsilon_{2}^{2}\over\ell^{2}(\lambda_{0}+\ell^{2}/4)}\right)+{\mathcal{O}}(\epsilon^{3}).$
(208)

Finally, evaluation of $Q_{2n-1}$ for any given $n$ is also straightforward
since $\lambda_{i}$ are known explicitly. As a result we obtain
$I_{k},h,Q_{2n-1}$ as functions of $\lambda_{0}$ and $\epsilon_{i}$. One can
then reverse-engineer coefficients in (13) such that it is satisfied.

##  Appendix E Spectrum of $Q_{2n-1}$ acting on primaries

In this appendix we outlined calculation of ${\sf Q}^{0}_{2n}$ (64) following
dorey2020geometric . Starting from the Schr$\ddot{\rm o}$dinger equation (60),
one introduces the following change of variables

$\displaystyle\Psi(x)=E^{l(l-3/2)/4\alpha}w^{(l-3/2)/4}y(w),\qquad x=E^{1\over
2\alpha}w^{1\over 2l},$ (209)

such that (60) becomes

$\displaystyle-\epsilon^{2}\partial_{w}^{2}y+Z(w)y=0,\qquad
Z(w)={w},\quad\epsilon=E^{-{\alpha+1\over 2\alpha}}.$ (210)

Taking $\epsilon$ as a formal small parameter this equation can be solved via
WKB expansion,

$\displaystyle y(w)=e^{{1\over\epsilon}S(w)},\qquad-\epsilon
S^{\prime\prime}-S^{\prime 2}+Z=0,\qquad
S(w)=\sum_{n=0}^{\infty}\epsilon^{n}S_{n}.$ (211)

The resulting Riccati equation can be rewritten as the iterative relation to
find $S^{\prime}_{n}$ with $S_{0}^{\prime}=-\sqrt{Z(w)}$. It is more
convenient for what follows to make another change of variables
$z=w^{\alpha/(l+1/2)}$ and introduce the polynomial ansatz

$\displaystyle S^{\prime}_{n}=i{\alpha\over
2l+1}z^{1-{l+1/2\over\alpha}}\tilde{S}_{n},\qquad\tilde{S}_{n}=\sum_{k=0}^{n}i^{n}c_{k}^{(n)}z^{-k+(n-1)(1-1/2\alpha)}(1-z)^{k-(3n-1)/2}$
(212)

The Ricatti equation rewritten in terms of $c_{k}^{(n)}$ gives rise to (62),
which can be used together with (63), to iteratively find $c^{(n)}_{k}$. The
first few $c_{k}^{(n)}$ read

$\displaystyle c_{0}^{(2)}=\frac{5}{8}\alpha,\quad
c_{1}^{(2)}=\frac{1}{4}(2\alpha-1),\quad
c_{2}^{(2)}=-\frac{1}{8\alpha}(4u^{2}\alpha^{2}-1).$ (213) $\displaystyle
c_{0}^{(3)}=-\frac{15\alpha^{2}}{8},\quad
c_{1}^{(3)}=-\frac{9\alpha^{2}}{4}+\frac{9\alpha}{8},\quad
c_{2}^{(3)}=\frac{1}{2}\alpha^{2}\left(u^{2}-1\right)+\frac{3\alpha}{4}-\frac{3}{8},\quad
c_{3}^{(3)}=-\frac{1}{8\alpha}(4u^{2}\alpha^{2}-1).$ (214)

To obtain ${\sf Q}^{0}_{2n-1}$ one needs to integrate $\tilde{S}_{2n}(w(z))$
over a Pochhammer contour $\gamma_{P}$,

$\displaystyle{\sf Q}^{0}_{2n-1}$ $\displaystyle=$
$\displaystyle(-1)^{n}\frac{\Gamma\left(\frac{3}{2}-n-\frac{2n-1}{2\alpha}\right)}{\sqrt{\pi}\Gamma\left(1-\frac{2n-1}{2\alpha}\right)}\frac{(2n-1)\Gamma(n+1)}{4^{n}(\alpha+1)^{n}}\check{I}_{2n-1}(\alpha,\hat{l}),$
(215) $\displaystyle\check{I}_{2n-1}$ $\displaystyle=$
$\displaystyle\frac{1}{2(1-e^{\frac{-i\pi(2n-1)}{\alpha}})}\int_{\gamma_{P}}dz\,\tilde{S}_{2n}(z).$
(216)

This integral can be evaluated using,

$\displaystyle(1-e^{2\pi ia})(1-e^{2\pi
ib})B(a,b)=\int_{\gamma_{P}}dz\,z^{a-1}(1-z)^{b-1},$ (217)

where $B(a.b)$ is the Euler beta function. Combining everything together
yields (64).

Evaluating ${\sf Q}^{0}_{2n-1}$ explicitly, using computer algebra to solve
for $c_{k}^{(n)}$ iteratively, for small and moderate $n$ is an easy task. To
obtain $1/c$ expansion of $\lambda^{0}_{2n-1}$ for arbitrary $n$ requires
knowing corresponding $c_{k}^{(n)}$ in $1/c$ expansion, i.e. in the limit of
large $\alpha$. This proved to be a difficult task. We obtained first three
non-trivial terms of $\lambda^{0}_{2n-1}$ in $1/\tilde{c}$ expansion (65),
with the first two terms (66,67) in closed analytical form. Functions $y_{i}$
and $\zeta_{i}$ there are defined as follows

$\displaystyle y_{1}(j)$ $\displaystyle=$
$\displaystyle\sum_{\ell=0}^{j}{1\over 2\ell+1},$ (218) $\displaystyle
y_{2}(j)$ $\displaystyle=$
$\displaystyle\sum_{\ell=0}^{j}{1\over(2\ell+1)^{2}},$ (219)
$\displaystyle\zeta_{2}(j)$ $\displaystyle=$
$\displaystyle\sum_{j_{1}+j_{2}=j}\zeta(-2j_{1}-1)\zeta(-2j_{2}-1),$ (220)
$\displaystyle\zeta_{3}(j)$ $\displaystyle=$
$\displaystyle\sum_{j_{1}+j_{2}+j_{3}=j}\zeta(-2j_{1}-1)\zeta(-2j_{2}-1)\zeta(-2j_{3}-1),$
(221)

where sum goes only over non-negative $j_{1},j_{2},j_{3}$. Third term (68) was
fixed up to one coefficient $p_{j}$, with the first several values for $0\leq
j\leq 17$ given below

$\displaystyle
p_{j}=\left(-\frac{31}{224},\frac{103}{576},-\frac{7883}{21120},\frac{868487}{748800},-\frac{505639}{100800},\frac{394694297}{13708800},-\frac{68117454019}{321753600},\frac{4929720750223}{2540160000},\right.$
$\displaystyle\left.-\frac{199232137825687}{9180864000},\frac{48745030162337923}{167650560000},-\frac{618684597383137}{134534400},\frac{7442737871872435019}{87783696000},\right.$
$\displaystyle\left.-\frac{1420749127340184137621}{788237049600},\frac{46636700018927407368821}{1065512448000},-\frac{198277953077778046100039}{164670105600},\right.$
$\displaystyle\left.\frac{21869843836862719834306038469}{587058612940800},-\frac{31428771773709445918185916879}{24404109649920},\right.$
$\displaystyle\left.\frac{4187283526052269558397574465940213}{84663488093184000},\dots\right).$

## References

  * (1) A. A. Belavin, A. M. Polyakov and A. B. Zamolodchikov, _Infinite conformal symmetry in two-dimensional quantum field theory_ , _Nuclear Physics B_ 241(2) (1984) 333.
  * (2) A. Zamolodchikov, _Conformal field theory and critical phenomena in two dimensional systems_ , vol. 10. CRC Press, 1989.
  * (3) T. Hartman, C. A. Keller and B. Stoica, _Universal spectrum of 2d conformal field theory in the large c limit_ , _Journal of High Energy Physics_ 2014 (2014) 1.
  * (4) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Integrable structure of conformal field theory, quantum kdv theory and thermodynamic bethe ansatz_ , _Communications in Mathematical Physics_ 177 (1996) 381.
  * (5) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Integrable structure of conformal field theory ii. q-operator and ddv equation_ , _Communications in Mathematical Physics_ 190 (1997) 247.
  * (6) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Integrable structure of conformal field theory iii. the yang–baxter relation_ , _Communications in mathematical physics_ 200 (1999) 297.
  * (7) M. Srednicki, _Chaos and quantum thermalization_ , _Physical Review E_ 50 (1994) 888.
  * (8) N. Lashkari, A. Dymarsky and H. Liu, _Eigenstate thermalization hypothesis in conformal field theory_ , _Journal of Statistical Mechanics: Theory and Experiment_ 2018 (2018) 033101.
  * (9) F.-L. Lin, H. Wang and J.-j. Zhang, _Thermality and excited state rényi entropy in two-dimensional cft_ , _Journal of High Energy Physics_ 2016 (2016) 116.
  * (10) P. Basu, D. Das, S. Datta and S. Pal, _Thermality of eigenstates in conformal field theories_ , _Physical Review E_ 96 (2017) 022149\.
  * (11) S. He, F.-L. Lin and J.-j. Zhang, _Dissimilarities of reduced density matrices and eigenstate thermalization hypothesis_ , _JHEP_ 12 (2017) 073 [1708.05090].
  * (12) S. He, F.-L. Lin and J.-j. Zhang, _Subsystem eigenstate thermalization hypothesis for entanglement entropy in CFT_ , _JHEP_ 08 (2017) 126 [1703.08724].
  * (13) N. Lashkari, A. Dymarsky and H. Liu, _Universality of quantum information in chaotic cfts_ , _Journal of High Energy Physics_ 2018 (2018) 70.
  * (14) W.-z. Guo, F.-L. Lin and J. Zhang, _Note on ETH of descendant states in 2D CFT_ , _Journal of High Energy Physics_ 2019 (2019) .
  * (15) A. Maloney, G. S. Ng, S. F. Ross and I. Tsiares, _Generalized gibbs ensemble and the statistics of kdv charges in 2d cft_ , _Journal of High Energy Physics_ 2019 (2019) .
  * (16) A. Dymarsky and K. Pavlenko, _Generalized gibbs ensemble of 2d cfts at large central charge in the thermodynamic limit_ , _Journal of High Energy Physics_ 2019 (2019) 98.
  * (17) A. Dymarsky and K. Pavlenko, _Generalized eigenstate thermalization hypothesis in 2d conformal field theories_ , _Physical review letters_ 123 (2019) 111602.
  * (18) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Higher-level eigenvalues of q-operators and schroedinger equation_ , _Adv.Theor.Math.Phys._ (2004) .
  * (19) A. V. Litvinov, _On spectrum of ilw hierarchy in conformal field theory_ , _Journal of High Energy Physics_ 2013 (2013) .
  * (20) E. Witten, _Coadjoint Orbits of the Virasoro Group_ , _Commun. Math. Phys._ 114 (1988) 1.
  * (21) E. M. Brehm and D. Das, _Korteweg–de vries characters in large central charge CFTs_ , _Physical Review D_ 101 (2020) .
  * (22) A. Dymarsky and K. Pavlenko, _Exact generalized partition function of 2d cfts at large central charge_ , _Journal of High Energy Physics_ 2019 (2019) 77.
  * (23) S. P. Novikov, _The periodic problem for the korteweg–de vries equation_ , _Funktsional’nyi Analiz i ego Prilozheniya_ 8 (1974) 54.
  * (24) A. Dymarsky and S. Sugishita, _KdV-charged black holes_ , _Journal of High Energy Physics_ 2020 (2020) .
  * (25) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Quantum field theories in finite volume: Excited state energies_ , _Nuclear Physics B_ 489 (1997) .
  * (26) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Spectral determinants for schrödinger equation and q-operators of conformal field theory_ , _Journal of Statistical Physics_ 102 (2001) 567.
  * (27) V. V. Bazhanov, S. L. Lukyanov and A. B. Zamolodchikov, _Higher-level eigenvalues of q-operators and schrödinger equation_ , _arXiv preprint hep-th/0307108_ (2003) .
  * (28) P. Dorey, C. Dunning, S. Negro and R. Tateo, _Geometric aspects of the ode/im correspondence_ , _Journal of Physics A: Mathematical and Theoretical_ 53 (2020) 223001.
  * (29) R. Conti and D. Masoero, _On solutions of the Bethe Ansatz for the Quantum KdV model_ , 2112.14625.
  * (30) A. Dymarsky, K. Pavlenko and D. Solovyev, _Zero modes of local operators in 2d CFT on a cylinder_ , _JHEP_ 07 (2020) 172.
  * (31) A. Maloney, G. S. Ng, S. F. Ross and I. Tsiares, _Thermal correlation functions of kdv charges in 2d cft_ , _Journal of High Energy Physics_ 2019 (2019) .
  * (32) J. de Boer and D. Engelhardt, _Remarks on thermalization in 2d cft_ , _Physical Review D_ 94 (2016) 126019.
  * (33) M. Rigol, V. Dunjko, V. Yurovsky and M. Olshanii, _Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons_ , _Physical review letters_ 98 (2007) 050405.
  * (34) A. Pérez, D. Tempo and R. Troncoso, _Boundary conditions for general relativity on ads3 and the kdv hierarchy_ , _Journal of High Energy Physics_ 2016 (2016) 103.
  * (35) M. Downing and G. M. T. Watts, _Free fermions, KdV charges, generalised Gibbs ensembles and modular transforms_ , 2111.13950.
  * (36) E. Ojeda and A. Pérez, _Boundary conditions for general relativity in three-dimensional spacetimes, integrable systems and the kdv/mkdv hierarchies_ , _JHEP_ 08 (2019) 079 [1906.11226].