ODE/IQFT correspondence for the generalized affine sl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathfrak{sl} $$\end{document}(2) Gaudin model

An integrable system is introduced, which is a generalization of the sl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathfrak{sl} $$\end{document}(2) quantum affine Gaudin model. Among other things, the Hamiltonians are constructed and their spectrum is calculated using the ODE/IQFT approach. The model fits into the framework of Yang-Baxter integrability. This opens a way for the systematic quantization of a large class of integrable non-linear sigma models. There may also be some interest in terms of Condensed Matter applications, as the theory can be thought of as a multiparametric generalization of the Kondo model.


Introduction
Suppose we are a given a set of quantum spins S (a) = S with a = 1, . . . , r. As was pointed out by Gaudin [1,2], the operators mutually commute for an arbitrary choice of the parameters {z a } r a=1 . They also commute with any projection of the total spin operator S = r a=1 S (a) and, furthermore, S 2 is linearly expressed through the Hamiltonians: z a H (a) + r a=1 j a (j a + 1) . (1. 3) It turns out that the problem of the simultaneous diagonalization of H (a) can be "solved" within the framework of the Bethe ansatz [1,2]. In this approach the energy E a of the a-th Hamiltonian is expressed as The integer M + takes all possible values from 0 to 2 r a=1 j a . The Gaudin model admits an almost straightforward generalization to any simple Lie algebra g (see section 13.2.2 in [2] and ref. [3]). The development of the mathematical apparatus of 2D Conformal Field Theory led to the idea that there should be a meaningful generalization to the case when the finite-dimensional Lie algebra is replaced by an affine Kac-Moody algebra g. Then the Hilbert space would be built out of Verma modules for an JHEP09(2021)201 algebra of extended conformal symmetry. According to general principles of integrability in CFT [4], the diagonalization problem would be formulated for an infinite set {I s } of socalled local Integrals of Motion (IM), which depend on the arbitrary parameters {z a } r a=1 . These would mutually commute, while by local what is meant is that where T s+1 is a chiral local field of Lorentz spin s + 1. An important step in exploring the affine Gaudin model was made by Feigin and Frenkel in ref. [5]. The key idea came from an interesting link between the original Gaudin model described above and a certain class of linear differential equations of the form [2,6,7] − ∂ 2 z + t 0 (z) Ψ = 0 with t 0 (z) = r a=1 j a (j a + 1) (1.7) The ODE possesses r regular singular points at z = z a . If we further assume that r a=1 E a = 0 (1. 8) then z = ∞ is also a regular singularity so that (1.7) is a Fuchsian differential equation. A remarkable phenomenon occurs when the residues {E a } r a=1 in the potential t 0 (z) coincide with the set of energies corresponding to some common eigenvector of the Hamiltonians H (a) . In this case all the singular points at z = z a turn out to be apparent. 1 The roots x ( while (1.12)

JHEP09(2021)201
Thus there is a link between the spectrum of the Gaudin Hamiltonians and a class of differential equations possessing certain monodromy properties. This provides, perhaps, one of the simplest illustrations of a broad phenomena, known as the ODE/IQFT correspondence [8][9][10][11]. 2 In ref. [5] Feigin and Frenkel introduce the Hamiltonians, which can be interpreted as an "affinization" of H (a) (1.2). They are built from r independent copies of the affine Kac-Moody sl ka (2) Feigin and Frenkel put forward the conjecture that the spectrum of these operators would be encoded in a class of differential equations that generalizes (1.7), though they did not explain exactly how the spectrum would be extracted from the ODEs. The last point was clarified in refs. [12,13].
In this work we introduce and study the model, which possesses an infinite set of mutually commuting local integrals of motion. The simplest ones, the "Hamiltonians", are expressed in terms of the sl(2) Kac-Moody currents and Virasoro field (1.14) as H (a) gen = 2π 0 du 2π (1.17) 2 The term ODE/IM correspondence, where IM stands for integrable model, is also used in the literature.
However, we find that the latter tends to be heavily focused on formal properties of ODEs at the expense of a clear study of the quantum theory which, in our opinion, is the most interesting part of the relation.
Recently, there has been a great flux of works that introduce and discuss different sorts of "correspondences", which are essentially ODE/IM. In order to emphasize that our main motivation is the study of integrable quantum field theory, rather than abstract "integrable models", we use the term ODE/IQFT.

JHEP09(2021)201
Notice that 1 − β 2 β 2 K r a=1 H (a) gen = r a=1 2π 0 du 2π G (a) , (1.18) which can be thought of as the affine counterpart of eq. (1.3). The operators (1.16) depend on the parameter β. The Hamiltonians of the affine Gaudin model (1.15) are obtained through a certain limiting procedure, which includes taking β → 1 − . We formulate the ODE/IQFT correspondence for the model and explain how the spectrum of H (a) gen for arbitrary β ∈ (0, 1) can be extracted from the differential equations. The theory will be referred to as the Generalized Affine Gaudin Model (GAGM).
The GAGM fits within the framework of the standard Yang-Baxter integrability. In particular, the Hamiltonians H (a) gen are part of a large commuting family which, as usual, involves the quantum transfer-matrices and Baxter Q-operators. These are explicitly constructed along the lines of the BLZ approach [14][15][16]. It is proposed that the GAGM governs the critical behaviour of a lattice system, which in the simplest case coincides with the inhomogeneous six-vertex model introduced by Baxter [17]. The local Boltzmann weights are contained in the R-matrix that is the trigonometric solution of the Yang-Baxter equation. The anisotropy parameter entering into the R-matrix, commonly referred to as q, is related to the parameter β in (1.16) as (1. 19) In the limit β → 1 − the trigonometric R-matrix becomes the rational one. The paper is organized as follows. Sections 2-5 contain no new material and their purpose is to illustrate the ideas of the BLZ approach. Its main ingredient is a realization of the Borel subalgebra of the quantum algebra U q sl (2) in terms of the vertex operators. We use the example from ref. [18] as it contains the essential blueprints for the construction of the commuting family of operators in the GAGM as well as the ODE/IQFT correspondence. The realization of the Borel subalgebra of U q sl (2) , which gives rise to the commuting family for the GAGM, is presented in section 6. The next section is a central one and contains a detailed discussion of the ODE/IQFT correspondence for the model. Some comments concerning the literature are also presented therein. The way the spectrum of the local and non-local integrals of motion are extracted from the ODE is described in section 8. Sections 9 and 10 deal with some specific cases and provide an illustration of the rather abstract ideas that preceded them. The Hamiltonians H (a) gen (1.16) are deduced in section 11. Also considered are certain limits of the model such as the isotropic, classical as well as the limit, which yields the Hamiltonians of the affine Gaudin model (1.15). Finally, we sketch how the GAGM appears in the scaling limit of the Baxter-type statistical systems in section 12.

Quantum transfer-matrices
The algebraic structure underlying the Yang-Baxter relation was clarified within the theory of quasi-triangular Hopf algebras by Drinfeld [19]. A basic example is when the Hopf JHEP09(2021)201 algebra is U q ( g) -the quantum deformation of the universal enveloping algebra of the affine algebra [19,20]. The central rôle is played by the universal R-matrix, which lies in the tensor product U q ( g) ⊗ U q ( g) and satisfies the relation R 12 R 13 R 23 = R 23 R 13 R 12 . (2.1) An important feature of R is that it is decomposed as stand for the Borel subalgebras of U q ( g). In this paper we restrict to the case g = sl (2). Let us consider the evaluation homomorphism of U q ( g) to the loop algebra U q (g)[λ, λ −1 ] and specify a finite dimensional matrix representation π of U q (g). In the case under consideration the Borel subalgebra U q ( b + ) is generated by four elements, {y 0 , y 1 , h 0 , h 1 }, and its evaluation homomorphism is defined by Here h, e ± are the generators of U q sl(2) , subject to the commutation relations is a U q ( b − )-valued (2 +1)×(2 +1) matrix whose entries depend on an auxiliary parameter λ. In turn the formal algebraic relation (2.1) becomes the Yang-Baxter algebra R , λ 1 /λ 2 L (λ 1 ) ⊗ 1 1 ⊗ L (λ 2 ) = 1 ⊗ L (λ 2 ) L (λ 1 ) ⊗ 1 R , λ 1 /λ 2 (2.5) with R 1 , 2 (λ 1 /λ 2 ) = π 1 (λ 1 ) ⊗ π 2 (λ 2 ) [R] .
As an immediate consequence the operators τ (λ) = Tr q 1 2 h h 0 L (λ) , (2.6) usually referred to as the transfer-matrices, obey the commutativity condition τ (λ), τ (λ ) = 0 . (2.7) With the expression for the universal R-matrix given in [21], one can obtain L (λ) as a formal series expansion in powers of the spectral parameter λ. The first few terms read as 2 (2.8)

JHEP09(2021)201
Here and below, abusing notation, we do not distinguish between the formal generators of U q (sl 2 ) and their (2 + 1) × (2 + 1) matrices in a finite dimensional representation π . Also, [n] q ≡ (q n − q −n )/(q − q −1 ). The expression in the square brackets contains the elements x 0 , x 1 ∈ U q ( b − ), which obey the quantum Serre relations There are two remaining generators h 0 , h 1 satisfying (2.10) Since h 0 + h 1 is a central element, for our purposes and without loss of generality we have set it to be zero.
Up to this point, there was no need to specify a representation of U q ( b − ) -the Yang-Baxter relation (2.5) holds identically provided (2.9), (2.10) are true. An important case is when the generators x 0 , x 1 are realized as integrals over the vertex operators (2.11) The latter are required to satisfy the braiding relation along with the quasiperiodicity condition Here the operator Ω obeys and can be identified with As was pointed out in the work [16], using the braiding relations (2.12) it is possible to express monomials built from the generators x 0 and x 1 in terms of the ordered integrals This way, the formal power series (2.8) can be brought to the form The latter is recognized as the path ordered exponent The first realization of the generators x 0 , x 1 in terms of the vertex operators was proposed in refs. [14][15][16]. It reads as (2.19) where φ = φ(u) stands for the chiral Bose field satisfying the Operator Product Expansion (OPE) The parameter β is related to q, entering into the braiding relation (2.12), as q = e iπβ 2 . It is taken to lie in the domain If the field φ is assumed to be quasiperiodic it can be expanded in a Fourier series of the form 3 In turn the vertices satisfy the quasiperiodicity condition (2.13) with Ω = e 4πiβâ 0 . It is easy to see that the operators τ (λ) (2.6) commute with the zero mode momenta, and hence act invariantly in the Fock space F P . The latter is generated by the action of the creation operatorsâ n (n = −1, −2, . . .) on the Fock vacuum | P :â n | P = 0 (n = 1, 2, . . .) ,â 0 | P = P | P . (2.25) 3 Parafermionic realization of U q ( b − ) In ref. [18] a generalization of (2.19) was proposed, which goes along the following lines. Suppose we are given the algebra sl k (2) with central charge k = 1, 2, . . . . The Kac-Moody currents obey the OPEs

JHEP09(2021)201
As is well known [22], they admit a realization in terms of the Z k parafermionic fields ψ ± and the chiral Bose field φ (2.20): The parafermionic fields are quasiperiodic withΩ k being the operator of the Z k charge: where again β ∈ (0, 1). Note that the vertex operators can be interpreted as a oneparameter deformation of the currents J ± from the case β = 1. The braiding relations (2.12) and quasiperiodicity condition (2.13) are satisfied with Let us now briefly discuss the diagonalization problem for the transfer matrices τ (λ) (2.6) with the vertex operators as in (3.5). To this end some basic facts concerning representations of the algebra of parafermionic currents are needed. The theory was developed by Fateev and Zamolodchikov in ref. [22], in the construction of the CFTs describing the multicritical points of the Z k statistical systems [23] (certain generalizations of the Z 2 invariant Ising model). The algebra contains a set of non-local currents {ψ n } k−1 n=1 with conformal dimensions Their defining OPEs are invariant w.r.t. Z k transformations ψ m → ω 2ma ψ m with a = 0, 1, . . . , k − 1 and ω = e − 2πi k . All the parafermionic currents can be generated through the OPE of the currents with the lowest conformal dimensions which carry the Z k -charges +2 and −2, respectively (see eq. (3.4)).
The OPE of the fundamental parafermions is of special interest. It has the form

JHEP09(2021)201
All the fields occuring in the expansion are Z k neutral. The field W 2 , having Lorentz spin 2, generates the Virasoro algebra with central charge Further terms in the OPE (3.9) involve a set of local fields W 3 , W 4 , . . . , labeled by their value of the Lorentz spin. They form the WA k−1 algebra introduced in refs. [24,25] with the special value of the central charge c = c k (3.10). The chiral component of the Hilbert space of the Z k CFT can be decomposed into irreps V (k) j of the chiral algebra. Here the subscript j stands for the highest weight of the irrep with highest weight vector σ (k) j having conformal dimension (3.12) It will be further assumed that and | σ (k) k 2 −j possess the same non-vanishing conformal dimension. However they are distinguished by their value of the Z k charge. In fact the whole space V j is naturally splitted on the invariant subspaces of the operatorΩ k [22]: (3.14) The lowest possible conformal dimension in the subspace V (k) j,m is given by possess a structure of the highest weight irrep of the WA k−1 algebra. From the formal algebraic point of view such an irrep is obtained by factorizing a highest weight module over submodules of "null-vectors". In the case under consideration the highest weight can be thought of as a pair of eigenvalues (∆, w) of the mutually commuting operators corresponding to the highest weight vector. In the case of V (k) j,m with m = 2j, 2j − 2, . . . , −2j the conformal dimension ∆ is given by (3.15), while w = w j,m with

JHEP09(2021)201
For m = 2j + 2, . . . , 2k − 2j − 2, the highest weight is ∆ k 2 −j,k−m , w k 2 −j,k−m . The primary state (the highest weight vector) w.r.t. the WA k−1 algebra will be denoted by σ Let us return to the transfer-matrices τ (λ) built from the vertex operators First of all it is straightforward to show that they commute withâ 0 andΩ k : This implies that τ (λ) acts invariantly in the space V where The eigenvalues of this operator are given by with L a non-negative integer. As such I 1 naturally introduces a grading in V (k) j,m ⊗ F P . The dimensions of each level eigenspaces characterized by given L = 0, 1, . . . is finite. Due to the commutativity condition (3.19) the transfer-matrix acts invariantly in each level subspace of V (k) j,m ⊗ F P . An immediate consequence is that the primary state σ (k) j,m ⊗ |P is an eigenvector of the transfer-matrix.
One can address the diagonalization problem of τ (λ) in the level subspaces of V (k) j,m ⊗F P . It turns out that the transfer-matrices are part of a larger commuting family. The latter includes the operators a ± (λ) ∈ End V (k) j,m ⊗ F P : τ (λ), a ± (λ) = a + (λ), a − (λ) = 0 , (3.22) which satisfy the Baxter type relation The definition of a ± (λ) is similar to that given by eqs. (2.18), (2.6) and involves the vertex operators V ± (3.5) and Ω (3.6). However, the path-ordered exponent now contains the generators of the q-oscillator algebra E ± and H: (3.24) Let ρ ± be representations of this algebra such that the traces

Local IM
Among the operators which commute with the transfer-matrix, a special rôle belongs to the local Integrals of Motion (IM) [4]. These are a set of mutually commuting operators which can be written in the form with T s+1 being a chiral local density of integer Lorentz spin s + 1. Remarkably, for a given choice of vertex operators, there exists a purely algebraic procedure which, in principle, allows one to explicitly build the local IM. Later, a generalization of the vertex operators k will be proposed, which gives rise to new commuting families involving τ (λ), a ± (λ) as well as the corresponding sets of local IM. Here, for future references, we illustrate the construction of {I s } for the basic case. Let L (s) be the linear space of chiral local fields of Lorentz spin s built out of the Bose field φ and the fundamental parafermions ψ ± . For given positive integer s it is a finite dimensional space. We choose one of the vertices, say V + , and consider the linear subspace W (s) ⊂ L (s) made up of local fields X s such that the singular part of the OPE X s (u)V + (v) is a total derivative in v: In the physical slang, one says that the fields X s commute with the "screening charge" Suppose X s and Y s both commute with Q. By construction any local field which appears in the OPE X s (u)Y s (v) also commutes with the screening charge. Hence the direct sum ⊕ s≥1 W (s) possesses the structure of an operator algebra. Starting from the seminal work of A.B. Zamolodchikov [24], algebras of such type are referred to as the W -algebras.

JHEP09(2021)201
In the case under consideration the space L (1) contains only the field ∂φ, while W (1) = ∅. The space L (2) is spanned by (∂φ) 2 , ∂ 2 φ and the W 2 current occuring in the OPE of the fundamental parafermionic fields ψ ± (3.9). It is easy to check that, up to an overall multiplicative factor, there is a single spin 2 field commuting with the screening charge, is proportional to the derivative ∂X 2 . The space W (4) contains ∂ 2 X 2 and (X 2 ) 2 . The latter is a spin 4 local field that is the first regular term in the OPE: It turns out that dim W (4) = 3, i.e., together with the descendants of X 2 there is an extra field X 4 . Its construction can be simplified if one takes advantage of the bosonization formulae for the parafermion currents [28][29][30] This involves two chiral Bose fields with α(u)α(v) = − 1 2 log(u−v)+O(1) and similarly for γ. Then the field X 4 is a certain differential polynomial built from ∂φ, ∂α, ∂γ. In principal one can proceed further and explicitly describe the higher spin components W (s) . The original definition of the parafermion algebra requires that k is a positive integer. Nevertheless it can be treated as an arbitrary complex number in the bosonization formulae (4.6). This makes it possible to introduce the W -algebra associated with the screening charge Q (4.3), for arbitrary k as was done in the unpublished work of V.A. Fateev as well as in refs. [26,27]. In the last paper ⊕ ∞ s=1 W (s) was referred to as the corner-brane W -algebra. The second vertex operator V − comes in to play when obtaining the local IM from the fields in the W -algebra. One searches for local fields T s+1 ∈ W (s+1) such that the OPE of T s+1 and V − possesses the following structure Note that, if T s+1 exists, it is defined up to the addition of a total derivative ∂X s (∀X s ∈ W (s) ). This ambiguity has no affect on I s , which is expressed through an integral as in (4.1). There is also an ambiguity in the overall multiplicative normalization of T s+1 , which is carried over to the local IM. Following the arguments of ref. [16], it is expected that I s commutes with the transfer-matrices and a ± (λ): τ (λ), I s = a ± (λ), I s = 0 .

JHEP09(2021)201
In the case under consideration the operator I s exists and is unique (up to overall normalization) only for odd s = 1, 3, 5, . . . . For example, for s = 1 the density T 2 = X 2 from (4.4), while the explicit formula for I 3 was originally obtained in refs. [31,32]. Further facts on this commuting family of local IM can be found in [34].

Prototype example of the ODE/IQFT correspondence
The most effective way for computing the spectrum of the operators τ (λ), a ± (λ) and j,m ⊗ F P is provided by the ODE/IQFT correspondence. As part of our study of the new commuting family, a class of ODEs for the eigenstates will be proposed. In order to prepare for that discussion it would be useful to demonstrate the approach on the eigenvalues of a + (λ) corresponding to the primary state | σ (k) j,m ⊗ P . Let us start with the simplest case when j = m = 0. According to the work [18], one should consider the linear differential equation Here ξ > 0 while k is assumed to be a positive integer. For m(A) ≥ 0 one can introduce the Jost solution, which asymptotically approaches to a plane wave It turns out that as a function of A it is meromorphic. This allows one to unambiguously define Θ (←) A (x) for any complex A except for a discrete, pure imaginary set of values, where it develops simple poles. The function Θ (←) −A (x) is another solution to the ODE (5.1), which is linearly independent from Θ (←) A (x) when A = 0. Since the potential in (5.1) grows rapidly for large positive x, the ODE admits a solution which decays at x → +∞. We denote it by Θ (→) and specify its overall normalization through the asymptotic condition Θ (→) e −F (y)+o (1) as y = x + log(µ) → +∞ , where F (y) = 1 4 (1 + ξ)k y ( 2 F 1 stands for the Gauss hypergeometric function) and The solution Θ (→) may be introduced for any complex value of A 2 and is an entire function of this parameter. Consider the connection coefficients for the ODE (5.1), which are expressed via the Wronskian, (5.6)

JHEP09(2021)201
Here the overall factor has been chosen to ensure the existence of the limit The connection coefficient W is a meromorphic function of A. For given A it turns out to be an entire function of µ and, therefore, admits a convergent Taylor series expansion. This is related to the fact that in the variable y = x + log(µ), the ODE takes the form The term δU = (e (1+ξ)y + µ e ξy ) k − e (1+ξ)ky analytically depends on µ and can be treated perturbatively so long as ξ > 0, k = 1, 2, . . . . In ref. [18] it was shown that the ratio W (µ)/W (0) coincides with the eigenvalue of a + (λ) corresponding to the primary state j,m ⊗ |P with j = m = 0, provided the following identification of the parameters is made: and also Recall that the vertices V ± are non-local fields with fractional conformal dimensions. In writing the µ − λ relation (5.9) we assume that the phase ambiguity of V ± is specified through the following condition imposed on the OPE We now turn to the eigenvalues corresponding to the states σ (k) j,m ⊗ |P with j and m not necessarily zero. To the best of our knowledge the ODE/IQFT correspondence for this case has not yet been discussed in the literature despite that it follows immediately from the results of the works [27,34]. The generalization of the ODE (5.1) reads as Together with ξ and k this equation involves three extra parameters. With the specialization B 2 = A 2 and C 2 = 1 4 , it boils down to (5.1). As with the previous case one can consider the Wronskian W (µ) (5.6). However an important difference now is that W (µ) turns out to be an entire function of log(µ) rather than µ. This is related to the fact that when the ODE is rewritten in the form (5.8), the term δU is no longer an entire function of µ.
For the analysis of eq. (5.12) it is convenient to perform the change of variables

JHEP09(2021)201
bringing the ODE to the form For generic values of the parameters it possesses three singular points at z = 0, −1, ∞. Let us explore the condition when the singularity at z = −1 becomes apparent. First we consider the case κ = 0 when (5.14) belongs to the Fuchsian class. The test for the apparent singularity is well known and can be found, for example, in ref. [35]. Suppose we are given the equation and V (z) admits a Laurent expansion of the form The singularity at z = w is apparent if and only if v 0 = j(j + 1) with j = 0, 1 2 , 1, . . . (5.17) and Here the polynomials F j are defined through the determinant In the simplest cases (5.18) reads explicitly as Applying the above condition to (5.14) with κ = 0 one concludes that the singularity at z = −1 is apparent provided the parameters A, B and C satisfy the conditions

JHEP09(2021)201
Suppose now that κ = 0. For k = 1, 2, 3, . . . (5.22) the term ∝ κ has a zero of order k at z = −1. Hence its presence will not affect the conditions (5.17) and (5.18) as long as When z = −1 is an apparent singularity we define the normalized connection coefficients in the frame (z, Ψ) (5.13) as Then W j,m,A admits a power series expansion in µ. 4 It is possible to show that It also confirms eq. (5.9). This way one identifies D j,m,A (µ) with the eigenvalue of a + (λ) corresponding to the state σ Regarding the other members of the commuting family, the eigenvalue of a − (λ) for the state σ (k) j,m ⊗ |P coincides with D j,−m,−A (µ). The eigenvalues of the transfer-matrices τ (λ) are also given by certain connection coefficients for the ODE (5.14). Those of the local IM can be extracted from the large µ behaviour of D j,m,A within the standard WKB technique. All this is well known and widely discussed in the literature and will not be elaborated on here.

JHEP09(2021)201
where φ a and ψ (a) ± are r independent copies of the chiral Bose and Z ka parafermionic fields, respectively. Also the positive integer K stands for The definition (6.1) has been arranged so that the braiding relations ± and φ a , respectively, while are mutually local fields. Notice that the field ϕ satisfies the OPEs We will assume that the non-local fields V which is similar to (5.11). Let V ja be the space of representation of the parafermionic algebra generated by the fundamental parafermions ψ Pa is the space of representation for the Heisenberg algebra generated of the Fourier modes of the field ∂φ a with P a being the value of the corresponding zero mode momentum. Introduce the space where and s is an arbitrary number. The operators V and hence invariantly within the direct sum In other words, for any given P 0 and j 1 , . . . , j r , where we use the shortcut notation Each of the components in the linear decomposition (6.13) is an eigenspace for the zero mode momentum of the field ϕ (6.5): with the corresponding eigenvalue An important property is that the non-local fields V (a) ± ∈ End H j obey the quasiperiodicity condition

JHEP09(2021)201
with the operator valued factor being independent of a = 1, 2, . . . , r. The relations (6.3) and (6.19) imply that the vertex operators − , (6.20) satisfy the conditions (2.12)-(2.14) with The parameters {z a } r a=1 entering into the definition of V − may be arbitrary. The operators are expected to obey the quantum Serre relations (2.9). Together with they provide a realization of the generators of the Borel subalgebra U q ( b − ). 5 In all likelihood, it can be understood as coming from r consecutive applications of the comultiplication to the operators x 0 , x 1 with r = 1.
As soon as a representation of the Borel subalgebra U q ( b − ) is specified, an infinite family of commuting operators can be constructed via the general definitions (2.6) for the transfer-matrices and (3.26) for a ± (λ). They commute with the zero-mode momentum of the field ϕ and the operatorÛ (6.18): As a result, they act invariantly in the subspaces H j (6.13) with fixed value of these operators. The latter will be denoted by H j,m,P ⊂ H j : Here the integer m is defined modulo 2K. We employ the conventions that

JHEP09(2021)201
As follows from (6.17), the value of the zero-mode momentum P must obey the relation Explicitly, the decomposition of H j,m,P is given by where each of the components in the direct sum is described by formula (6.10). The summation is taken over the set In view of eqs. (6.17), (6.18) and (6.28), H j,m,P is a common eigenspace for the operatorŝ a 0 andÛ . In turn Similar to the case r = 1 it is expected that the commuting family of operators contains an infinite set of local IM for any positive integer r. As was already discussed, their construction starts by exploring the W -algebra built out of the local fields commuting with the screening charge du V + . It is easy to check that the algebra includes the field where the spin 2 local fields W appear in the first non-trivial term in the OPE ψ , see eq. (3.9). The field T 2 forms a closed subalgebra, similar to (4.5), with the central charge The study of this W -algebra for r > 1 is a considerably more difficult task than that for the corner-brane W -algebra corresponding to r = 1. Some results concerning the case when all k a = 1 are presented in sections 9 and 10. Having considered in detail various particular cases, we conclude that there exists a non-trivial W -algebra for any r = 1, 2, . . . . It will be denoted by W central charge of the Virasoro subalgebra. 6 We conjecture that, out of the fields from the W -algebra, one can construct a set of local IM such that where T 2 is the local density defined in eq. (6.33).
The space H j,m,P is a highest weight module of W (c,r) k . As usual, a grading is introduced by the operator I 1 . Its spectrum is bounded from below, while the eigenvalues are given by I The representations H j,m,P are labeled by the set j = (j 1 , . . . , j r ) with j a = 0, 1 2 , 1, . . . , 1 2 k a and the integer m ∼ m + 2K, which we take to lie in the interval −2J, −2J + 2, . . . 2K − 2J − 2. The latter can be further subdivided into two subsets Notice that the transformation j a →ǰ a and m →m witȟ maps m from the interval (ii) to the interval (i) with J replaced byJ = r a=1ǰ a . In the next section it will be proposed that the spaces H j,m,P and Hˇj ,m,P should be treated as equivalent representations of the algebra W (c,r) k : H j,m,P ∼ = Hˇj ,m,P . (6.39) We'll mainly focus on case (i), and comment on case (ii) as required.

ODE for the W -primary states
Provided that the integer m is restricted as in (6.37a), the spectrum of the operator I 1 in H j,m,P is given by Similar as for the corner-brane W -algebra, W (c,r) k can be defined for arbitrary values of ka if one uses the bosonization formulae (4.6) for each copy of the parafermionic currents ψ

JHEP09(2021)201
The dimensions of the subspace of primary states, H (0) j,m,P , can be read off from the formula . When m = ±2 r a=1 j a the corresponding eigenspaces are one-dimensional and spanned by the ground states The path-ordered integral formulae for τ (λ) and a ± (λ) allow one to express their matrix elements as a convergent power series in λ 2 . The corresponding expansion coefficients are given by the well defined multifold contour integrals (for further details see [16]). In the simplest set-ups, it is possible to calculate them directly. For example for the vacuum states (7.3), one can show that a + (λ) e j,±2J,P = a (vac,±) + (λ) e j,±2J,P : log a (vac,±) 2j a z a (7.5) and However, the computation of the higher order coefficients H 2 , H 3 , . . ., even for the ground states, turns out to be highly cumbersome. A practical way of studying the spectrum of τ (λ) and a ± (λ) is through the ODE/IQFT correspondence for the commuting family, which is proposed below. Generalizing eq. (5.14) we consider the ODE where ξ > 0 and k = (k 1 , . . . k r ) is a set of positive integers. As before the singularities at z = z a (a = 1, . . . , r) are required to be apparent. Repeating the same line of arguments

JHEP09(2021)201
as for the case r = 1 one arrives at the following conditions imposed on the parameters. An immediate one is that Then for some given (half-)integers (j 1 , . . . , j r ) and fixed value of A, the set γ = (γ 1 , . . . , γ r ) should solve the system of algebraic equations Here v in the vicinity of z = z a : The polynomials F j (v 1 , . . . , v 2j+1 ) are given by the determinant (5.19).
Numerical work shows that, if A is generic, the algebraic system (7.9) possesses r a=1 (2j a + 1) solutions. These are splitted on the classes labeled by an integer m such that where m is restricted as in eq. (6.37a), i.e., m = −2J, −2J−2, . . . , 2J−2, 2J. The connection coefficients are defined similar to the case r = 1. Namely, They are treated as a function of Let γ be a solution of (7.9). We expect that there exists a state labeled by this set such that it is a simultaneous eigenstate for the operators of the commuting family with a ± (λ) e j,m,P (γ) = D j,±m,±A (µ|γ) e j,m,P (γ) .
The parameters on the ODE and the field theory side are related as The states e j,m,P (γ) form a basis in H (0) j,m,P . The support for the proposed ODE/IQFT correspondence is based on the arguments that were originally developed in the works [10,11]. Following ref. [10], it is straightforward to derive a set of functional relations for various connection coefficients. The most fundamental of these is the so-called quantum Wronskian relation It holds true as long as all the singularities at z = z a are apparent. Furthermore in this case the connection coefficients turn out to be entire functions of µ.
On the other hand, all the eigenvalues of a ± (λ) obey the quantum Wronskian relation, which comes from the operator valued relation (3.28) with = 0. This coincides with (7.19) upon making the identification (7.18) as well as ξ = β 2 1−β 2 and µ ∝ λ 2 . The precise µ − λ relation (7.17) can be established via a first order perturbative calculation in µ of D j,m,A for (7.7) with all j a = m = γ a = 0 and comparing the result with the corresponding specialization of eq. (7.5).
For given A and j = (j 1 , . . . , j r ), there exists two solutions of the algebraic system (7.9) which are particularly simple: where The ODE (7.7) in this case takes the form Here, together with the notations (7.21), we use The integer m (7.12) for the solutions γ (vac,±) a coincides with ±2J. Numerical work shows that they correspond to the ground states e j,±2J,P (7.3). Note that, in view of eq. (7.18), (7.24)

JHEP09(2021)201
It is important to keep in mind that the connection coefficients D j,±m,±A are defined with m ranging from −2J to +2J as in eq. (6.37a). In order to cover the case (6.37b), we introduceǰ a ,m andǍ via the formulaě and consider the ODE Here the set {γ a } solves the system of equations (7.9) with A and j a replaced by their "checked" counterparts. In turn, eq. (7.12) is swapped for Then the ODE/IQFT correspondence for the primary states is formulated as a ± (λ) e j,m,P (γ) = Dǰ ,±m,±Ǎ (µ|γ) e j,m,P (γ) .
and therefore the quantum Wronskian relation for the connection coefficients Dǰ ,±m,±Ǎ is identical to (7.19) as long as the parameter P is kept fixed. This is the first piece of support for the isomorphism (6.39).

Relation to the Bethe ansatz equations for the sl(2) Gaudin model
There is an alternative description of the sets γ, which solve the algebraic equations (7.9) going back to the works [2,[5][6][7]. Assuming that t 0 (7.10) is given, consider the Riccati equation for the unknown function f 0 = f 0 (z), so that t 0 is the Miura transform of f 0 . One can search for a solution using the ansatz  (7.34) and The requirement that r m = 0 yields an algebraic system for the auxiliary parameters x (+) m : Note that for M + = 0, formula (7.34) gives back the vacuum solution γ (vac,+) a (7.20). Moreover it is straightforward, using eqs. (7.36) and (7.34), to calculate the sum r a=1 z a γ a . This yields (7.12) with m being related to the non-negative integer M + as Having at hand a solution of eq. (7.36), the set γ obtained via (7.34) would automatically obey the conditions (7.9), which guarantee that the singularities at z = z a of the ODE are apparent. Our numerical work suggests that for generic A there is a one-to-one correspondence between x (+) solving (7.36) and γ. However, for some specific values of A the correspondence breaks down. For example upon taking A = − i 2 , eqs. (7.34), (7.36) imply that r a=1 γ a = 0. On the other hand, for A 2 = − 1 4 , the system (7.9) admits solutions with r a=1 γ a = 0. The ODE depends on A 2 rather than A. However eqs. (7.34) and (7.36) are not invariant if the sign of A is flipped. For this reason the set γ solving (7.9) may be alternatively expressed as where S (0) is a spin operator that acts in a highest weight infinite dimensional representation (Verma module) of sl (2). The corresponding Casimir would be given by The diagonalization problem can be considered in the finite dimensional eigenspaces of the operator S 3 , which commutes with the Hamiltonians. If instead of the highest weight, one takes the lowest weight representation for S (0) with lowest weight −1 − j 0 , the spectral problem for the Hamiltonians (7.41) would lead to the Bethe ansatz equations (7.39).

Excited states ODE
In ref. [34], developing the ideas from [11], the ODE/IQFT correspondence was proposed for the highest state irreps of the pillow brane W -algebra. Here, following the construction from that work, we extend the correspondence to all the eigenstates of the commuting family of operators. These form a basis in the highest weight irrep of W (c,r) k . Such basic states e j,m,P (γ; w) ∈ H (L) j,m,P (7.43) would be distinguished by the sets (γ; w) = (γ 1 , . . . , γ r , w 1 , . . . , w L ), which are solutions of a certain algebraic system. The latter has already been discussed for L = 0. The system of equations imposed on (γ; w) for general L = 0, 1, 2, . . . is described as follows.
The corresponding residues are given by the sets γ = (γ 1 , . . . γ r ) and Γ = (Γ 1 , . . . , Γ L ). through the Laurent expansion of t L (z) in the vicinity of z = z a and z = w α , respectively: Assuming that (z 1 , . . . , z r ) and (w 1 , . . . , w L ) are given, the residues γ and Γ are determined through the solution of the coupled system of r + L equations: Here the polynomials F j (v 1 , . . . , v 2j+1 ) are defined via the determinant (5.19). The meaning of these equations should be obvious at this point: if all 2j a + 1 are positive integers they form the full set of conditions that all the singularities of the Fuchsian differential equation where t L (z) is given by (7.44), while We impose that all the singularities except z = 0, ∞ for any value of the parameter κ are apparent. Eqs. (7.47) guarantee this property for κ = 0 so that they are necessary conditions. Furthermore, as was already discussed, for generic κ the admissible values of j a must be restricted as in (7.8): Also it turns out that the positions of the apparent singularities w α may not be chosen at will. Instead they should satisfy L extra conditions (for details see [34]) or explicitly This expresses the set of residues Γ through w. Thus (7.47) becomes a system of r + L equations imposed on the r + L variables (γ 1 , . . . , γ r ; w 1 , . . . , w L ). Similar to the case of L = 0, its solutions are splitted on the classes labeled by an integer M such that j a (j a + 1) .

JHEP09(2021)201
This integer takes the values with some M max . Numerical work suggests that An alternative description is provided if one considers, rather than t L (z), the solution of the Riccati equation t L = f 2 L −∂ z f L . As explained in section 7.2, it leads one to introduce the auxiliary sets . These would parameterize the residues γ a as together with (7.52) forms a closed system of L + M ± equations, which determine the sets x (±) = (x (±) 1 , . . . , x (±) M + ) and w = (w 1 , . . . , w L ). Finally γ = (γ 1 , . . . γ r ) is obtained via eq. (7.56).
In order to formulate a precise conjecture regarding the ODE/IQFT correspondence, there is an important issue that needs to be addressed. According to eqs. (7.53)-(7.55), the algebraic system (7.47) admits solutions (γ; w) with |M| > 2J. However, there exists a map (γ; w) → (γ,w) such that the transformed set satisfies the equations similar to (7.47) with the parameters j a , L and A replaced bỹ

JHEP09(2021)201
For M > 2J the transformation is described as follows. Suppose that w α and γ a are given. The setw α is the solution of the equations: which is uniquely specified by the extra conditions Oncew is found, the residuesγ are obtained from If one now considers the change of variables For the other case with M < −2J, the sign of A andÃ in eqs. (7.62)-(7.66) needs to be flipped. An important point is that since Ψ and Ψ are related through the action of a first order differential operator, the monodromic properties of the ODEs (7.48) and (7.67) are the same. In consequence, the corresponding connections coefficients for these differential equations coincide. Thus, although there exists the solutions sets (γ; w) with |M| > 2J, by performing the above transformation, the value of |M| can be reduced. By successive applications if necessary, we expect that the integer M may always be brought to the interval M = −2J, −2J + 2, . . . , 2J .

JHEP09(2021)201
The connection coefficients D j,±m,±A (µ | γ, w) are introduced similarly as in eq. (7.13). They turn out to be entire functions of for any given choice of (γ, w). Since the singularities w α that were introduced are apparent, the quantum Wronskian relation (7.19) is still satisfied. Interpreting D j,±m,±A (µ | γ, w) as an eigenvalue of the operator a ± (λ) for a certain state in the level subspace H (L) j,m,P , the zero-mode momentum P would be related to A as in eq. (7.18). With this set-up, we conjecture that the connection coefficients coincide with certain eigenvalues of a ± (λ): j,m,P . Of course, the ξ − β and µ − λ relations remain unchanged and are given by (7.17).
For the second case in (6.37) we make the identification The conjectured ODE/IQFT correspondence is described as Hereǰ,m andǍ are defined in (7.25). It is important to keep in mind that (γ,w), which label the states in H (L) j,m,P , are solutions of the algebraic system (7.47) and satisfy (7.53) with A, j a replaced by their "checked" counterparts and M =m.
It is expected that the eigenstates of a + (λ) form a basis in the representation H j,m,P , which is fully specified by its eigenvalues. The conjecture (7.73) implies that the spectrum of a + (λ) in the space H j,m,P with m = 2J, . . . , 2K − 2J − 2 coincides with its spectrum in Hǰ ,m,P . This leads us to propose the isomorphism of these two representations of the W (c,r) k algebra, see eq. (6.39).

Some comments concerning the literature
Isotropic limit. In ref. [14] an expression for the vacuum eigenvalue of the transfermatrix τ 1 2 (λ) was presented for the case r = 1. It coincides with the grand canonical partition function of a 1D plasma of alternating charges. The latter, since the seminal work of Anderson and Yuval [39], is well known to be the partition function of the one-channel anisotropic Kondo model, where the impurity has spin 1 2 [40]. This allowed the techniques of ref. [14] and its further developments, including the ODE/IQFT correspondence, to be transferred to different variants of the quantum impurity problem such as the quantum dot [41] and the Coqblin-Schrieffer model [42]. In the work [18], the partition function of the multichannel anisotropic Kondo model was expressed in terms of the solution of the ODE (5.1). Renewed interest in the Kondo problem came from refs. [43,44], where some

JHEP09(2021)201
previously obtained facts concerning the isotropic case were rediscovered. Here we explain the relation between the differential equation (7.7) and the one that was considered in those works.
The isotropic limit corresponds to taking β → 1 − . It is a complicated problem to perform the limit starting with the path-ordered exponent (2.18), which in the context of the impurity problem is interpreted as an imaginary time evolution operator in the interaction picture. One of the main advantages of the ODE/IQFT correspondence is that it is well adapted for exploring such a limit. The corresponding technique was developed in the works [41,42] and can be straightforwardly applied to the ODE (7.48).
Let us make a uniform shift of the variable z and the positions of the apparent singularities and rewrite the ODE (7.48) using the parameters Of course such an innocent change of variables and redefinition of the parameters should not affect the properties of the ODE. We set ε = 2 ξK (7.76) and consider the limit ξ → +∞ keeping κ and A fixed. This corresponds to β = ξ 1+ξ → 1 − . Taking the isotropic limit in the differential equation (7.48), one obtains (7.78) and the limiting form of (7.52) reads as The sets γ and w still satisfy the algebraic system (7.47). Of course, now v (z). In the isotropic limit the relation (7.18) becomes For the "vacuum" case corresponding to the states e j,±2J,P (7.3), one should set L = 0 and choose γ a to be The ODE (7.77) boils down to This equation with r = k = 1 and j 1 = 0 originally appeared in refs. [41,42] in the description of the isotropic Kondo model. It is worth mentioning that the parameter 2iP can be identified with H/T , where H is an external local magnetic field acting on the impurity spin, while T stands for the temperature. Also κ = T K /T with T K being the Kondo temperature.
The ODE (7.83) with r > 1 and arbitrary k a was introduced in the recent work [43]. Note that in the case A (vac,±) = 0 it corresponds to the chiral primary states for the ⊗ r a=1 sl ka (2) WZW model. Indeed, as follows from (7.1) with P = m j a , the value of the local IM I 1 is given by Gaudin limit. The ordinary differential equation for the isotropic case (7.77)-(7.78) still admits an interesting limit, which can be interpreted as a double scaling limit of the original system. 7 One performs the change of variables z = δ y (7.85) and takes δ → 0, keeping fixed the combinations Further assuming that A = 0, i.e., P = m 2 √ K (7.87) the differential equation becomes and v α are solutions of an algebraic system, which ensures that the points y = y a , v α are apparent singularities of the ODE. Of course, this system can be obtained via a limit of (7.47). The differential equation (7.88) was originally proposed by Feigin and Frenkel [5] in their study of the affine Gaudin model.
8 Large µ asymptotic expansion of D j,m,A

Leading behaviour
The operators a ± (λ), which depend continuously on the spectral parameter λ, are generating functions of the family of commuting operators. In this regard, of special interest is the large λ expansion of a ± (λ) [15]. With the ODE/IQFT correspondence at hand, it can be explored via the study of the connection coefficients D j,±m,±A (µ) at large µ. Below we restrict our attention to D j,m,A .
A brief inspection of the ODE shows that the connection coefficients depend on µ and the parameters {z a } r a=1 only through the combinations z a µ. Thus a multiplication of all the z a by the same factor can be absorbed into a redefinition of the spectral parameter. Focusing on the case when z a = 0 for any a, without loss of generality, one can set The asymptotic coefficients are given by and R j,m,A = (−1) Here the parameters ε and δ should be taken to be sufficiently small, such that z a / ∈ D for all a = 1, . . . , r. The contour C lies inside D and is the image of a Pochhammer loop on the Riemann sphere under stereographic projection. The homotopy class of the loop is schematically shown in figure 1. It winds around the south and north poles of the sphere, which are mapped to z = 0 and z = ∞, respectively.
The asymptotic coefficient R j,m,A (8.4) is interpreted as an eigenvalue of the so-called reflection operator. 8 For the case r = 1 and k = 1 it was discussed in refs. [45][46][47][48]. 9 It also appears in various applications of the ODE/IQFT correspondence.

Eigenvalues of the local IM in the GAGM
The omitted terms in (8.2), denoted by o(1), encode the eigenvalues of the local and socalled dual non-local integrals of motion. It is expected that D j,m,A (µ) admits a systematic asymptotic expansion of the form where the eigenvalues of the local and dual non-local IM appear in the coefficients of the formal power series log(B) = o(1) and log(X) = o(1), respectively. The series log(B) can be studied within the standard WKB technique. Carrying over the analysis in refs. [33,34], one finds The asymptotic formula (8.2), when combined with the quantum Wronskian relation (7.19), yields an interesting identity 9 In ref. [47], the auxiliary variables x (±) m are not employed and the eigenvalues of the reflection operator are expressed solely in terms of the apparent singularities wα. This results in more cumbersome formulae.

JHEP09(2021)201
There is a simple recursion procedure to generate the densities U 2n for any n. For instance, Also in writing the numerical coefficients in the sum (8.8), it was assumed that the densities are normalized by the condition where the omitted terms contain derivatives and lower powers of t L (7.44). Let us focus on the first nontrivial integral q 1 . It turns out that the particular choice of the integration contour C is not essential. For this reason we consider where C is an arbitrary closed contour. Independently on the choice of C, the following relation holds true The coefficients in the sum read explicitly as Here we use the notations while In spite of the somewhat cumbersome expressions, the proof of the identity (8.13) is elementary. Namely, taking into account eq. (7.51), it is straightforward to show that Integrating both sides of (8.17) over the closed contour yields (8.13).

JHEP09(2021)201
In the context of the ODE/IQFT correspondence, q (C) 1 (8.12) should be interpreted as the eigenvalue of a certain local IM. Formula (8.13) shows that this quantity is a linear combination of I (a) 1 with a = 1, 2, . . . r, where the coefficients encode all of the dependence on the integration contour. Thus I Here a = 1, . . . , r but the summation index "b" in the first sum runs from 0 to r and Also d 0 , defined by eq. (8.16), can be rewritten as It is instructive to compare (8.22) with the expression (1.4) for the energies E a in the original Gaudin model. In principle, similar computations can be performed for the higher local IM. Upon obtaining the expression for the density U 2n , one should establish the identity which is the analogue of (8.17). The eigenvalues of the local IM would be read off from the expansion coefficients in the sum. The calculations are purely algebraic and straightforward, but rather cumbersome. For the integral I (a) 3 , they were performed only for the case r = 1 in ref. [34].

Some comments on the eigenvalues of the dual nonlocal IM
The factor X in the asymptotic expansion (8.7) is a formal power series of the form The coefficients X n are interpreted as the eigenvalues of certain integrals of motion. However, contrary to I (a) 2n−1 the operators can not be represented as integrals over the local chiral densities. Following the terminology of ref. [15] we refer to these operators as the "dual nonlocal IM".
The appearance of X(µ) in the large µ expansion of the connection coefficient is related to the fact that the applicability of the WKB approximation to the ODE (7.44)-(7.52) breaks down in the vicinity of z = 0. This makes it difficult to provide a mathematically rigorous justification of the asymptotic series (8.7) and, in turn, to calculate the expansion coefficients X n . Nevertheless, the presence of X(µ) in (8.7) can be argued for similarly as it was done in ref. [10] for the case r = k = 1.
Consider the simplest variant of the ODE (7.48), (7.44): Taking into account the convention (8.1), the change of variables z = ζ −1 , Ψ(z) = ζ −1Ψ (ζ) brings the differential equation to the form The latter looks similar to the original ODE with the parameters substituted as ξ → −1 − ξ, z a → z −1 a . In the parametric domain −1 < ξ < 0, this formal substitution can be interpreted as a duality transformation, which allows one to relate the properties of the solutions in the vicinity of z = ∞ and z = 0.
The subject of our current interest is the domain with ξ > 0. It is mapped to ξ < −1 under the reflection ξ → −1 − ξ, which thus can not be interpreted as a duality transformation. However, since it relates the basic solutions Θ (←) ±A , defined through their behaviour at the regular singular point z = 0, and Θ (→) , which is the fast decaying solution at the irregular singularity z = ∞, it leads to an important property for the connection coefficients. Namely, while D j,m,A (µ) possesses the convergent Taylor series expansion in µ, its large µ asymptotic involves the formal power series X(µ) in the "dual" parameter µ: In this work, we forgo a systematic study of the eigenvalues of the dual nonlocal IM. An illustration for a specific example with r = 2 and k 1 = k 2 = 1 is provided in section 9.3. Finally we note that the eigenvalues of the transfer-matrices τ (λ) can be expressed through D j,m,A (µ) by using the set of well known functional relations. Hence, provided that JHEP09(2021)201 the large µ asymptotic expansion of the connection coefficients is available in the whole complex µ plane, one can straightforwardly derive the large λ expansion of the eigenvalues of τ (λ). The latter are of special interest in the context of the quantum impurity problem, where it permits a systematic study of the infrared fixed points in the theory. 9 Example with r = 2 and k 1 = k 2 = 1 The local IM for the GAGM have not been sufficiently studied in the general setup. The case r = k = 1 has been explored in detail in the context of the quantum KdV theory [14][15][16], while some results are known for r = 1 and k > 1. The simplest case beyond the quantum KdV model is when r is a positive integer but with all k a = 1. Then the representation of U q ( b − ) discussed in section 6 does not involve the parafermionic fields, leading to significant simplifications. This section is devoted to an analysis of the case r = 2 and k 1 = k 2 = 1.

W -algebra and the first local IMs
Taking r = 2 and k 1 = k 2 = 1 in eqs. (6.20) and (6.1), the vertex operators V ± become The exponent involving the Bose field θ can be "fermionized" as where χ 1,2 are a pair of chiral Majorana fermion fields, satisfying the OPE Then formula (9.1) becomes The construction of the local IM follows the procedure outlined in section 4. First we note that the OPEs of V + with the spin 3 2 and spin 2 fields, obey the condition (4.2) provided the parameter ρ is taken to be

JHEP09(2021)201
The fields S and G form the N = 1 super Virasoro algebra, whose commutation relations are encoded in the singular part of the operator product expansions with the central charge c sV = 3 2 − 6ρ 2 . The Majorana fermion χ 2 also obeys the condition (4.2) simply because the OPE χ 2 (u)V + (v) is not singular at u = v. The space of local fields satisfying (4.2) can be generated by the spin 3 2 field S and spin 1 2 field χ 2 . Among such fields we will focus on those having integer Lorentz spin, which form a closed operator algebra. The latter is a W -algebra that will be denoted by W (1,1) . In section 4, the linear subspace of chiral local fields of Lorentz spin s commuting with the screening charge du V + was denoted by W (s) . In the case at hand W (2) is a linear span of three fields G, χ 2 S and χ 2 ∂χ 2 . Out of these satisfy the second requirement (4.7) involving the vertex operator V − . There is a shortcut way of checking the last statement. In the construction of the densities (9.9), one can alternatively start by considering the space of local fields satisfying the condition (4.2), where V + is replaced by V − . Parameterizing z 1,2 via the complex numbers s and α as z 1 = is e −iα , z 2 = −is e iα (9.11) (if one adopts the convention (8.1), then s 2 = 1), the vertices take the form Then such local fields would be generated by the 3 2 and 1 2 spin fields S = ρ∂ + i∂ϕ χ 2 ,χ 1 = cos(α) χ 1 − sin(α) χ 2 . (9.14) Again we introduce the W -algebra, W The characteristic property of the densities entering into the local IM is that they are invariant under the reflection, up to a total derivative: Note that this immediately implies that the local IM commute with the reflection operator: 10 It is easy to check that the fields (9.9) possess the property (9.16). Below we'll argue that the local IM are linearly expressed in terms of I 1 and I 1 , whose eigenvalues are given by eq. (8.14) with r = 2, k 1 = k 2 = 1. Namely, In fact the first relation follows immediately. It was already mentioned (see eq. (8.21) and below) that its r.h.s. coincides with the local IM I 1 defined by (6.36) with the local density T 2 in (6.33). The latter, specialized to the case r = 2, k 1 = k 2 = 1, should be compared with the integrand in the first line of eq. (9.18), taking into account that One can proceed further and construct the higher spin densities T s+1 . It is easy to see that the linear subspace W (3) ⊂ W (c,2) 1 includes, together with the derivatives of the fields from W (2) , a single local field ∂χ 2 S. The latter does not satisfy the condition (9.16). The space W (4) , apart from total derivatives, contains five composite fields ∂SS, G 2 , χ 2 GS, ∂ 2 χ 2 S and G 2 χ 2 . The first two of them, together with the spin 7 2 field GS are defined through the regular terms in the OPE (9.8). The spin 4 local field G 2 χ 2 is the first regular term in the OPE G χ 2 (u)G χ 2 (v) as u → v. The explicit expressions in terms of the Bose Again, the eigenvalues of certain linear combinations of these operators are expected to coincide with the coefficients I

Irreps of the W -algebra
The space H j,m,P was introduced in section 6 and conjectured to be a highest weight representation of the algebra of local fields W (c,r) k . It is labeled by P , which is the value ofâ 0 -the zero mode momentum of ϕ (see eq. (6.26)), the set j = (j 1 , . . . , j r ) with j a = 0, 1 2 , . . . , 1 2 k a and the integer m ∼ m + 2K such that Moreover, the representations H j,m,P and Hǰ ,m,P witȟ k a are equivalent. In the case when r = 2, k 1 = k 2 = 1 we are left with four possibilities. The irrep H j,m,P corresponding to j 1 = j 2 = 1 2 and m = ±2 is the same as the one with j 1 = j 2 = m = 0. As a shortcut, we'll denote it by H 0,0 . Similarly the case with j = ( 1 2 , 0) and m = +1 is equivalent to the one where j = (0, 1 2 ) and m = −1. The corresponding space will be denoted by H 1,−1 . The notation H 1,+1 will stand for the irrep H j,m,P with j = ( 1 2 , 0) and m = −1, which is the same as j = (0, 1 2 ) and m = +1. Finally H j,m,P with j 1 = j 2 = 1 2 and m = 0 will be referred to as H 2,0 . Speaking in general, the W -primary states are the ones which are annihilated by  a = 1, . . . , r). As follows from eq. (7.2) the representations H 0,0 , H 1,−1 and H 1,+1 all contain only one linearly independent state, which we denote by e 0,0 , e 1,−1 and e 1,+1 , respectively. As for H 2,0 , there are two primary states e (±) 2,0 . A summary of the notation, together with the corresponding eigenvalues of the local IM I   Neveu-Schwarz sector. In this sector the fermion fields are antiperiodic, so that their Fourier mode expansions are given by The NS fermionic module (Fock space) F NS is generated by the action of the "creation" operators f Since all the fields of the W -algebra W Here we use the notation The latter coincide with the Fourier coefficients of the current Ramond sector. In the Ramond sector the Majorana fermions are periodic fields, so that The Ramond vacuum states |± R are annihilated by f (a) n with n > 0 and form a two dimensional representation for the algebra of zero modes

JHEP09(2021)201
Let σ A (A = x, y, z) be the Pauli matrices, such that σ z |± R = ± |± R . One can set The Ramond fermionic module F R , generated by the action of the creation operators f commute with (−1) F , which allows one to identify In turn, for the primary states The local IM (9.18) in the Ramond sector, written in terms of the creation/annihilation operators, read as and With these expressions it is easy to see that 1 come from a comparison of the eigenvalues quoted in the above formula and eq. (9.30) with the last two columns of table 1. One needs to also take into account that τ = −ρ z 1 +z 2 z 1 −z 2 and ρ =  for the primary states of the W (c,2) 1 algebra. These were computed directly using the expression (9.21) for the local IM as well as the formulae (9.30), (9.31) for the states e 0,0 , e (±) 2,0 and (9.41) for e 1,±1 . Here we use the notation g = 4P 2 − ρ 2 + τ 2 .
completeness, we present in table 2 the eigenvalues of the operators I The identifications immediately yield formulae for the characters The characters are the generating functions for the dimensions of the level subspaces. In view of the ODE/IQFT correspondence discussed in section 7, these formulae provide highly nontrivial predictions concerning the number of solutions of the algebraic system (7.47)-(7.51) on the apparent singularities of the ODE for r = 2 and k 1 = k 2 = 1.

Dual nonlocal IM
In section 8.3, we briefly mentioned a certain "duality transformation" and traced its appearance at the level of the ODE. On the quantum field theory side, the simplest instance where it was originally observed was in the work of Schmid [49] on the instanton calculus in dissipative quantum mechanics. Around ten years later, the duality was used in the computation of the conductance in a fractional quantum Hall system [50] and also appeared in the study of the quantum KdV theory [15]. The latter is the special case of the generalized

JHEP09(2021)201
affine Gaudin model with r = k = 1. For r = 2 and k 1 = k 2 = 1 the duality is manifest in essentially the same way, as will be discussed below. 11 An examination of the explicit formulae for the local IM (9.18) and (9.21) shows that they remain unchanged under the transformation so long as one simultaneously swaps the parameters The defining property of the corresponding local densities T (a) 2n is that they satisfy the OPEs of the form with the vertices from eq. (9.5). The invariance of the local IM under the "duality" transformation immediately implies that the similar OPEs hold true with V ± replaced by the "dual" vertices which are obtained from the original ones via the substitutions (9.46), (9.47). Introduce the operators (9.50) These would satisfy the Serre relations (2.9), but with q replaced by A realization of the Borel subalgebra U q ( b − ) is provided by x 0 , x 1 , as well as h 0 , which is defined such thatqh withÛ from eq. (6.26) (U 2 = +1 and U 2 = −1 in the Neveu-Schwarz and Ramond sector, respectively). Consequently, one can construct the dual operator a ± ( λ) from the vertices V ± and Ω 1 2 via the formulae, which are only notationally different to (3.26) and (3.27). 11 Some hints of the appearance of this duality were observed in the study of the Z2 invariant inhomogeneous six-vertex model in ref. [51]. The latter corresponds to the case z1 = −z2 = i. The relation between the GAGM and the lattice system is discussed in section 12 of this work.

JHEP09(2021)201
Notice that the transformation (9.46) flips the sign of ϕ and, therefore, its zero mode momentumâ 0 . As a result, the operators a ± would be given by the trace over a representation ρ ± of the q -oscillator algebra, subject to the requirement A comparison with the similar condition (3.25) for the representation ρ ± , appearing in the construction of a ± , suggests that under the duality transformation ρ ± → ρ ∓ .
The operator a ± ( λ) possesses the formal power series expansion in the "dual" spectral parameter: Based on the fact that the residue in the OPE T (a) 2n (u) V ± (v) is a total derivative, it is possible to argue, following the lines of [16], that the operators H (±) n commute with the local IM. They are referred to as the dual nonlocal IM.
Similar as in the case of the quantum KdV, one expects the eigenvalues of the dual nonlocal IM to appear in the large µ asymptotics of the connection coefficients. For instance, the formal series X(µ), which enters into the asymptotic formula (8.7) for D j,+m,+A (µ), coincides with the eigenvalue of a − ( λ), provided that the expansion parameters µ and λ are properly related. This way, formula (8.7) may be promoted to the operator relation: Here R is the reflection operator (9.15), whose eigenvalues coincide with R j,m,A (8.4), while B(µ), with eigenvalue B(µ) (8.8), involves only the local integrals of motion. Recall that µ ∝ λ 2 according to eq. (7.17) which, specialized to the case at hand, reads as (9.57) The relation between λ and λ can be established in the following way. Introduce µ through the formula which is just the "dual" version of (9.57). On the other hand, the analysis of the ODE performed in section 8.3 suggests that µ = µ −β −2 (see (8.29)). Combining this with the above equations one finds (9.59)

W -currents of Lorentz spin 2 and 3
When all k a = 1, the vertex V + defined by eqs. (6.20) and (6.1) takes the form Here {φ a } r a=1 is a set of independent chiral Bose fields subject to the OPE In what follows O ab will stand for the exponential operators which are local fields of Lorentz spin 2.

JHEP09(2021)201
Consider the r fields as well as the 1 2 (r − 1) r fields One can check that they commute with the screening charge du V + (u). In other words, X a , Y ab belong to the space W (2) . We argue that any spin two W -current can be expressed as a linear combination of the fields X a and Y ab . This way they form a basis in W (2) . The last statement is supported by the observation that the OPE of X a and Y ab generates no spin 2 fields, which are linearly independent from the basic ones. For example, a straightforward calculation yields the following result and C ab are some constants whose explicit form is unessential. This shows that the spin 2 fields generated in the OPE X a (u)X b (v) are linearly expressed in terms of the basic ones. The less singular term in (10.8) involves the spin 3 local fields They also commute with du V + (u). The explicit formula for Z ab reads as

JHEP09(2021)201
The singular part of the OPE of X a and Y bc is given by All the fields appearing in the singular parts of this OPE are linearly expressed in terms of the basic spin 2 fields X a , Y ab as well as the spin 3 fields Z ab . Finally one should consider the OPE involving the Y -fields only. Its singular part possesses the following structure Notice that so that the number of independent components of Z bdc is equal to 1 2 r(r − 1)(r − 2). The above analysis suggests that X a , Y ab form a basis in W (2) , while any spin 3 field from W (3) may be represented as a linear combination of Z ab , Z abc as well as the derivatives ∂X a and ∂Y ab . One can continue the process of generating the higher Lorentz JHEP09(2021)201 spin s = 4, 5, . . . fields, which would belong to the linear spaces W (s) . An important property of the operator algebra W (c,r) 1 = ⊕ ∞ s=2 W (s) is that it contains the spin 2 field (10.19) which forms the Virasoro subalgebra with central charge

Local IM
The integrals of motion I (a) 1 are built from the densities belonging to the linear space W (2) . They can be expressed as a linear combination of the basic fields with some numerical coefficients: The coefficients may be fixed using the method discussed in section 9.1, which is based on the notion of the reflection operator. Introduce R, defined via its action on the W -currents whereX a andỸ ab are obtained from the fields X a (10.5) and Y ab (10.7) by means of the formal substitution The constants C a are defined as the solution of the linear system, so that under the replacement (10.24), Then the densities entering into the local IM are those, which are invariant under the reflection up to a total derivative: The following r spin 2 local fields

The isotropic limit
The limit was already discussed in the context of the ODE in section (7.4). As prescribed by eq. (7.74) one must first perform the substitution and then take β → 1 − . For the local IM I (a) 1 , we define Consider the local densities (10.29). Following the above procedure with K = r in eqs. (11.1), (11.2), one obtains (11.3) Recall that these operators act invariantly in H j,m,P (6.29). Provided that (11.4) the latter can be realized as a subspace of the tensor product of r integrable representations of the Kac-Moody algebra at level one. It is easy to see that the densities for the local integrals of motion can be expressed in terms of the currents J (11.6) The notation G (a) = (∂φ a ) 2 stands for the Virasoro field associated with the k a = 1 Kac-Moody algebra. Formula (11.5) admits an immediate generalization for arbitrary positive integers k a . One assumes that the local densities are expressed as a quadratic combination of the Kac-Moody currents, subject to the OPE (1.13) The local IM should be invariant w.r.t. the global sl(2) symmetry that occurs in the isotropic limit. Also, taking into account the expression (8.22) for the eigenvalues of the IM, specialized to ξ → +∞, one arrives at the formula It is instructive to consider the limit when all the k a tend to infinity. One writes k a = ν a K and takes K → ∞ while keeping ν a fixed. The latter automatically satisfy r a=1 ν a = 1 . (11.10) This can be interpreted as the classical limit, with K being identified with the inverse Planck constant: The currents j . (11.13) In turn, for the integrals of motion (11.7), . (11.14) Notice that the term G (a) outside the sum in (11.7) is absent from the classical expression. It can be interpreted as an effect of renormalization (quantum counterterm) which appears already at the first perturbative order.

The Gaudin limit
To take the Gaudin limit, following the discussion in section 7.4, the parameters z a should be rescaled as z a = δ y a with δ → 0. Defining the Gaudin Hamiltonians as one obtains (11.16) Note that the spectrum of the first few lowest integrals of motion for the affine Gaudin model associated to sl(N ) was studied in the work [13], while for the case of an arbitrary Lie algebra some conjectures are formulated in ref. [12]. Also, the classical limit of H (a) G results in a similar expression as in the r.h.s. of (11.14) with z a replaced by y a .

General case
For arbitrary positive integer k a and β ∈ (0, 1) we propose the following formula for the Hamiltonians of the GAGM,  (11.19) Recall that the spin 2 field W The result was already presented in the introduction, see eq. (1.16). Taking the limit β → 1 − , the combination 1 gen coincides with the r.h.s. of (11.7). Another consistency check is that the eigenvalue of the operators (11.18), computed on the primary states e j,±2J,P (7.3), agrees with the prediction (8.22) that comes from the ODE side. One can show by direct computation that H (a) gen with a = 1, 2, . . . , r mutually commute.
The generalized affine Gaudin model, being a multiparametric integrable system, admits various interesting limits. In the previous subsection the isotropic limit was discussed as well as the classical one. Let us emphasize that the latter was performed only after taking β → 1 − . It turns out that a meaningful classical limit exists for any fixed value of β and z a . Indeed, setting k a = ν a K and j with g (a) = η AB j and ν 1 + ν 2 + . . . + ν r = 1. Needless to say that these classical observables mutually Poisson commute w.r.t. the Poisson bracket (11.13).

Baxter-type statistical systems
Behind the integrability of the generalized affine Gaudin model are the algebraic structures, which are inherited from the quasi-triangular Hopf algebra U q sl(2) and the universal Rmatrix R. The development of quantum groups, that encompasses these notions, was inspired by the works of Baxter on exactly soluble lattice models. The R-matrix, encoding the Boltzmann weights of the statistical system, comes from taking a particular realization of R. Namely, the finite dimensional matrix where π stands for the 2 + 1 dimensional representation of the U q sl(2) algebra, satisfies the Yang-Baxter equation. Let's suppose that 1 = 1 2 . Denoting the matrix (12.1) is given explicitly by [52] with some overall factor r(λ), which cancels out in the Yang-Baxter equation and is not essential for our purposes. We set r(λ) = 1 and define It is a 2 × 2 matrix, whose elements act in a 2 + 1 dimensional linear space. The lattice transfer-matrix is constructed from such building blocks. First introduce the monodromy matrix (12.6) Its entries act in the "physical space", which is the tensor product (12.7) Taking the trace over the two dimensional "auxiliary space", one obtains the transfer-matrix T(ζ) = Tr ω σ z M (ζ) . (12.8)

JHEP09(2021)201
Here the parameter ω can be arbitrarily chosen. It will be assumed to be a unimodular number, so that ω 2 = e 2πik (12.9) with some real k. Note that the similarity transformation, which is performed in (12.4), has no effect on the trace in the definition (12.8). It was done to make apparent that T(ζ) is a polynomial in ζ = −λ 2 of order ζ N , normalized according to Also, as ζ → ∞ one has Here S z stands for the z projection of the total spin operator, (12.11) which commutes with the transfer-matrix. The Yang-Baxter equation guarantees that T(ζ) commutes with itself for different values of the spectral parameter ζ. The diagonalization problem can be solved using whatever version of the Bethe ansatz approach (e.g., quantum inverse scattering method). In the sector with given value of S z , the Bethe ansatz equations read as [17,53,54] where M = N J=1 J − S z .

Main conjecture concerning the scaling
Being a purely algebraic procedure, the diagonalization of T(ζ) can be performed for arbitrary values of the inhomogeneities η J and the positive integers 2 J . However, to make connection with the field theory where a fundamental symmetry is translational invariance, we'll focus on the lattice with N being divisible by some integer r, N = rL , (12.13) while the parameters η J , J satisfy the r -site periodicity conditions 1, 2, . . . , N ) . (12.14) The lattice model turns out to be critical when q is a unimodular number, |q| = 1. With a suitably defined scaling limit where L → ∞, the universal properties are described by a CFT. However, different types of critical behaviour occur depending on the domain of the parameter q. We conjecture that with an appropriate definition of the scaling limit, and a

JHEP09(2021)201
proper identification of the parameters, the lattice transfer-matrix T(ζ) becomes the operator τ 1 2 (λ) (2.6) as N → ∞. The scaling limit of the Baxter Q-operator (whose eigenvalues are polynomials in ζ with zeroes being the roots of the Bethe ansatz equations (12.12)) yields a + (λ). 12 The anisotropy parameter q should belong to the domain 12.15) and can be written as Also, the positive integers k a from the GAGM are identified as k a = 2 a . (12.17) In taking the scaling limit, one should send ζ → 0 such that the combination ζ N 2 remains fixed. Up to an overall constant, it coincides with the spectral parameter λ 2 entering into τ 1 2 (λ) and a ± (λ). The r inhomogeneities η a are related to the r parameters z a though, in general, the relation is highly non-trivial. 13 It is important to keep in mind that the scaling limit is defined not for the full space of states V N (12.7), but only a certain class of so-called low energy states. A study of the scaling limit, among other things, requires an analysis of the algebraic structures underlying the lattice system. In the case when all the a = 1 2 or, equivalently, k a = 1 the matrix T(ζ) coincides with the one-row transfer-matrix of the original Baxter inhomogeneous six-vertex model [17] subject to quasi-periodic boundary conditions. The recent paper [55] contains a summary of the known results concerning the algebraic aspects of this model. A comprehensive study of the critical behaviour of the lattice system is beyond the scope of this work. Here we would merely like to provide an illustration of how the GAGM appears in the scaling limit. As such we'll stick to the case a = 1 2 and rely on the results of [55]. For the reader's convenience, we follow the notations of that work.
In the case of the homogeneous six-vertex model the transfer-matrix commutes with the spin 1 2 Heisenberg XXZ Hamiltonian. A similar important property holds true for the model, where the parameters η J satisfy the periodicity conditions (12.14) with any r ≥ 1. Namely, the commuting family contains spin chain Hamiltonians that are given by a sum of terms, each of which is built out of local spin operators supported on r + 1 consecutive sites of the lattice. There are r such Hamiltonians and, in terms of the one-row transfer-matrix, they are expressed as (see eq. (6.11) in ref. [55]) 1, 2, . . . , r ) . (12.18) 12 There are two Baxter Q-operators Q±(ζ). The zeroes of the eigenvalues of Q+ satisfy (12.12), while those of Q− obey the system, which is obtained from (12.12) via the substitutions S z → −S z and ω → ω −1 (for details see, e.g., ref. [55]). In what follows we always assume that S z ≥ 0, and discuss the zeroes of the eigenvalues of Q+, rather than Q−. 13 It may also be that some restrictions on the domain of the inhomogeneities ηJ need to be imposed in order to make connection with the GAGM. This point requires further investigations.

JHEP09(2021)201
For the definition of the scaling limit, a central rôle belongs to the sum H (a) , (12.19) which is essentially the logarithmic derivative of a r row transfer-matrix. The ground state of this Hamiltonian serves as the reference state from which the energy of the excited states is counted. In performing the scaling limit one takes N → ∞ but considers only the class of states, whose excitation energy over the ground state energy is sufficiently low. Then as N → ∞ the low energy spectrum organizes into the conformal towers, which are classified w.r.t. to the algebra of extended conformal symmetry. In the case at hand, it is expected that the latter coincides with W . The two W -algebras in the tensor product are isomorphic and describe the left and right chiralities of the underlying CFT. In turn, for the low energy spectrum of H (12.19) at large N , one would have . (12.20) Here I 1 is the local IM in the generalized affine Gaudin model, defined via (6.36), (6.33).
The operator I 1 , corresponding to the other chirality, is described by the similar equations. The first term in the r.h.s. is proportional to the identity operator and so does not depend on the particular low energy state. It contains the constant e ∞ , which stands for the specific bulk energy. The explicit form of e ∞ is not significant in the context of the field theory. The positive constant v F , usually referred to as the Fermi velocity, depends on the overall normalization of the Hamiltonian. The asymptotic formula (12.20) involves a contribution from both the left and right chiralities. In order to study the scaling limit of a low energy state, it is useful to employ the relations which involve each chirality separately. Among these are the sum rules for the Bethe roots. Let's illustrate them on the example of the inhomogeneous six-vertex model with anisotropy parameter q = − e iπ 2 (β 2 −1) , subject to the two site periodicity condition. We focus on the case when η 1 and η 2 are unimodular numbers and, without loss of generality, set Also, it will be assumed that The number of lattice sites N should be even. Suppose for now that N is divisible by 4 and S z = 0. Then the ground state of the Hamiltonian H (12.19) is a singlet. Together with ζ n it is convenient to use The roots θ (±) n with n N or N/4 − n N will be referred to as the left and right edge roots, respectively, while the rest will be called the bulk roots. The latter, as N → ∞, become densely distributed along the lines m(θ) = ± 1 2 . Namely, introducing In contrast the edge roots develop a scaling behaviour. It turns out that, keeping either n or N/4 − n fixed as N → ∞, the following limits exist A consequence of our general proposal is that, up to an overall factor, s (±) n coincide with the zeroes of the spectral determinant D j,m,A (µ) for the vacuum ODE (7.7) with As usual, β 2 = ξ ξ+1 and A = i β −1 −β √ K P − 1 2 β m while, for the case under consideration, The ODE involves the parameters z 1 , z 2 and we are free to set z 1 z 2 = 1. Provided the inhomgeneities are restricted as in eqs. (12.21) and (12.22), they turn out to be unimodular numbers and it is convenient to swap them for x, With such reality conditions imposed on z 1 and z 2 , the properties of the zeroes of the spectral determinant are described in appendix B. In particular, they are split into two groups µ   with α 0 (x) from (B.9). Comparing the above with eq. (12.28) one finds that α 0 (x) = . (12.33) This allows one to relate the inhomogeneities of the lattice system η 1 and η 2 with z 1 and z 2 . Notice that the domains ∈ π 2 (1 − β 2 ), π 2 and ∈ π 2 , π 2 (1 + β 2 ) map to x ∈ [ 0, 1) and x ∈ (−1, 0 ], respectively. Moreover, a comparison of (12.29) with the asymptotic formula (B.12) yields s (±) n = f 1−β 2 µ (±) n . (12.34) Here f = f (x), given by (B.11), is a real and positive number.
Relation (12.34) can be checked using the sum rules. For any given solution to the Bethe ansatz equations (12.12) consider the finite sum For the class of low energy states it turns out that the following limit exists . (12.36) Without going into details, we just mention that the existence of the limit comes from the Bethe ansatz equations. Our main conjecture implies that for a low energy state, h (∞) 1 appears in the first term of the Taylor expansion of the corresponding spectral determinant.
In view of (9.60), as well as the µ − λ relation specialized to K = 2 (see eq. (9.57)), the latter can be expressed in terms of the eigenvalues of the non-local integrals of motion H (+) n : .
is the eigenvalue of the non-local IM corresponding to the second chirality. The sum rules (12.38) and (12.40) can be applied to the Bethe roots corresponding to the ground state of the Hamiltonian H in the sector S z = 0 with even N/2. We found that in the scaling limit it becomes the statē The numerical data in support of the identification is presented in figure 2.

JHEP09(2021)201
Consider the formula describing the scaling limit of the Hamiltonian (12.20) specialized to the ground state. The extensive part of the energy is determined by the density of the bulk roots (12.26). A standard computation results in the expression (12.43) As was already mentioned, the Fermi velocity depends on the overall normalization of the Hamiltonian. With the definition (12.18) and (12.19), it turns out that v F does not depend on the inhomogeneities and is given by From eqs. (12.42) and (12.20), one has the prediction The latter has been confirmed via numerical work (see figure 2). For odd N/2 with S z = 0 the ground state is doubly degenerate. Again, using the sum rules, one finds that the states in the scaling limit can be identified with (some details can be found in figure 3) . (12.46) In this case, There is a doublet of first excited states, whose energy in the scaling limit is described by the same formula. Clearly, they correspond tō In making the above identifications it was assumed that 0 < x < 1. In the domain −1 < x < 0, the picture is reversed. The ground states flow to (12.48), while the first excited states become (12.46). For even N/2 and S z = 0 the first excited state is also two-fold degenerate. Applying the sum rules, one finds that in the scaling limit they becomē  As for the energy, one has N/2 − even, S z = 0 : lim (12.51) When studying the scaling limit, it is apparent how to assign an N dependence to the ground state and perhaps the first few excited states in the disjoint sectors of the Hilbert space. These sectors, in the case under consideration, are distinguished by the value of the quantum number S z and the parity of N/2. However, for a general lattice system there are difficulties in introducing the N -dependence, i.e., the Renormalization Group (RG) flow for an individual stationary state. Of course, since the Hilbert space is not isomorphic for different lattice sizes, the problem only makes sense for the low energy part of the spectrum. But even then, forming individual RG flow trajectories is not a trivial task. For integrable spin chains, there is a procedure for assigning an N dependence to a Bethe state that relies explicitly on the Bethe ansatz equations (for a discussion see, e.g., ref. [56]).
In the general set-up with S z ≥ 0 and both even or odd N/2, the low energy stationary states in the scaling limit can be classified according to the irreps of the W

JHEP09(2021)201
The admissible values of P NS andP NS are given by (12.53) where − 1 2 ≤ k < 1 2 ; w = 0, ±1, ±2, . . . . (12.54) In turn, the Ramond states are described as ej ,m,P R (γ;w) ⊗ e j,m,P R (γ; w) withj, j ∈ 1 2 , 0 , 0, 1 2 , m = 1 . (12.55) Here (12.56) while 0 ≤ k < 1 ; w = 0, ±1, ±2, . . . . (12.57) The following comment is in order here. The transfer-matrix and the Hamiltonian essentially depend only on the unimodular parameter ω 2 (12.9). Taking the logarithm results in log(ω 2 ) = 2πi (k + w) (12.58) with w ∈ Z. The parameter k can be brought to lie in any segment of unit length, which is usually referred to as the "first Brillouin zone". Formulae (12.54) and (12.57) mean that for the Neveu-Schwarz and Ramond sectors, the first Brillouin zone should be chosen in a different way, to be − 1 2 , 1 2 and [ 0, 1), respectively. The integer w (the winding number), which labels the bands of the spectrum, is a quantum number that appears only in the scaling limit. It is worth reiterating that for finite N this quantum number is not well defined.
Let |Ψ N be a trajectory describing the RG flow of some low energy Bethe state. For fixed N it is an eigenstate of the matrices H (a) (12.18) and we denote by E with v F = 2 (1 + ξ), e ∞ as in (12.43) and I  values of P ,P are given by eqs. (12.53) and (12.56) in the Neveu-Schwarz and Ramond sectors, respectively. Notice that the relation (12.59) does not depend on the inhomogeneity parameters η a . The inhomogeneous six-vertex model with r = 2 and η 1 = −η 2 = i possesses an additional global Z 2 symmetry. This model was studied in the work [51]. Among other things, based on the modular invariance of the CFT partition function (for k = 0), the authors discuss the gluing conditions for the right and left chiralities. Their results are also applicable to the more general case with the inhomogeneities as in eqs. (12.21) and (12.22).
Contrary to the eigenvalues of the operator H (12.19), i.e., E N + E N the scaling part of the eigenvalues of each H (a) individually contain an explicit dependence on the inhomogeneities. One can obtain the asymptotic formula where ξ = β 2 1−β 2 and f = f (x) is given by (B.11). Note that , ∂ x f (0) = 0 . (12.62) An illustration for the Bethe states, which flow to the tensor product of the Ramond primary states (12.49), is provided in figure 4.

JHEP09(2021)201
To determine the scaling limit of |Ψ N one could employ formulae (12.59) and (12.61). However, unlike the sum rules, they contain the contribution of the eigenvalues of both right and left IM I . Nevertheless, their computation is a much simpler task than that of the eigenvalues of the non-local IM. Perhaps the most convenient way of identifying the scaling limit of an RG trajectory involves the reflection operator. Its eigenvalues admit a closed expression (8.4), while the construction of the eigenstates is a straightforward algebraic procedure. On the other hand, the eigenvalues of the left and right reflection operators appear separately in the large N limit of certain products over the Bethe roots. Some related discussion can be found in section 11 of ref. [56]. 14

Outlook
Among the lessons we took away from carrying out this project is that the generalized affine Gaudin model fits within the standard framework of Yang-Baxter integrability. As a result it can be studied using conventional techniques such as the quantum inverse scattering method and its variants. This has exciting implications, since the model shares the same integrability hierarchy as a variety of physically interesting systems. The latter includes a large class of classically integrable Non-Linear Sigma Models (NLSM). For example, in refs. [57][58][59] the classical affine Gaudin model (11.14) and its generalizations to Lie algebras of higher rank are used in assembling integrable NLSM. The results of our work provide an avenue for the systematic quantization of such theories within the ODE/IQFT correspondence, in the spirit of ref. [36]. Another application lies in the domain of Condensed Matter Physics. The GAGM may be understood as an integrable multiparametric generalization of the Kondo model. Finally, we hope that the output of this work would be useful in the context of the representation theory of infinite dimensional algebras.

JHEP09(2021)201
Comparing the above equation with the WKB formula (A.14), one finds The and R j,m,A is given by (8.4). The asymptotic coefficient q −1 is expressed in terms of an integral, which converges only for 0 < (1 + ξ)K < 2. Keeping in mind that K is a positive integer, it can be literally applied to the case K = 1 and 0 < ξ < 1. Nevertheless, the integral for q −1 can be analytically continued from the domain of convergency. This is done by replacing the integration over the positive real semiaxis by that over a certain contour integral. To illustrate the procedure, consider the case r = 1 when The integration contour C is the image of C ζ under the conformal map (A.25). Alternatively it can be described as follows. Consider the Riemann sphere under stereographic projection such that the south and north poles are mapped to z = 0 and z = ∞, respectively. Then C is the image of a Pochhammer loop on the sphere that winds around the two poles. The homotopy class of that loop is schematically depicted in figure 1. Notice that ν 0 and ν ∞ are the exponents, which appear in the leading behaviour of P(z) at z = 0 and z = ∞: For the general case with r > 1, the analytic continuation of q −1 is given by equations (A.26) and (A.27). However one should take care that the integration contour C does not wind around the points z = z a . This can be achieved by taking C to lie inside the domain D, which is the union of |z| > ε −1 , |z| < ε and the wedge | arg(z)| < δ such that z a ∈ D. Also in the case r > 1 the branch of P(z) can not be chosen to be real for positive real z. Instead it is sufficient to require that P(z) asymptotically approaches to positive real values as z → +∞. Under the above conditions formulae (A. 26 A breakdown of the leading asymptotic formula for large complex µ is related to the lines along which the zeroes of the entire function D j,m,A (µ) accumulate asymptotically. In the case under consideration the positions of the Stokes lines significantly depend on the domain of x. Below we take x 2 to be a real number.