Kondo line defects and affine Gaudin models

We describe the relation between integrable Kondo problems in products of chiral SU(2) WZW models and affine SU(2) Gaudin models. We propose a full ODE/IM solution of the spectral problem for these models.

These are families of mutually commuting line defects parameterized by a conformal symmetry-breaking scale e θ ≡ λ −1 .

JHEP01(2022)175
We will begin by a quick review of the affine Gaudin model in section 2. We then define λ-opers with singularities of trivial monodromy and derive the affine Bethe equations in section 3. We analyze the Stokes data at large and small λ in sections 5 and 6, respectively, and compare it with direct calculations for the Kondo defects. In the large λ regime, the Stokes data are obtained with the help of (exact) WKB analysis [28][29][30][31][32][33][34][35][36][37][38]. The examples under study have some unusual features that require us to generalize the standard WKB analysis in order to evaluate the complete collection of the Stokes data. In order not to clutter the main body, we will only quote the results in section 6 while leaving the detailed review of the WKB analysis and our generalizations in the appendix B.
We will also see that the construction automatically includes integrable defects in coset . We discuss briefly some alternative semiclassical limits in section 7.
Although we focus on SU (2) examples in the main body of the paper, we expect the results to extend to general affine ADE Lie algebras and will comment on that in section 8.
We conclude the paper with a list of open questions in section 9.

Affine Gaudin models, classical and quantum
The affine Gaudin model, first studied in [21], is a somewhat conjectural integrable system which quantizes the classical affine Gaudin model, in a manner analogous to the relation between the classical and quantum KdV integrable systems [39][40][41]. 1 The classical affine Gaudin model is defined by a collection of Poisson-commuting Hamiltonians built from classical currents J a i with Kac-Moody Poisson brackets. The latter are organized into a Lax matrix where the z i are couplings and we use the auxiliary 1-form where k i are the levels for the currents J a i . The Lax matrix is used to define both the Poisson-commuting transfer matrices and families of local Hamiltonian densities H (n) u (σ) labelled by exponents n of the affine Kac-Moody algebra and zeros ζ u of the twist function (2.2). In classical types, the H (n) u are given by specific homogeneous polynomials [42][43][44][45] in Tr ϕ(z)L a (z; σ)t a r z=ζu , (2.4) of total degree n + 1 in the currents J a i .

JHEP01(2022)175
One of our main proposals is that the correct quantization of the affine Gaudin model involves Kondo line defects defined by coupling a spin to a collection of quantum Kac-Moody currents J a i . These Kondo lines are defined in the UV in the same manner as the classical transfer matrices: (2.5) The couplings g i [z] need to be renormalized and acquire a scale dependence. It was conjectured in [19] that the RG flow factors through a flow of the spectral parameter, so that the spectral parameter may be identified with the (complexified) renormalization scale e θ via dimensional transmutation. In other words, the RG flow defines a one-parameter family of commuting line defects. The specific functional form of the couplings g i [z] along the commuting family depends on the chosen renormalization scheme. The RG flow is physically rich and depends sensitively on the spin of the auxiliary sl 2 generators t a and on the relative UV couplings. The endpoints is some IR-free line defect whose nature can be predicted with the help of the ODE/IM correspondence. For special choices of spin and couplings, the endpoint is an irrelevant deformation of a single identity line defect. Such a deformation must take the form for some collection of bulk quasi-primaries O (n) u (σ) of dimension n + 1. Now we denote with u the choice of UV line defect flowing to the identity line. We will see that this generalizes naturally the choice of a zero ζ u for ϕ(z) above. A special property of such deformations by bulk chiral currents is that there exists a renormalization scheme where the path-ordered exponential becomes effectively integration along separate contours. Therefore the path-ordered exponential is essentially Abelian, and reduces to the exponential of the zero-modes of the O

Opers, λ-opers and affine opers
In this section, for simplicity, we specialize to the case of sl 2 . The generalisation to sl 3 will be discussed in section 8. The main objective of this section is to introduce the family of Schrödinger operators which provides the conjectural full ODE/IM solution of the Kondo defects spectrum problem: where P (x) = e 2αx N a=1 (x − z a ) ka and t(x) is an auxiliary meromorphic classical stress tensor which will be determined by a solution of the Bethe equations.

JHEP01(2022)175
We will motivate some of our definitions in analogy to the well-known correspondence between the (non-affine) Gaudin model and opers with singularities of trivial monodromy [47][48][49]. The latter is one of the most basic manifestations of the Geometric Langlands correspondence and can be investigated with the help of supersymmetric gauge theory [50,51]. It would be very nice to give a similar derivation of the ODE/IM proposal based on gauge theory or string theory constructions.

sl 2 opers
An sl 2 oper is a complexified Schrödinger operator with a natural transformation law under a change of coordinate This implies that t(x) transforms as a classical stress tensor We will always consider sl 2 opers for which t(x) is a rational function and only allow coordinate transformations which preserve this property. We will typically denote a solution/flat section of the Schrödinger equation as ψ(x): (3.5) and the (constant) Wronskian of two solutions as The data of an sl 2 oper (3.2) is equivalent to that of a flat connection More generally, see for instance [52], an sl 2 oper can be described as a flat connection of the form where a(x) and b(x) are rational functions, modulo gauge transformations by unipotent upper-triangular matrices whose entries are rational functions. We can fix the gauge invariance completely by bringing (3.8) to its unique canonical form (3.7) with stress tensor given by (3.9) JHEP01(2022)175

sl 2 λ-opers
An sl 2 λ-oper, or simply λ-oper, is a complexified Schrödinger operator with standard dependence on a quantization parameter , here denoted as λ, namely (3.10) The coordinate transformations act in the same way as for an sl 2 oper, so that t(x) is a classical stress tensor and P (x) is a quadratic differential. We will always work with λ-opers for which t(x) is a rational function on C but allow P (x) to be a more general analytic function, typically with branch points or an essential singularity at infinity, but whose logarithmic derivative is a rational function. The data of a λ-oper can be encoded in a flat connection of the form More generally, an sl 2 λ-oper is a flat connection where a(x) and b(x) are rational functions, modulo gauge transformations by uppertriangular matrices of the form 1 v(x)λ 0 1 , (3.13) for some rational function v(x). Every λ-oper admits a unique canonical form as in (3.11) with stress tensor as in (3.9). Of course, we could equally describe a λ-oper using a flat connection of the form . (3.14) This leads to another (equivalent) way of describing λ-opers, namely as a flat connection (3.15) modulo gauge transformations by lower-triangular matrices of the form (3.16)

Miura opers and singularities of trivial monodromy
A Miura sl 2 oper, or simply Miura oper for short, is a connection of the form (3.8) with b(x) = 0, namely . (3.17)

JHEP01(2022)175
Its equivalence class modulo gauge transformations by unipotent upper-triangular matrices defines an oper, with stress tensor t(x) = a(x) 2 + ∂ x a(x), which we refer to as the oper underlying (3.17). The corresponding Schrödinger operator (3.2) factorises as (∂ x + a(x))(∂ x − a(x)) and has an obvious solution ψ(x) = e a(x)dx which is an eigenline of the monodromy around each singularity of t(x).
It is useful to allow the Miura oper to have apparent singularities where the monodromy eigenline ψ(x) has a simple zero but where t(x) is regular. At such an apparent singularity, a(x) behaves as . (3.18) Another important type of singularity is one of the form for a non-negative half integer l. Then 3.20) and the Schrödinger operator (3.2) has a local solution with ±1 monodromy around z, of the form ψ(x) ∼ (x − z) −l + · · · (3.21) The Miura condition then gives a second local solution with ±1 monodromy around z, of the form ψ (x) ∼ (x − z) l+1 + · · · (3. 22) It follows that the monodromy of any flat section around z must be ±1. Therefore z is a regular singularity of trivial monodromy for t(x). One can also see this by noting that the Miura oper (3.17) is gauge equivalent to the connection where r(x) = a(x) + l x−z , which is manifestly regular at z for non-negative l. A quick discussion of the term "trivial monodromy" is in order here. If l is allowed to be half-integral, we have to consider the monodromy as living in PSL (2), so that ±1 is a trivial monodromy. If l is restricted to be integral, then we can take the monodromy to be valued in SL (2), and will still be trivial. This binary choice is a manifestation of Geometric Langlands duality: PSL (2) opers are dual to the SL(2) Gaudin model, and viceversa.

Miura λ-opers and singularities of trivial monodromy
A Miura sl 2 λ-oper, or simply Miura λ-oper, is a connection of the form . (3.24)

JHEP01(2022)175
This is of the general form (3.12) and therefore a Miura λ-oper defines a λ-oper with stress tensor t(x) = a + (x) 2 + ∂ x a + (x). We refer to this as the λ-oper underlying (3.24). It can be described as a complex Schrödinger operator Crucially, the connection (3.24) is locally gauge equivalent (in PSL (2)) to a connection of the following alternative form where a + (x) + a − (x) = − 1 2 ∂xP (x) P (x) . We refer to this connection as the dual of the Miura λ-oper (3.24). Since it is of the general form (3.14) it also defines a λ-oper, which we call the dual λ-oper underlying (3.24), with stress tensort(x) = a − (x) 2 + ∂ x a − (x). The latter can also be described as a complex Schrödinger operator of the same form as in (3.25) with a + (x) → a − (x).
Since the Miura λ-oper (3.24) and its dual (3.26) are gauge equivalent, we therefore identify a crucial property of λ-opers: the pair of λ-opers with stress tensors built from a ± (x), i.e. the λ-oper and the dual λ-oper underlying a given Miura λ-oper, have the same monodromy (in PSL(2), unless P (x) is a perfect square).
If at some generic point w we have then it follows by the above arguments for singularities of Miura opers that the Miura λ-oper built from a + (x) has an apparent singularity at w while the other Miura λ-oper built from a − (x) has a regular singularity at w, which must necessarily have trivial monodromy (in PSL (2)). The same argument applies if at a point w we have with the roles of the two Miura λ-opers (3.24) and (3.26) interchanged. If at a zero z of order k of P (x) we have As long as 0 ≤ l ≤ k 2 , the pair of Miura λ-opers both have trivial monodromy around z. Indeed, the Miura λ-oper (3.24) is gauge equivalent to the connection where we wrote P (x) = (x − z) k q(x) with q(z) = 0 and r(x) = a + (x) + l x−z , which is manifestly regular at z when 0 ≤ l ≤ k 2 . A quick discussion of the term "trivial monodromy" is again in order here. If l is allowed to be half-integral and k integral, we have to consider the monodromy as living in PSL(2), so that ±1 is a trivial monodromy and gauge transformations can have a sign ambiguity. If l is restricted to be integral and k even, then we can take the monodromy to be valued in SL (2). This binary choice is presumably a manifestation of an affine Geometric Langlands duality: PSL(2) λ-opers are dual to the affine SL(2) Gaudin model, and viceversa.

Opers with singularities of trivial monodromy and Bethe equations
For a general oper, the condition for a regular singularity to have trivial monodromy is an intricate polynomial constraint on the coefficients of the expansion of t(x) near the regular singularity.
Given a Miura oper on C with a rank 1 irregular singularity at infinity and whose other singularities are all regular with trivial monodromy, we can write The condition that each w i is an apparent singularity reduces to the Bethe equations These are the Bethe equations for a sl 2 quantum Gaudin model with sites of spectral parameters z a , supporting sl 2 irreps of dimension 2l a + 1.
We call the overall residue of a(x) at infinity the weight at infinity of the Miura oper. Since the underlying oper has trivial monodromy at all the z a and w i we refer to it as an oper with singularities of trivial monodromy. The eigenvalues of the quantum Gaudin Hamiltonians can be extracted from the expression of t(x).

λ-opers with singularities of trivial monodromy and affine Bethe equations
We are interested in the class of Miura λ-opers on C for which a ± (x) take the same form and satisfy the Bethe equations

JHEP01(2022)175
These ensure that the (dual) Miura λ-oper built from a + (x) (resp. a − (x)) has apparent singularities at each w i (resp. w i ). The weight at infinity of the Miura λ-oper is the pair of residues of a ± (x) at infinity. These are the Bethe equations for an affine sl 2 quantum Gaudin model with sites of spectral parameters z a , supporting sl 2 Weyl representations induced from sl 2 irreps of dimension 2l a + 1, for WZW current algebras of level k a .
The λ-oper (resp. the dual λ-oper) underlying a given Miura λ-oper has regular singularities with trivial monodromies at the zeroes z a of P (x) as well as at w i (resp. w i ), for all values of λ. We therefore refer to the λ-oper and its dual as a pair of λ-opers with singularities of trivial monodromy. They have interesting Stokes data at x = ∞ which we call the monodromy data of the pair of λ-opers.
The eigenvalues of the affine sl 2 quantum Gaudin model transfer matrices, as well as the quantum local Hamiltonians, can be extracted from the Stokes data in a manner described in the remainder of the paper.

Weyl reflections
Given some sl 2 oper with trivial monodromy t(x) = a(x) 2 + ∂ x a(x) and α = 0, there must be two canonical solutions which at infinity behave like e ±αx times some analytic functions. The one behaving as e −αx is the Miura eigenline, with logarithmic derivative a(x). The other gives a second rational solutionã( with opposite weight at infinity to a(x). This gives an action of the Weyl group of sl 2 on the collection of Miura opers with the same stress tensor t(x). In particular, it acts as a Weyl transformation on the weight at infinity of the Miura oper.
In terms of connections of the form (3.17), the above transformation a(x) →ã(x) is implemented as a gauge transformation by a unipotent upper-triangular matrix which preserves the Miura form of the connection. Explicitly, a gauge transformation of the Miura oper (3.17) by If we have a Miura λ-oper with trivial monodromy, such that α ± are sufficiently generic, then we have two transformations, α + → −α + or α − → −α − , which map it to a different Miura λ-oper, with the same P (x) and the same monodromy data, as either one of the λ-opers is fixed by the transformations.

JHEP01(2022)175
This produces a new pair of Miura λ-opers with a + (x) →ã provided that the functions f ± (x) are (rational) solutions of the Riccati equation These two reflections can be iterated to generate an interesting group: the Weyl group of sl 2 . It acts as a Weyl transformation on the weight at infinity of the Miura affine oper. More precisely, repeated reflections act as and similarly on the weight at infinity.

Conjectural count of Bethe solutions
For generic values of α, the relation between opers with singularities of trivial monodromy and the Gaudin model suggests that the number of solutions of the Bethe equations should coincide with the graded dimension of the Gaudin Hilbert space, which is the product of sl 2 irreps of dimension 2l a + 1, graded by total weight. A priori, this statement is rather not obvious. We expect a similar statement for the affine opers with singularities of trivial monodromy: for generic values of α the number of solutions of the Bethe equations should coincide with the graded dimension of the affine Gaudin Hilbert space, which is the product of sl 2 Weyl representations induced from sl 2 irreps of dimension 2l a + 1, for WZW current algebras of level k a .

Special values of α ±
The Weyl reflection is not well-defined for an oper with singularities of trivial monodromy when α = 0, essentially because there isn't a canonical choice of a second solution. Any choice of solution will do, so we really get a CP 1 family of opers with singularities of trivial monodromy. Only one of these solutions is special, in the sense that it decreases at infinity faster than the others, and will thus have a special weight.
In the dual Gaudin model, we are turning off a parameter which breaks the global sl 2 symmetry. A whole sl 2 irrep of eigenstates is represented by a single special oper with singularities of trivial monodromy.
Something similar happens if α + = n(α + + α − ) for any integer n. One of the Weyl reflections in the chain breaks down, and instead we get a continuous family of solutions. In the dual affine Gaudin model, we are restoring one of the sl 2 subgroups of the total affine sl 2 symmetry. A whole sl 2 irrep of eigenstates is represented by a single special affine oper with singularities of trivial monodromy.
If we set both α ± = 0, the whole Weyl chain breaks down and we get an intricate continuous family of solutions. In the dual affine Gaudin model, we are restoring the whole total affine sl 2 symmetry. Essentially, the transfer matrices commute with the total Kac-Moody currents and thus secretly live in the coset CFT.

WKB expansion and quasi-canonical form
The Miura λ-oper (3.24) and its dual (3.26) are locally gauge equivalent to a connection of the more symmetric form . Following [20,21], we will refer to this as a Miura sl 2 oper. We can consider the equivalence class of such a connection under gauge transformations by matrices of the form for some formal power series where u n (x) and v ± n (x) are rational functions. This defines an affine sl 2 oper, or more precisely an sl 2 oper for the untwisted affine Kac-Moody algebra sl 2 associated with sl 2 . We are working in the loop realisation of sl 2 associated with the principal Z-gradation.
The relationship between the different opers described above is depicted in figure 1. Note that the notions of λ-oper and dual λ-oper are naturally associated with the two roots of the Dynkin diagram of sl 2 .
By allowing gauge transformations as in (3.42) one can bring the connection (3.41) to a quasi-canonical form [46] for some formal Laurent series Unlike the canonical form (3.7) of an sl 2 oper as in (3.8), however, the quasi-canonical form (3.44) of an affine sl 2 oper is not unique. Indeed, the quasi-canonical form is preserved JHEP01(2022)175 by residual gauge transformations of the form (3.42) with u(x; λ) = 0 and v − (x; λ) = v + (x; λ), the effect of which is Now the quasi-canonical form (3.44) can be transformed to and by a further gauge transformation we can bring it back to the form In particular, by comparing this with the expression (3.11) for the λ-oper underlying the Miura λ-oper we started with, we recognize the equation for the WKB momentum which is used to study the λ → 0 limit of the transport data of the λ-oper.
In particular, the contour integrals control the WKB asymptotics of certain transport coefficients. They are also known to match the eigenvalues of local integrals of motion for the affine Gaudin model [25,46,53].

Some Kac-Moody conventions
We follow the convention from [19]. Our normalization convention for the spin basis of sl 2 is which satisfy the relations The relations in the corresponding untwisted affine Kac-Moody algebra sl 2 read for n, m ∈ Z. Let |l, κ denote the ground state in the spin l module at level κ.

Action of the spectral flow
Spectral flow [54,55] is an automorphism of sl 2 given, for α ∈ R, by There is also an involutive automorphism induced by the Weyl group W(sl 2 ) = Z 2 which satisfy We therefore have Aut( sl 2 ) = R Z 2 . In particular, the even part is inner and corresponds to the affine Weyl group W( sl 2 ) = (2Z) Z 2 . Consequently, the induced action by U 2Z maps each integral highest weight representation into itself, whereas more general U α maps between the (twisted) modules. For example, (4.9)

The Bethe equations and Bethe vectors
Let us take the quadratic differential P (x) = e 2x x κ . This should correspond to an affine Gaudin model with a single site, i.e. an integrable Kondo problem in the SU(2) κ WZW model. A generic state in the spin l module of the SU(2) κ WZW model which is singular under the zero-mode sl 2 subalgebra can be described by a pair of Miura λ-opers of the form with the Bethe roots w i and w i satisfying the Bethe equations It is useful to denote the set of all Bethe roots {w i } ∪ {w i } collectively as {t i }. To each Bethe root w i we associate the lowering operator F w i = J − 0 in sl 2 and to each Bethe root w i the lowering operator F w i = J + −1 . By analogy with the finite-dimensional case [56] and based on the expression for the Bethe vector in affine Gaudin models with regular singularities [57],

JHEP01(2022)175
we then conjecture that, associated with each solution of the Bethe equations (4.11) with m Bethe roots, there is a corresponding Bethe vector in the spin l module given by . (4.12) This state is singular under the zero-mode sl 2 subalgebra. Moreover, its Virasoro level is equal to the number #w of Bethe roots w i and its spin is l + #w − #w.
The leading non-trivial term in the UV expansion of the transfer matrix is proportional to the zero-mode quadratic Casimir T (2) = J a 0 J a 0 . Since (4.12) has definite spin it is an eigenvector of T (2) with eigenvalue 2(l + #w − #w)(l + 1 + #w − #w). From the subleading term in the UV expansion we obtain the operator We checked that this non-trivial operator is indeed diagonalized by the examples of Bethe vectors in subsection 4.3. It would be very nice to derive the Bethe equations directly from the diagonalization of T (3) with the Bethe vector ansatz (4.12). Moreover, we conjecture that the expectation value of T (3) in the generic eigenstate (4.12) is (4.14) A similar conjecture was made in [22] for the eigenvalue of the first non-local integral of motion of the affine sl 2 Gaudin model describing quantum KdV theory as a coset CFT. We also expect from the UV expansion of the corresponding λ-oper in appendix A.2 that the eigenvalues of all the higher order UV expansion coefficients T (n) of the transfer matrix are given by supersymmetric polynomials in the Bethe roots {w i } ∪ {w j }. This was conjectured in [22] for the higher non-local charges of quantum KdV theory.
The expression for the next subleading term T (4) is much more complicated. We did check in appendix A that the Bethe vector expectation value of the UV expansion of transfer matrices matches the Stokes data of the corresponding λ-opers.
The multichannel case can be treated similarly. For the product of WZW model we take the quadratic differential P (x) = e 2x a (x − z a ) κa . We can then describe a generic state in the tensor product of spin l a modules which is singular under the total zero-mode sl 2 subalgebra using a pair of Miura λ-opers of the form

JHEP01(2022)175
where w i and w i are the Bethe roots satisfying the Bethe equations If we denote the set of all Bethe roots collectively as {t i } then we conjecture that the corresponding Bethe vector in the tensor product of spin l i modules is given by where the sum is over all partitions of the set {t i } into N ordered subsets (t a,1 , . . . t a,ma ) with m 1 + . . . + m N = m. It has Virasoro level #w and spin a l a + #w − #w.
In the two-point case, the leading term in the UV expansion of the transfer matrix is proportional to the total zero-mode quadratic Casimir (J a 0,1 + J a 0,2 ) 2 with eigenvalue 2(l 1 + l 2 + #w − #w)(l 1 + l 2 + 1 + #w − #w) on the Bethe vector (4.17). The next term in the UV expansion is proportional, in a suitable renormalization scheme, to We conjecture that its expectation value in the eigenstate (4.17) is given by the supersymmetric polynomial in the Bethe roots

Examples
In this subsection we study the solutions to the Bethe equation (4.11) in the vacuum and spin 1 2 WZW modules.
In particular, when κ = 0 it is a descendant of the singular vector J + −1 |0, 0 . The situation when κ = −1 is more subtle since we see that the state J + −2 |0, −1 in the spin 0 module of level −1 is not described by a solution of the Bethe ansatz. In the limit κ → −1 of a solution of the Bethe equations for κ = −1, all the Bethe roots collide with the origin so that the Miura λ-oper (4.10) becomes This is no longer of the form (4.10) since the residues at the origin are not given by the pair (−l, − κ 2 + l) = (0, 1 2 ) corresponding to the highest weight of the vacuum module at level κ = −1. However, the residues of (4.19) do correspond to a shifted Weyl reflection of this highest weight. Indeed, the simple roots of sl 2 act by shifted Weyl reflections on the residues at the origin as · · · ←→ (0, 1 2 ) ←→ (1, − 1 2 ) ←→ (−1, 3 2 ) ←→ · · · Therefore the pair (4.19) describes a generalised Miura λ-oper in the sense of [20,49]. We conjecture that the state J + −2 |0, −1 is described by this generalised Miura λ-oper. This is checked to fourth order in the UV expansion in appendix A.
• For κ = 1 we have just one spin 3 2 Bethe vector • For κ = − 1 2 we also have just one spin 3 2 Bethe vector • For κ = −2 there are no spin 3 2 Bethe vectors. #w = 2, #w = 2: there are again two inequivalent families of solutions: (i) The first family is valid for κ = −2, − 7 3 and the corresponding Bethe vector is proportional to This vanishes when κ = −2. On the other hand, when κ = − 7 3 we find which has zero norm.
(ii) The second family is valid for κ = 1 with Bethe vector proportional to
• For κ = 1 we have just one spin 1 2 Bethe vector we also have just one spin 1 2 Bethe vector

The UV expansion
Kondo line defects in i SU(2) k i WZW models are defined as the trace of the path ordered exponentialT where t a are the generators of the Lie algebra su(2) and the trace is taken in an su(2) representation R, labeled by its dimension n from now on. The factor i in front of the integral is √ −1. The WZW currents are denoted as J a i (σ) for each SU(2) factor. The integration is along the compact direction.
Following [19], we adopt the convention that the physical RG flows start from asymptotically free defects and the couplings grow in the positive real direction approaching the infrared. Therefore the UV expansion is concerned with small positive g i . Perturbatively in g i , we can expand the exponential where eachT (N ) n depend on the set {g i }, and perform the loop integral. In doing so, one carries out a careful and lengthy renormalization procedure since the currents don't commute with each other. This was done [19] following the prescription given in [17]. We are interested in the expectation value of the Kondo line operator in a generic state | , with denoting the list of quantum numbers of the state. Details can be found in appendix A. If the twist α + is nonzero, the Kondo line defect comes with a twist,

JHEP01(2022)175
Therefore the leading order is simply given by the character in the representation R of dimension n, namely Note that in order for (5.3) to make sense, the integrand inside the path ordered exponential has to be single-valued. Therefore, the inclusion of the nonzero twist in the trace forces us to work with the twisted affine Lie algebra. It is easy to see that this twisted affine algebra is precisely the one we get by acting with U α + defined in (4.6b) on the untwisted affine algebra sl 2 .
It was proposed 2 in [19] that the expectation values in a state | of such a Kondo line defect will coincide with the transport coefficients of the Miura λ-oper, where the quadratic differential is is constructed from solutions of the Bethe equations for the state | . The corresponding Schrödinger equation reads The UV perturbative expansion is available whenever the Stokes data for the Schrödinger equation becomes close to the Stokes data for the simpler equation In order to study the UV asymptotics (0 < g i 1), we rewrite this as where e 2θ = 1 is identified with an RG scale. When t(x) = 0, the Miura λ-oper describes the vacuum state whereas non-trivial t(x) = a + (x) 2 + ∂ x a + (x) that we build from the solution to the Bethe equations (3.34) describes a more general state | .
Let ψ(x; θ) be the unique small solution, whose precise meaning is defined in the next section, that decays exponentially fast along the ray of large real positive x + θ. The normalization is chosen to match asymptotically the WKB series in (6.5). By the ODE/IM correspondence, we identify As the transfer function T n; (θ) is independent of x, we can evaluate the Wronskian in a convenient region where the explicit form of ψ is accessible. For the case of SU(2) k chiral WZW model, this region is [19] for details). This allows for a systematic expansion in powers of the g i .
In the same vein as in [23,65], we can also define Q-functions, essentially as the coefficients of ψ(x) in an expansion at large negative x. If the twist α + in a + (x) is zero,

JHEP01(2022)175
then the construction of Q-functions is given in [19] and reviewed in appendix A. As a result, (5.8) becomes expressible as a quantum Wronskian a form that is familiar from the integrability literature. This turns out to be a useful tool in the perturbative calculations as well. We can find the expression for Q andQ to sufficient order in g by comparing to a direct perturbative evaluation of (5.7). For the ground state in a generic spin l module, the above claim has been verified in [19]. In this paper we will present a few examples of the claim for excited states in appendix A.

The coset scaling limit
We can scale the variable x in the λ-oper to get a slightly different parameterization fixed will bring this to a λ-oper which is naturally associated to integrable lines in a coset model.
Physically, we are sending the g i to infinity while keeping their ratios fixed and adjusting the RG flow scale. We expect the Kondo lines RG flow in that limit to admit an intermediate regime where the line defects become effectively transparent to the overall WZW currents, so that they can be identified with defect lines in a coset model. It would be nice to explore this limit more carefully.

The IR expansion
The Kondo line defects L j [θ], where j is the spin that labels the SU(2) representation and θ is the spectral parameter, are asymptotically free line defects defined in a product of WZW models. They have a non-trivial, possibly non-perturbative RG flow which can be explored by looking at their action on the circle Hilbert space with the help of the ODE/IM correspondence. In the UV, the action is given perturbatively by the corresponding operator T j , defined in (5.1).
In the IR, the defects will flow to conformally-invariant defects. Because of their chiral nature, in the IR they will commute with both holomorphic and anti-holomorphic stress tensor and define topological line defects. A rich CFT such as the product of WZW models can have a very large variety of topological line defects, which commute with the stress tensor but with little else.
The ODE/IM correspondence, though, gives immediate evidence that the IR limit of Kondo defects should be more special than that, and commute with all the Kac-Moody currents. Indeed, we will see that the far IR limit of the ODE/IM solution is controlled by a WKB leading answer which depends very little on the details of t(z), up to the choice JHEP01(2022)175 of l i . In particular, they are blind to the details of the Bethe roots, which control which current descendants one is taking expectation values on.
Line defects which commute with the whole current algebra of the product of WZW models are referred to as Verlinde line operators. They are labeled by the Kac labels, i.e. same as current algebra primary operator. They are products of individual Verlinde lines in each WZW factor.
Denoted by L j , with j = 0, 1 2 , 1, . . . , k 2 , their expectation value in the vacuum state, often called the quantum dimension, is given by In the ground state of spin l, or in any descendant of that, their expectation value is .
If we go to the IR, but not to the infinitely far IR, the Kondo line defects will be described as IR free deformations of some sums of products of Verlinde lines. The deformation can involve any SU(2)-invariant local operators supported on the Verlinde lines. For generic Verlinde lines, there are many such operators, looking like descendants of chiral primary fields of various spins. We will see that the subleading WKB corrections do generically involve fractional powers of the scale e θ which can be explained by the conformal dimension of these operators.
As we mentioned in section 2, something special happens when the far IR defect is the identity line, or some other Verlinde line which does not support non-trivial chiral primaries. In such a situation, the IR free deformation must involve the integral along the line of SU(2)-invariant bulk chiral operators, starting with the stress tensor. The expectation values of these deformed identity lines are simply the exponential of the zero-modes of these bulk operators, which behave as local Hamiltonians for the affine Gaudin model.
We will see below that the ODE/IM correspondence predicts such IR destiny for line defects associated with pairs of Stokes sectors which are joined by a generic WKB line. 3 The number of such pairs is precisely the same as the number of zeroes for ϕ(z), as expected from the classical affine Gaudin model. The WKB expansion of these reproduce the expectation values of the quantum local Hamiltonians.

Vacuum state t(x) = 0
Let us first focus on the case of single SU(2) and t(x) = 0. The Schrödinger equation takes the form In this section, we are interested in the IR limit λ −1 = e θ → ∞, where Voros/GMN-style WKB analysis is applicable. The analysis is essentially the same as in [19], except that we are also interested in the sub-leading terms in the λ expansion.

JHEP01(2022)175
We will briefly review the analysis and leave details in appendix B. One starts by reading off the quadratic differential P (x)dx 2 = e 2x (1 + gx) k dx 2 , which has a zero of order k at x 0 = −1/g and an exponential singularity at infinity. For any angle ϑ ∈ R/2πZ, ϑ-WKB lines are curves in the complex plane where where ∂ t is the tangent vector of the curve. One such line passes through any point in the x plane. Generic WKB lines go to positive infinity in both directions, joining two Stokes sectors there. Special WKB lines hit a zero of P (x) such as x = − 1 g or flow to negative infinity. 4 The union of special WKB lines is called WKB diagram/spectral network. ϑ is chosen such that e θ lies in the half plane centered on e iϑ , where the WKB approximation gives the correct e −θ → 0 asymptotics. The structure of the spectral network governs which solutions of the Schrödinger equation have a specific WKB asymptotic expansion.
In our current example, the structure of the WKB diagram is shown in figure 2. Special WKB lines go towards positive infinity along the positive real x + θ + iπn, n ∈ Z direction. There are k + 2 of them that are connected to the order k zero x 0 = − 1 g . The remaining special WKB lines go towards negative infinity with imaginary part shifted by ± π 2 k. Next, one needs to find a set of solutions, referred to as small solutions, which decrease exponentially fast along the Stokes lines towards positive infinity (and thus along WKB lines asymptoting to them). In particular, we define ψ 0 (x) to be the small solution that decreases fast along the line of large real positive x + θ and agrees with the WKB asymptotics along this line where S asym (x, e θ ) is the primitive of the WKB momentum, given by an asymptotic series in large x and small e −θ . Although the WKB momentum is uniquely defined, its primitive needs a choice of integration constant. We choose the leading term to be Different choices clearly lead to different normalization of the Wronskians, namely T functions (5.8). This choice has the nice property later on that the exponent of (6.8) is zero at the leading order. In a practical calculation involving subleading terms, one also needs to make a choice for every order in e −θ . Then for all n ∈ Z, we have a small solution ψ n (x) ≡ ψ 0 (x; θ + iπn) along the large positive real x + θ + iπn direction.

JHEP01(2022)175
Next, we want to use the WKB network to evaluate the WKB asymptotics of the Wronskians. In the standard Voros/GMN-style WKB analysis, one studies Wronskians between two of the small solutions joined by a generic WKB line. These Wronskians are controlled by the contour integral along the WKB line of the WKB one form whose leading term is P (x)dx. This collection of Wronskians is incomplete, though, unless all zeroes of P (x) are simple.
As has been developed in [19] and appendix B, WKB analysis can be generalized to study non-simple zeros. Roughly speaking, one also needs the information around the matching regions, which, in the current example, are the order k zero x 0 = − 1 g , and the large negative x. Correspondingly, one can derive WKB asymptotics for Wronskians between two of the small solutions joined by a special WKB line to the same zero, or to negative infinity.
Following from the WKB diagram shown in figure 2 let us suppose, for convenience, that the numbering of the special WKB lines that are connected to the zero at x = − 1 g is n 0 , n 0 + 1, n 0 + 2, . . . , n 0 + k + 1. The precise value of n 0 depends on the parity of k and Im θ, which are given in [19] and are not important to us. There are three different scenarios: • Wronskians between two of the k + 2 small solutions whose special WKB lines are connected to the zero, namely i(ψ n , ψ n ) for n 0 ≤ n < n ≤ n 0 + k + 1.
• Wronskians between small solutions whose special WKB lines are connected to large negative matching region (to be made precise below), namely i(ψ n , ψ n ) for n < n ≤ n 0 or n 0 + k + 1 ≤ n < n or n ≤ n 0 < n 0 + k + 1 ≤ n .
• The remaining ones can be related to the first two scenarios using Plücker formula.
In particular, to deal with the first case, it is important to study the local behavior around the zero x 0 = − 1 g . Locally around x 0 , a zero of order k, the stress tensor should take the form y k + . . . , with y being the coordinate in the local coordinate system. Indeed, one can always find the coordinate transformation x → y(x) such that stress tensor takes the form of y k + a k−2 γ k y k−2 + · · · + a j γ 2+j y j + · · · + a 0 γ 2 (6.7) where γ = e −θ 2 k+2 . Importantly there are k −1 coefficients a j = a j + . . . that are uniquely fixed in γ k+2 asymptotics. One can find a set of nice solutions A k;i (y) to this local problem and evaluate the Wronskian perturbatively. The general procedure to do this is described in appendix B. On the other hand, Wronskians between small solutions ψ n (x) are equal to the Wronskians between the corresponding local solutions A k;i (y) with a careful treatment on the normalization of the solutions. We will only quote the result here, leaving the details in appendix B, whose leading term is given by the quantum dimension defined as .
Subleading terms are computable order by order in γ. See appendix B for the general prescription. There are two exceptions k = 1, 2 where we can calculate Wronskians exactly. The important part is that the corrections to the Wronskians come in as integer power of γ but start from γ 2 order. In the second scenario, the special WKB lines 'meet' at the large negative x. It turns out, for some suitably chosen x −∞ , a shift of the coordinate x → δ = x − x −∞ will transform the quadratic differential into The details are given in appendix B. Here, what matters to us is that x −∞ has a large negative real part and in the IR limit θ → ∞ we have The coupling g eff (θ) is defined by the relation and goes to 0 in the IR limit θ → ∞. This is precisely the effective coupling for the infrared free line defect, whose physical meaning will be given below in section 6.3. For now, we only need the fact that g eff (θ) → 0 as θ → ∞. Therefore, we can study the Wronskians of solutions in g eff expansion.
In the leading order, the local solutions are given by Bessel functions. Therefore by means of Bessel function identities and with normalization factors carefully taken into account, the results are whenever n ≤ n 0 and n ≤ n 0 , whenever n ≥ n 0 + k + 1 and n ≥ n 0 + k + 1, and 15) whenever n ≤ n 0 and n ≥ n 0 + k + 1.
In the third scenario, we can just use Plücker formula to reduce to the previous two scenarios. Details can be found in [19].
Recall that the leading order of the second term agrees with the one from the UV expansion and intuitively just counts the number of spacing between different special WKB JHEP01(2022)175 lines at the left hand side of the special WKB diagram, see e.g. figure 2. The exponential factor is the non-perturbative ground state energy shift.
As we discussed in section 2 and at the beginning of this section, the vevs of local integrals of motion for the affine Gaudin model are given by the Wronskians which correspond to generic WKB lines. In the example at hand, there is only one such line depicted by the dotted burgundy line in figure 2, corresponding to L k 2 in the infrared. Indeed, the L k 2 line does not support nontrivial chiral WZW primaries. Since the corresponding Wronskian is controlled by the contour integral of the WKB momentum along the generic WKB line, it doesn't involve the local analysis around the zero or negative infinity, hence it is simply organized by odd powers of e −θ .

t(x) = 0
When t(x) = a + (x) 2 + ∂a + (x) = 0, the evaluation of the Wronskians via WKB analysis is basically the same as the previous section except for the following modifications.
In the first scenario of the previous section, namely, for the Wronskians of two solutions whose special WKB lines are connected at the zero x = − 1 g , the local coordinate system in general has an additional piece, compared to (6.7) where −l is the residue of a + (x) at the zero. This will change the leading order of (6.8) to be The nonzero regular part of a + (x) has smaller impact. It will change the coefficients a j , the details of the map x → y(x), and therefore the higher order corrections. But importantly, the corrections are still organized by integer powers of γ.
In the second scenario, where two special WKB lines are connected at the negative infinity, we can again go to the coordinate in δ = x − x −∞ where the quadratic differential reads We are then in a situation very similar to the UV expansion. Therefore, the higher order corrections of the Wronskians come in powers of g eff .

Physical interpretation
According to the ODE/IM correspondence (5.8), which we repeat here, the expectation value of the Kondo line operator in the state | is given by where n = 2j + 1. We showed in section 5 that the leading term in the UV expansion is given by the dimension n of the representation, T n; ∼ n + . . . .

JHEP01(2022)175
We now provide the physical implication of the IR expansion we evaluated using WKB analysis previously in this section. The IR expansion of T n; (θ) reviews an interesting infrared structure. The leading order has been demonstrated in [19]. We will review briefly now and explain how the structure of higher order corrections we obtained in this section fits in the paradigm.
Depending on the imaginary part of θ, and whether 0 ≤ 2j ≤ k or 2j > k, the RG flow takes the Kondo line operator L j [θ] to different IR line operators.
Firstly, if θ is real, or more precisely, valued in a strip around the real θ axis of width about 5 (n − k − 1)π, we have the physical RG flow: 6 Second, if we increase the imaginary part of θ either positively or negatively, there is an interesting sequence of transitions starting from | Im θ| ∼ (n−k−1)π 2 , the edge of the strip mentioned above. Every time | Im θ| increases by π, we trade one unit of spin for the topological defect with one unit of spin for the internal degrees of freedom. More precisely, After s decreases to zero, i.e. when | Im θ| is large enough, we will have the circular RG flows where L j flows to L IR j . Since e θ labels the RG scale and θ → ∞ is the deep infrared, we have just demonstrated that the leading term of T n; (θ) simply tells us which infrared line defects we flow to starting from L j defined in the UV. Correspondingly, the far IR destiny of UV line defects can be determined from some simple combinatorics from the topology of WKB network.
Subleading terms are obviously due to the deformation that brings us away from the deep IR. More precisely, we need to look at the RG flow a bit away from the deep IR. In the case of the RG flow that flows to the Kondo line defect L IR j , an infrared free line defect, the effective coupling g eff (θ) is negative and becomes smaller in the IR. Therefore, the corrections are expected to come in powers of g eff , which is small and negative.
On the other hand, topological lines in the infrared are not free and appear in the RG flow as strongly coupled infrared fixed points. 8 Nevertheless, we know a lot about the local operators that are supported on the lines. Among them, it is the least irrelevant operator that contributes the corrections the closest to the deep IR fixed point [12]. On a spin 0 < j < k 2 Verlinde line in chiral WZW, it is the WZW descendant of the spin 1 operator of dimension 2 k+2 , denoted by J a −1 φ a . It is a Virasoro primary of scaling dimension 1 + 2 k+2 and has zero one-point function. Then a simple dimensional analysis implies that there will be corrections in powers of γ = e −θ( 2 k+2 ) , starting from γ 2 . Another obvious candidate J a J a has dimension 2, thus contributing corrections in integer powers of = γ k+2 2 , which is 5 This is true up to ±π/2, depending on the parity of k and n. 6 Recall that e θ labels the RG scale, so the physical RG flow corresponds to real θ. 7 This terminology is based on the intuitive physical picture that Kondo defect disappears in the IR because magnetic impurity spin is screened by the bulk fermions. See e.g. [14] for more details. more irrelevant. Note that the spin 1 operator φ a is not supported on the Verlinde line L k 2 , obtained in the exact screening case. This is precisely what we found in the WKB analysis around (6.7) and (6.8).

Some comments about the semiclassical limit
The RG flow of the Kondo defects is often non-perturbative. An important exception occurs when the levels k i are large: in an appropriate RG scheme the couplings remain small all the way to the IR. It is useful to illustrate this in the N = 1 case.
We can start from the quadratic differential λ −2 e 2x x k dx 2 . The RG prescription we used before for the couplings sets

JHEP01(2022)175
so that the quadratic differential becomes e 2x− 2 g (gx) k dx 2 which can be mapped by a translation to e 2x (1 + gx) k dx 2 and treated perturbatively. In these conventions, g is small as λ 1 in the UV, but flows all the way to infinity in the IR as λ 1. When k 1, these RG conventions are not very good. For example, even if g ∼ k −1 is small, (1 + gx) k ∼ e gkx is not close to 1. Consider, instead, a different definition of g where we employ both a translation and a scale transformation of x to arrive to some e 2α(g)x (1 + gx) k dx 2 . We can map this back to e 2x− 2α(g) g g k x k α(g) −k−2 dx 2 and thus to If we select α(g) = 1 − k 2 g, then e 2α(g)x (1 + gx) k dx 2 ∼ e 2x dx 2 up to corrections of order g 2 k, which are small even if g ∼ k −1 .
This seems a more reliable way to deal with the perturbative RG flow. Now we have which flows from g = 0 to g = 2 k from the UV to the IR, so that the coupling is perturbatively small all the way.
This RG scheme also seems appropriate to make contact with the classical affine Gaudin model in a k → ∞ semiclassical limit. In general, define This definition is inspired by the classical affine Gaudin Lax matrix (2.1). Then e 2x ϕ(z) can be mapped to and then to e −2z ϕ(z) 2 so that we have a z RG flow controlled by and fixed points at z ∼ z i where g i is finite and small, while all other g j vanish.
8 Generalizations: sl 3 The discussion in section 3 of the affine sl 2 opers, the Bethe equations and WKB solutions can be easily generalized to higher rank Lie algebras. We demonstrate the case of sl 3 in this section. We leave a proper discussion of the corresponding Kondo defects and ODE/IM solutions to future work.

Basic definitions
An sl 3 oper is a complexified third order differential operator with a natural transformation law under a change of coordinate We will work with sl 3 opers for which both t 2 (x) and t 3 (x) are rational functions. The data of an sl 3 oper is equivalent to that of a flat connection More generally, an sl 3 oper can be described as a flat connection of the form An sl 3 λ-oper is a complexified third order differential operator with a particular dependence on the auxiliary complex parameter λ of the form where t 2 (x) and t 3 (x) are rational functions. Equivalently, we can describe this as a flat connection More generally, an sl 3 λ-oper is defined as a connection of the form where a(x),ã(x), b(x),b(x) and c(x) are rational functions, modulo gauge transformations by upper-triangular matrices of the form where v(x),ṽ(x) and w(x) are rational functions. Every sl 3 λ-oper has a unique canonical form as in (8.7). If we conjugate the connection (8.8) by cyclic permutation matrices then we obtain two alternative formulations of the differential operator (8.6), leading to two alternative (but equivalent) formulations of sl 3 λ-opers. Specifically, an sl 3 λ-oper can equally be described as a flat connection of the form Equivalently, an sl 3 λ-oper can also be described as a connection of the form The unique canonical form of an sl 3 λ-oper in the second description (8.10) is given by and that of an sl 3 λ-oper in the third description (8.12) reads (8.15) JHEP01(2022)175

Miura λ-opers and singularities of trivial monodromy
A Miura sl 3 oper is a connection of the form with a(x) andã(x) rational. Since it is of the general form in (8.4) it defines an sl 3 oper which corresponds to the differential operator There are two types of apparent singularities, corresponding to the two nodes of the Dynkin diagram of sl 3 . These can be points w where so that, in particular, the constant term of a(x) + 2ã(x) vanishes.
If at a singularity z the Miura sl 3 oper is of the form for some non-negative integers n 1 and n 2 , then z is a regular singularity of the sl 3 oper of trivial monodromy. Indeed, one can bring the Miura sl 3 oper to the form x−z , which are regular at z. A Miura sl 3 λ-oper is a connection of the form where a 1 (x) andã 1 (x) are rational functions. This is of the general form (8.8) and so a Miura sl 3 λ-oper defines an sl 3 λ-oper with We refer to this as the sl 3 λ-oper underlying (8.21). It can be described as a third order differential operator of the form

JHEP01(2022)175
There are two other gauge equivalent ways of presenting the same Miura sl 3 λ-oper as in (8.21), namely 3P (x) . The Miura sl 3 λ-oper (8.24) is of the particular form (8.10) so it defines a second sl 3 λ-oper. Likewise, the Miura sl 3 λ-oper (8.25) is of the form (8.12) and thus it also defines a third sl 3 λ-oper. Crucially, all three sl 3 λ-opers share the same monodromy data since they are gauge equivalent. This triality generalises the duality of sl 2 λ-opers associated with a given Miura sl 2 λ-oper discussed in section 3.4.
There are three types of apparent singularities, corresponding to the three nodes of the Dynkin diagram of sl 3 . In particular, we can have the same types of singularities as for a Miura sl 3 oper in (8.17), namely points w where, cf. (8.17), so that a 1 (x) −ã 1 (x) has vanishing constant term, or points w where, cf. (8.18), so that a 1 (x) + 2ã 1 (x) has no constant term. Both of these singularities are absent from the sl 3 λ-oper underlying (8.21). The third type of apparent singularity is at points w where If P (x) has a zero of order k at a singularity z of the Miura sl 3 λ-oper with a 1 (x) = − 1 3

JHEP01(2022)175
then provided n 1 , n 2 ≥ 0 and n 1 + n 2 ≤ k, the underlying sl 3 λ-oper has trivial monodromy. Indeed, one can bring (8.21) to the form x−z , which are clearly regular at z. The behaviour of the second Miura sl 3 λ-oper (8.24) at z is given by while the third Miura sl 3 λ-oper (8.25) behaves as

λ-opers with singularities of trivial monodromy and affine Bethe equations
A Miura sl 3 oper on C with a rank 1 irregular singularity at infinity and whose other singularities are all regular with trivial monodromy is of the form where the apparent singularities w i and w i satisfy the Bethe equations − a n 1,a We are interested in the case of a Miura sl 3 λ-oper with a rank 1 irregular singularity at infinity and whose other singularities are all regular with trivial monodromy. This can be written as With P (x) = e (α 1 +α 2 +α 3 )x a (x − z a ) ka , the second Miura sl 3 λ-oper then reads

JHEP01(2022)175
The condition that w i and w i are apparent singularities for the first Miura sl 3 λ-oper (8.37) and that w i are apparent singularities for the second Miura sl 3 λ-oper (8.38) leads to the Bethe equations − a n 1,a

WKB expansion and quasi-canonical form
The three Miura sl 3 λ-opers (8.21), (8.24) and (8.25) are locally gauge equivalent to a connection of the more symmetric form . We refer to (8.42) as a Miura sl 3 oper. The underlying sl 3 oper, or affine sl 3 oper, is then defined as its equivalence class under gauge transformations by matrices of the form where the various functions have the following formal power series expansions with u n (x),ũ n (x), v ± n (x),ṽ ± n (x) and w ± n (x) rational functions. As we will show below, the sl 3 oper controls the WKB asymptotics of the λ-oper (8.6) underlying the Miura sl 3 λ-oper (8.21); see figure 3.
Recall first that an sl 3 oper can be brought to the quasi-canonical form [46]

JHEP01(2022)175
Miura  where the coefficients are given by the formal Laurent series In the case of the sl 3 oper underlying (8.42), the first few orders explicitly read where t 2 (x) and t 3 (x) are given by (8.22). Just as in the sl 2 case considered in section 3.10, the quasi-canonical form (8.44) of an affine sl 3 oper is not unique. Indeed, it is preserved by residual gauge transformations of the form (8.43) with u(x; λ) =ũ(x; λ) = 0 and v ± (x; λ) =ṽ ± (x; λ) = w ∓ (x; λ), the effect of which is to transform the quasi-canonical form as We look for flat sections of the sl 3 λ-oper (8.6), i.e. solutions of the third order differential equation in the form of the WKB ansatz

JHEP01(2022)175
where the normalization factor, fixed by requiring the Wronskian of the three solutions ψ 1 , ψ 2 and ψ 3 to be 1, is given by Working perturbatively in λ we find WKB momenta of the form The first few coefficients of the expansion (8.49a) are given explicitly by and those of the expansion (8.49b) read 9P (x) 4 3 , (8.51a) From this we find that the pair of WKB momenta p 1 (x; λ) and p 2 (x; λ) are related to the coefficients of the quasi-canonical form (8.44) by for some functions f 1 (x; λ) and f 2 (x; λ).

JHEP01(2022)175 9 Future directions and open problems
1. Higher rank. It would be natural to extend our work to other Lie algebras [66][67][68]. The definition of the Kondo defects is valid for any Lie algebra, with an important caveat: the matrices t a do not have to be generators of the Lie algebra, they only need to transform in the adjoint representation of the global symmetry group. The RG flow will thus involve extra couplings, controlling the specific choice of t a . 4d Chern-Simons gauge theory predicts integrability for choices of couplings related to representations of the Yangian. It would be nice to understand how these structures manifest themselves in the affine Gaudin description. The definition of affine opers with singularities of trivial monodromy should be possible for general gauge groups and straighforward in type A. The precise correspondence between the Stokes data and the UV labels of Kondo defects is less obvious. The IR WKB analysis is still possible, but will likely require some more refined technology such as spectral networks [38].

Non-integral levels.
In sections 5 and 6, we made the assumption that all the level k i in are non-negative integers, which corresponds to the WZW model for the product group i SU(2) k i . We can generalize to Kac-Moody algebras via replacing k i by κ i ∈ C. This forces us to study opers on a logarithmic covering space of the complex plane. (It reduces to a finite covering when κ i is rational.) Consequently, we have more Wronskians of small solutions to study, both between the small solutions on the same sheet and across different sheets. It would be interesting to work out some examples and understand their relationship with the expectation values of the local integrals of motions in the affine Gaudin model, which are conjecturally given by the integrals of the WKB momentum [46].

Bulk deformation. Throughout this article, we considered Kondo line defects in CFTs in
the bulk. It is well-known, however, that the bulk theory can be deformed in such a way that the integrable structure remains. Examples of such deformations are given by J aJ a in the WZW model and Φ 1,3 in minimal models. The transfer matricesT j [θ] can be naturally extended as well in such a way that commutativity and the fusion rules are still satisfied. Furthermore, there is evidence [69][70][71] that the ODE/IM correspondence can be generalized as well to the so-called 'massive ODE/IM correspondence'. It is natural to expect that the P -sinh-Gordon equation with an appropriate generalization of the potential P (x) = e 2x i (x − x i ) k i should correspond to the integrable Gross-Neveu model.

Coset limits and other models.
We only sketched the limiting procedure e 2x i (x−x i ) k i → i (x − x i ) k i mapping the Kondo problems to integrable defects in coset models. It would be nice to develop the relation further. Other choices of P (x) should be relevant for integrable line defects in other 2d CFTs. A dictionary may be developed along the lines of [72] from a 4d Chern-Simons gauge theory perspective, or equivalently [73], along the lines of [74] from the point of view of affine Gaudin models.

Non-isotropic generalization.
Here we only considered the case of Kondo problems which preserve the global SU(2) symmetry of the WZW models. Anisotropic Kondo problems are also very interesting and often still integrable. In [62], an ODE/IM solution was proposed for the ground states of the anisotropic defects in a single WZW model. It would be very interesting to find a full solution to the problem.

A Details of UV expansion
In this appendix, we verify the claim (5.8) for a few examples explicitly by working perturbatively in g to O(g 4 ). Next, based on these examples, we summarize a general recipe in the case of zero twist α + = 0. We then end this section with some remarks and checks with nonzero twist.

A.1 Examples of UV perturbative matching
Let us perform a perturbative UV analysis of the proposed excited state Schrödinger equations and compare them to SU(2) k WZW line defects evaluated between certain excited states. We use the Chevalley basis rather than the orthonormal basis used in [19]. To establish notation, we write the change of basis explicitly: We will consider the following states as examples: From the discussion in subsection 4.3.1, for generic values of k, they are described by Miura λ-opers with a (1) respectively. Note, however, that the state J + −2 |0, −1 is a special case of (A.4) since, as we conjectured in 4.3.1, it corresponds to the generalized Miura λ-oper

JHEP01(2022)175
The left hand side of the claim (5.8) is an obvious, alebit tedious, task, namely to compute the g expansion of the expectation value of the operatorT n (θ) in a certain state. This operator was calculated in the appendix E of [19]. Note that if we normalize the expectation values of the line defect as 9 then the ordinary form of the Hirota bilinear relations holds. 10 Finding the g expansion of the Wronskian on the right hand side of (5.8) is not obvious. In [19], we proposed an approach based on a combination of various limits. The basic idea consists of two steps: (1) analyze the solution of the Schrödinger equation in the regime 1/g −x and −x 0 separately and (2) match the solutions in the intermediate region . We can recycle many of the result from [19] as modifications of the Schrödinger equation by the potential terms do not contribute to the perturbative expansion of the differential equation at O(1) and O(g). Hence, we may use the direct evaluation result of the wavefunction in [19] given by at large negative x. Also, we introduce the effective coupling g eff (θ) defined by [19] g eff (θ) k/2 e −1/g eff (θ) = g k/2 e −1/g e θ (A.10) with the g expansion In this appendix only, we write the excited expectation values with angled brackets to emphasize that Tn(θ) r are quantities which are directly evaluated from the defect operator computed in [19]. The quantities we write as Tn(θ) in the perturbative analysis below correspond to the Wronskian result we obtain from the Schrödinger equation. 10 One could choose not to use the normalization 1 rk in (A.7), but the Hirota relations would then have (rk) 2 in place of 1 and the following perturbative analyses acquire extra normalization factors.

JHEP01(2022)175
A.1.1 J + −1 |0, k , one w , no w For the case with one w and no w, the excited state equation is We propose that this equation encodes the excited line defect expectation value T n (θ) 1 . The perturbative analysis as an expansion in g is done after shifting x → x + 1/g. The asymptotics of ψ can be determined by carefully analyzing the solutions of the Schrödinger equation in the regime 1/g −x 0. Doing so, the wavefunction can be parametrized as Before the shift x → x + 1/g, the dependence of ψ on g and θ combine into g eff (θ) such that ψ becomes a function of x and g eff only. The explicit g dependence in the parametrization of ψ above comes purely from the shift x → x + 1/g. As the general analysis in the next subsection will suggest, the difference #w − #w is responsible for the various factors present above. The overall coefficient 1/3 in the asymptotics of ψ arises via the combination (2(#w − #w) + 1) −1 , as does the overall power of g ±(#w −#w) eff in the expression of Q,Q. As usual, the T-function can be expressed as the quantum Wronskian The coefficients q i can be expressed in terms ofq i by imposing the condition T 1 = 1. The explicit expression of T n (θ) up to O(g 4 ) requires the knowledge ofq 0 andq 1 . Comparing the parameterization of ψ in terms of Q,Q with the direct perturbative evaluation up to O(g), we obtaiñ The full expression for T n (θ) to O(g 4 ) is then T n (θ) = n − g 2 1 3 π 2 n n 2 − 1 + g 3 1 12 π 2 n(n 2 − 1)(7k − 8θ − 8γ + log 256) − g 4 1 288 π 2 n(n 2 − 1)[195k 2 − 120k(5θ + 1 − 5 log 2) − 24γ(25k − 24θ + log 16777216) + 8π 2 (10 − 3n 2 ) + 288(θ − log 2) 2 + 288γ 2 ] + O(g 5 ).

A.2 A general perturbative prescription
One can easily see the pattern of the steps in the previous example and might wonder if it is possible to perform the perturbative analysis in a uniform way without specifying particular values for w i and w i , i.e. without solving the Bethe equations explicitly. Here in this section, we provide such a recipe for zero spin and zero twist.
In considering the Schrödinger equation in the regime 1/g −x 0, one first solves the equation with just the potential terms and then expands in g until one has O(x) contributions. That is, for (#w , #w) = (p, q) we take (after shifting x → x + 1/g) the g expansion of the solution to Notice that one exact solution of such an equation is Expanding up to the order shown is sufficient for p > q, but one needs to consider an O(g 2 ) term for p = q. Here, we concern ourselves with just the p > q cases as the p = q case is treated in a similar but more involved way.

JHEP01(2022)175
The solution above is that proportional toQ of the previous subsection. An exact second solution turns out to be more elusive, but this can be overcome by the fact that we only need the perturbative form of second solution. The second solution can then be determined by imposing that 24) i.e. the shifted Wronskian of the wavefunction, is normalized as Doing so, one finds that the parametrization of the wavefunction for the cases p > q becomes We see thatq 0 andq 1 , determined by comparison to the direct perturbative g expansion of the Schrödinger equation, will only depend on the difference of sums of w 's and w's. The analysis in this subsection can be generalized to other spin modules and those with twists.

A.3 Nonzero twist α + = 0
Mimicking the case of α + = 0, we proceed with the UV expansion as follows. We are interested in the ODE Consider the matching region 1 g −x 0. The first inequality means we are reduced to The solution is given by

B.1 General remarks and organizations
Given a sl 2 λ-oper, the goal of this appendix is to develop techniques to evaluate the Stokes data, namely the Wronskians between certain solutions (referred to as small solutions

JHEP01(2022)175
defined below) to the corresponding Schrödinger equation. Our formalism is based on the WKB analysis, which culminates in the Voros analysis [28][29][30][31][32][33][34][35][36][37][38]. We refer readers to the appendix of [19] and references therein for a review of WKB analysis and other related recent progress. Contrary to the common wisdom, we will find that in the general situation including the examples in this article, the standard WKB analysis will not provide a complete collection of Stokes data, which we briefly review below. The WKB analysis is concerned with the holomorphic differential equation defined on a Riemann surface by using the oper coordinate transformation between patches Here we start with brief definitions of two main players, Stokes diagrams and small solutions.
Other relevant notions will be explained wherever needed.
T (x) and P (x)dx 2 are referred to as the stress tensor and the quadratic differential respectively. We will briefly review the analysis. For any angle ϑ ∈ R/2πZ, ϑ-WKB lines are curves in the complex plane where where ∂ t is the tangent vector of the curve. One such line passes through any point in the x plane. a WKB curve is special if it ends on a zero of P or on an asymptotic region of exponentially fast decrease for P . The union of all special WKB lines is called a WKB diagram/spectral network. See footnote 4. For each singularity of P (x), we can define a set of small solutions. In particular, for each special WKB line coming into the singularity, we define a small solution to be the unique solution (up to normalization) that decays exponentially fast approaching the singularity along the corresponding special WKB line. We will choose the normalization such that asymptotically near the singularity, it matches with the WKB solutions ψ WKB ± (x). The definition of the WKB solutions and the associated WKB coordinate system will be given in section B.2.1.
We would like to find the Wronskians between small solutions. The standard GMN/Voros-style WKB analysis is interested in contour integrals along the generic WKB lines connecting every pairs of Stokes sectors, while, importantly, staying away from the zeros of the quadratic differential P (x). They are indeed all we need to provide a complete collection of Wronskians when all zeroes of the quadratic differential P (x) are simple. However, more generally, the collection is not complete and we also need information from the local analysis of the matching region.
An example where we need a generalized WKB analysis would be a polynomial oper with non-simple zeros. And the zeros are the matching region we need to understand. More precisely, for each zero of P (x), we can define a local coordinate system and nice local JHEP01(2022)175 solutions therein. Let's consider first t(x) = 0. For each zero x 0 of order n, we want to find a coordinate systemỹ(x) such that the oper of interest whereỹ(x 0 ) = 0 + O( 2 ). We also require the (possibly nonzero) subleading terms inỹ(x) are regular at x 0 .
Generically one can't findỹ(x) such that the stress tensor takes exactly the formỹ n 2 . Instead, it turns out the best one can do is y n 2 + a n−2ỹ n−2 + · · · + a 1ỹ + a 0 .
The constant coefficients can be determined order by order in , In particular when n = 1, namely, around a simple zero, there always exists a coordinate system such that the stress tensor takes the formỹ 2 . We can also write y = − 2 n+2ỹ , in which the stress tensor takes the form y n + a n−2 2n n+2 y n−2 + · · · + a j 2j+4 n+2 y j + · · · + a 0 4 n+2 . (B.8) We will discuss how to find such a local coordinate system and nice local solutions defined in there, as well as the cases with t(x) = 0, in subsection B.2.2.
Another example of the matching region is the negative infinity of the oper with exponential potential. The details on how to deal with such matching region is given in subsection B.7.3.
We will first describe three different coordinate systems in detail and the transformation between them in section B.2. Next in section B.3, we will explain why the WKB analysis away from the zeros is not enough and what information from the local analysis is crucial to give a complete collection of Wronskians/Stokes data. We then give a general recipe for evaluating the Wronskians by incorporating the missing local information while deferring the local perturbative analysis of extracting the local information to section B.5.
In order to verify our proposed recipe, we carry out numerical computations. The method of numerical implementation is given in section B.4.
Finally, we apply the generalized WKB analysis to the nice examples at hand including polynomial oper with non-simple zeros in section B.6 and the exponential case that describes the chiral WZW model in section B.7.

B.2 Coordinate systems
We would like to define interesting local coordinate systems with good → 0 asymptotics. There are a few useful coordinate systems we will use frequently in this paper: • the original one, denoted as x, where T = 1 2 P (x) + t(x), • the local coordinate around a zero, y orỹ, where T = y n + . . . orỹ n 2 + . . . , • the WKB coordinate s ab where T = 1 4 .

B.2.1 WKB coordinate systems
Near each (Stokes sector of a) singularity a of the quadratic differential P (x)dx 2 we can find a solution ψ WKB a which decreases exponentially fast approaching the singularity. It is given as a specific asymptotic expansion near the singularity, where s asy a (x) is a primitive of twice the WKB momentum p asy a (x) = 1 2 ∂ x s asy a (x), which satisfies the differential equation The one form p asy a (x)dx is also referred to as the WKB one form. If we write p asy a (x) in asymptotics as p asy a (x; ) = p asy −1 (x) + p asy 1 (x) + 3 p asy 3 (x) + · · · (B.11) the equation (B.10) can then be written as a recursive equation for p n (x). The leading term is given by p asy −1 (x) ∼ P (x). For a polynomial singularity at x = ∞, the asymptotic expansion s asy a (x) involves increasingly negative powers of x, with coefficients which are Laurent polynomials in . In order to fix the normalization at x = ∞ we only need to worry about powers of x greater than −1, so the last two terms are unimportant. For a singularity of odd degree, the expansion involves fractional powers of x and thus s asy a (x) can be chosen unambiguously to have no constant term. For a singularity of even degree, there will be a log x term and we will have to make some choice. Of course, sometimes there are some natural nonzero choices of the constant as well. See, for example, section B.6.1. For more generic cases, we will fix the constant terms of s asy a (x) on a case-by-case basis. Once we fix the constant term in s asy a (x), the normalization of the WKB solutions ψ WKB a is then fixed. This then further fixes the normalization of the small solutions, which, as we defined in section B.1, are normalized to match the WKB solutions asymptotically near the singularity.

JHEP01(2022)175
Although we defined s asy a (x) only as an asymptotic series, it is easy to produce actual functions which have such an asymptotic expansion: given any other solution 11 ψ * which is not proportional to ψ a , we could define Redefinitions of ψ * by multiples of ψ a would just shift e sa, * (x) by an x-independent constant.
For generic x, we sit along a generic WKB line associated to a pair of small solutions ψ a (x) and ψ b (x), with some normalization determined at the corresponding directions at infinity or singularities.
These solutions have good → 0 asymptotics in a certain sector of width π in the plane. So does their ratio, which defines an useful local coordinate Then so that the solutions ψ a and ψ b map to 1 and (ψ a , ψ b )z ab in the z ab coordinate. Or equivalently, the stress tensor in z ab coordinate is T (z ab ) = 0. Also, and thus the solutions ψ a and ψ b map to e − 1 2 s ab (x) and (ψ a , ψ b )e 1 2 s ab (x) in the s ab coordinate. Equivalently, the stress tensor in s ab coordinate is T (s ab ) = 1 4 s ab (x) can be determined in the following way. Similar to (B.10), the WKB momentum p ab (x) ≡ 1 2 ∂ x s ab satisfies the differential equation which can be solved recursively to find the WKB expansion of p(x) away from zeroes of P (x): It is easy to see that p n>0 (x) will involve increasingly negative powers of p −1 (x) = P (x). Therefore the zeros of P (x) remain the places where the WKB approximation breaks down, regardless of t(x).
In order to compute s ab from p ab , we need to fix the integration constants. This is done by expanding p(x) near the singularity and comparing term-by-term with s asy a (x). 11 Unless we state explicitly otherwise, by a solution, we mean an actual function, as opposed to a WKB solution, which is an asymptotic expansion.

JHEP01(2022)175
It is easy to express s ab as a regularized contour integral. We can write and sendx towards the singularity x a : Another useful observation is that and p ab + p ba = 0, so that with the integral taken along a path equivalent to the WKB line between x a and x b . The overall sign can be determined by computing the Wronskian explicitly:

B.2.2 Local coordinate system around zeros
Here in this section, we explain how one can find explicitly the coordinate system in which 1 2 P (x) around a zero of order n takes the form of y n 2 + a n−2ỹ n−2 + · · · + a 1ỹ + a 0 (B.23) where a i = a In this section, we will mostly useỹ so that only integer power of will appear. When we have nontrivial t(x) = a(x) 2 + ∂ x a(x), and a(x) has residue −l at x 0 , y n 2 + a n−2ỹ n−2 + · · · + a 1ỹ + a 0 + l(l + 1) y 2 (B. 25) or y n + a n−2 2n n+2 y n−2 + · · · + a j 2j+4 n+2 y j + · · · + a 0 4 n+2 + l(l + 1) In particular if t(x) is nontrivial but regular at the zero of interest x 0 , we should find the stress tensor takes the local form (B.6) and (B.8).
Below we will see thatỹ(x) can be uniquely determined with no ambiguity and that it is unavoidable to have the coefficients a i by providing two different ways of finding such y(x). By doing so we will explain why the coefficients a m are unavoidable. (1) In comparing the local coordinate and other coordinates, we need to have some parameters to adjust so that the cross ratios defined from the local solutions coincide with the same cross-ratio built from ψ n (x). (2) a m need to be specific values to make the coordinate transformatioñ y(x) non-divergent at the interested zero.

JHEP01(2022)175
Local coordinate transformation. We call this a local map since it is only valid in the neighborhood of the zero Suppose we are interested in the zero of order n, then apparently one can always do a shift in the coordinate such that where t 1 = t 2 = · · · = t n−1 = 0 and x = 0 is the order n zero of interest. By solving the following equation order by order in and x, we can determine the coordinate transformationỹ(x) and more importantly {a m } explicitly as functions of {t n , t n+1 , . . . , t N }.

Non-local coordinate transformation.
To state again, we want to find a coordinate transformationỹ(x) with respect to a zero x 0 of order n satisfying ∂ỹ ∂x . This equation can be solved recursively order by order in , and at each order we have a first order differential equation. The leading order equationỹ 2 0ỹ n (B.31) The integration constant is fixed by choosing the integration starting point from the interested zero x 0 of order n, such that locally around x 0 ,ỹ 0 (x) start from linear order in (x − x 0 ) and vanishes at x 0 . At the order of −1 , we have a homogeneous differential equation nỹ 1 y 0 + 2ỹ 0ỹ 1 = 0, which is solved byỹ Sinceỹ 0 (x) vanishes at the zero x 0 and we requireỹ(x) to be regular there, the only choice is to choose c 1 = 0, henceỹ 1 (x) = 0. The same is true for ever odd order in . At the order of 0 , we end up with the an inhomogeneous first order differential equation forỹ 2 (x). We fix the homogeneous part 1 y 0 (x) n/2 of the solution by choosing the lower limit of the integral at x 0 . This rendersỹ 2 (x) regular, and in general nonzero, at x 0 . For example Furthermore, one can also understand the necessity of the coefficients a i in (B.24) from the fact that we need them in order to haveỹ 2 (x) non-divergent at x 0 .

B.2.3 Relating two coordinate systems
Near a simple zero. There are three special WKB lines emanating from a simple zero. Therefore, near a simple zero, we can access three small solutions ψ a , ψ b , ψ c , decaying exponentially along the three special WKB lines respectively. Apparently there must be a linear relation among them, sometimes referred to as Plücker relations.
On the other hand, there exists a local coordinate system y associated to this simple zero, in which the stress tensor reads T (y) = y. The coordinate y satisfies the differential equation and can be expanded as → 0 to be y(x) = − 2 3 (ỹ 0 + 2ỹ 2 + 4ỹ 4 + . . . ). Then the nice local solutions in this local coordinate system are given by (B.86), In particular, we have This obviously satisfies also which gives us the relation between small solutions ψ • (x) and the local solutions Ai • (y(x)). For example, We can now relate y to s ab and to the other WKB coordinates nearby: The coordinate s ab (x) depend on the normalization of the local solutions Ai a (y) and small solutions ψ a (x). One can define an alternative local WKB coordinate, i.e. an alternative primitive of p ab (x), that is independent of the normalization of the solutions, given as follows

JHEP01(2022)175
If we expand the right hand side − Ai 1 (y(x)) Ai 0 (y(x)) in an asymptotic expansion at large y(x) using We can think about this as a regularized version of the integral of 2p ab from the zero to x.
More precisely, as all the ingredients of the definition above have good WKB asymptotics, we will show below that the expansion exactly coincide with the contour integral of 1 2 p ab from x to x along a path γ abc which winds around the zero while keeping away from it: namely, up to some multiple of i π 2 we have We should really keep track of factors of ±i in front of the exponents: We now show explicitly why (B.49) is true at the first two orders. Recall that the WKB momentum is The leading order is easy since P (x) is integrable at the zero and we can just pinch the contour to the zero The order is not as obvious

JHEP01(2022)175
is generically nonzero and regular as x → x 0 . Let's now try to understand how this realizes a primitive of the WKB one form at the order of , which is 4P P −5P 2 32P 5/2 . Apparently, it is not integrable at the zero x 0 because it is divergent there. However, we can rewrite it as Note that generically around the zero x 0 , the behavior ofỹ 0 (x) is So the first term in (B.56) is integrable, and we can shrink the contour γ abc (x) to the zero. The second term is badly divergent but a total derivative Therefore, we have found a good primitive of p 1 (x), which is the order part of p(x, ) This finishes the proof of (B.49) at the order of . This way of separating divergent total derivative from the less divergent part in (B.56) is very reminiscent of a more well-known way [75,76] to evaluate the contour integral, which goes as follows. Write P (x) = V (x) − E.
Notice that the first term is integrable at the zero and the second term is a total derivative.
Near a double zero. If the quadratic differential P (x) has a double zero at x 0 , then there are four special WKB lines emanating from x 0 , around which we can access four small solutions, which are denoted as ψ a , ψ b , ψ c , ψ d in the figure 4. We now have linear relations

JHEP01(2022)175
In short, we have a good coordinate in all four sectors. In particular, that means we could determine this way the asymptotic expansion of the cross-ration χ. On the other hand, that asymptotic expansion is already computable from a contour integral of p(x) on a contour wrapping around the double zero while keeping away from it.
We can now relate y to s ab and to the other WKB coordinates nearby: If we expand the left hand side in an asymptotic expansion at large y(x), we can obtain the expansion of e s abc (x) as follows. Asymptotically at large y, where we have parameterized the → 0 asymptotics y(x) = −1/2 (ỹ 0 (x) + 2ỹ 2 (x) + O( 4 )) and a( ) = a (0) + 2 a (2) + O( 4 ). We then obtain the expansion of s abc (x) Let's note a crucial point of (B.69). While (ψ c , ψ a ) cannot be computed by a naive WKB contour integral away from the zeroes, everything else in (B.69) can be in principle evaluated: (A 1 , A 2 ) is normalized to be −i; (A 2 , A 0 ) will be computed in (B.104); (ψ b , ψ c ) is controlled by a WKB contour integral. Therefore (B.69) provides an interesting prediction of (ψ c , ψ a ). So the relation above should really be written as Note that A 1 (y(x)) and A 0 (y(x)) have to be proportional to y (x)ψ b (x) and y (x)ψ a (x), respectively. Therefore to figure out (ψ c , ψ a ), which is x independent, we just need to figure out the constants of proportionality. We can do so by comparing the expansions and read out the x independent terms.

B.3 Recipe for evaluating Wronskians
3. If two solutions ψ n and ψ m are not connected via special WKB lines to the same zero, we can use Plücker relation to reduce to the previous case.
4. What remains is to figure out the Wronskians between local solutions (A k;a , A k,b ). This will be done via perturbation theory in section B.5.

B.4 Numerical implementation
Here in this section, we explain how we implement the numerics. In particular, given an ODE we would like to find the corresponding small solutions and evaluate the Wronskians between them. Let's first consider the case where P (x) is a polynomial of degree n and t(x) = 0. In this case the ODE is regular everywhere on the complex plane with n + 2 asymptotic direction towards the irregular singularity at infinity. Then the task would be to find the decaying solutions along each asymptotic direction. However, initial value problems are more natural in numerics, where one usually specifies the initial condition (the value of ψ and ∂ x ψ at a chosen initial point) and numerically integrate outward along a certain direction. An obvious way to proceed is the so-called shooting method, which reduces the boundary value problems to initial value problems and one adjusts the initial condition until the desired decaying asymptotics is reached. However, it turns out that, in practice, large x asymptotics is very sensitive to the initial condition at small x thus it is very hard to reach a decent accuracy.
Instead, we employ the inward integration approach where the boundary condition at a certain large value of x is provided by the chosen WKB solutions. Intuitively this works better for us because the unwanted dominant solution is suppressed by the inward integration. Indeed, as shown in figure 6, numerical error is greatly reduced this way. We will see more examples of the numerical calculation below.
It's not hard to imagine that dealing with small is challenging for numerics since it exponentially suppresses the solution. This can be easily resolved by a rescaling of the coordinate. For example, under a change of coordinate y = −2/5 x Therefore, the result only depends on the combination a − 2 5 . In the numerics we will study the latter and vary a. Wronskians and cross ratios will be invariant under the rescaling JHEP01(2022)175 of the coordinate. If one wants to study the wavefunctions, we can easily restore the dependence by going back to the original coordinate.
There are other numerical implementation methods available. 12 See, e.g. [77] for a recent study.

B.5 Solutions near zeros and their Wronskians
In general, given a local form around a zero of generic integer order k where the stress tensor is regularỹ Let's attempt to solve the ODE perturbatively in γ. At the leading order, the Schrödinger operator is just ∂ 2 y − y k . A set of nice solutions has been given in [19], which we now review. We choose the solution that decays along the positive real axis which takes the form with large y asymptotics We can produce more solutions by a rotation The inclusion of the factor e − πi k+2 a is such that asymptotically along the ray of e − 2πi Setting k = 1, we are reduced to the Airy functions To obtain higher order corrections, we parameterize the solution as where for now we suppressed the subscript. The differential equation ∂ 2 y A(y) = T local k (y)A(y) is expanded order by order in γ as We fix the normalization of the solutions by matching with the WKB asymptotics, which is uniquely defined as 1 ±2∂ y S e ∓S (B.89) where S = 2 k+2 y k+2 2 + . . . is a (fractional) power series of y, with no constant term. When the zero of even order, there is also log y, which we choose to be the principal branch. For example, We can now give the prescription for determining the solutions order by order in γ. The small A k;a (y), which decays along the ray of exp(− 2πi k+2 a), is given by k;a (y) + γA (1) k;a (y) + γ 2 A k;a (y) + . . . (B.91)

JHEP01(2022)175
A k;a can be obtained recursively in the following way. We already defined the leading order A (0) k;a above. Obviously WKB solutions (B.89) expand in power of γ starting from γ 2 , so we need to choose A (1) k;a = 0. And for each A (n) k;a with n ≥ 2, we need to solve an inhomogeneous ODE with two integration constants to fix, one of which is fixed by requiring the solution to decay along the chosen direction and the other one is fixed to match the normalization of the WKB solution (B.89). In practice, as seen in the examples below, this is achieved by choosing the lower limit of the integration to be at infinity.
Given two small solutions A k;a = A k;a +. . . and In particular, we can verify that (A k;a , A k;a+1 ) = (A k;a+1 ) = −i with no higher order corrections. More interesting ones are the Wronskians between non-consecutive small solutions. We will give concrete examples below for the double zero and cubic zero, since it is trivial for simple zero.
If, on the other hand, we are interested in the local coordinate system around a zero that is a singularity of trivial monodromy, we would havẽ k+2 y k−2 + · · · + a j 2j+4 k+2 y j + · · · + a 0 4 k+2 + l(l + 1) where t(x) = a(x) 2 + ∂ x a(x), and a(x) has residue −l at x 0 . The perturbative solutions can be found in a similar procedure as above except that we need to start with different solutions at the leading order. Since the Schrödinger operator at the leading order is ∂ 2 − y k − l(l+1) y 2 , the solutions that agree with the asymptotics (B.80) and (B.83) are given by k,l;0 (e 2πi k+2 a y) (B.96) whose Wronskians are given by .

(B.97)
One sanity check is to look at i(A k;a+k+2 ). One would like this to be zero since there is a unique decaying solution along a certain ray and therefore two must be proportional to each other. This is indeed mostly true since 2l + 1 ∈ Z. It fails when 2l + 1 is an integer multiple of k + 2, where we have This is another manifestation of the requirement that 2l ≤ k.

B.5.1 Example: double zero
Consider the ODE where a = a (0) + 2 a (2) + 4 a (4) + . . . . To find perturbative solution, we rewrite using γ = 1/2 andỹ = γy, and we have On the other hand, the ODE 13 (B.99) can be solved exactly by using Their Wronskians can be evaluated easily. For example, i(Ã n ,Ã n+1 ) = 1, The cross ratio which exactly coincides with the contour integral exp p(x, )dx around this double zero of the WKB momentum which has a pole

B.5.2 Example: cubic zero
In the local coordinate system around a cubic zero, T (ỹ) =ỹ a = a (0) + 2 a (2) + 4 a (4) + . . . (B.108) 13 The reason we solve (B.99) instead of the one in y coordinate is because we have 2 in the former, so that we can apply the trick → e −iπn to find other solutions. Due to the same reason, we don't shift → e −iπn in the Jacobian part −1/4 if one wants to go back to y coordinate.

JHEP01(2022)175
or equivalently the stress tensor is y 3 +bγ 3 y +aγ 2 with γ = 2/5 . The leading order solutions are defined in (B.79) and (B.81). Since here in this section we only flesh out details for three solutions, for convenience we write φ 0 = A (0) 3;−1 . Let's now solve the ODE perturbatively using the prescription described in B.5. Let's denote the small solution along the ray e −i2πn/5 as ψ(x) − 1 2 P (x). We don't show the ∆ error for the WKB asymptotic solution since the error is too big. The legend of coloring is shared in both diagrams.

JHEP01(2022)175
We will choose the constant A = 0 in this article to normalize the WKB solutions. With this normalization, the large x asymptotics will take the form of only power of x without any constant 14 x ∞ x 2 − 2dx ∼ We can also find this result via going to the local coordinate system. Essentially one needs to figure out the constant of proportionality in (∂ x y(x, )) −1/2 A(y(x)) ∝ ψ(x). To this end, we look at the large x asymptotics of both sides. By definition, the large x asymptotics of ψ(x) is given by 14 There are some other natural choices as well. For example, we chose A = − log √ 2e in [19] such that the regularized integral  There is one simple zero, one double zero and five asymptotic directions. The Stokes diagram is shown in figure 7.
In the local coordinate system around the double zero, therefore we have to choose a (0) = 7i 64 . And since typically y 0 (x) ∼ αx + . . . , y 2 (x) will be nonzero at the zero x = 0.

(B.132)
But this is not problematic because the nontrivial factor actually cancels in the Wronskian. Therefore, they give the same i(ψ −1 , ψ 1 ) = 1 + . . . , as expected. This Wronskian actually has nontrivial corrections, which we will present below. First, of course, we can integrate the WKB momentum along the generic WKB line One can obtain the same result in local coordinate system around the simple zero and the double zero.

B.7.2 Matching around the zero
When t(x) = 0, after shifting the coordinate, the stress tensor looks like In the local coordinate around the zero, y k + a k−2 γ k y k−2 + · · · + a j γ 2+j y j + · · · + a 0 γ 2 (B.146)  Figure 10. In the top figure, g is assumed to be a some order 1 constant, independent of θ. x −∞ (θ) is farther away from − 1 g as θ → ∞. δ is the local variation around x −∞ (θ) that is complex. So it doesn't have to be on the real axis. In the bottom figures, the red line is an example of g eff discussed in this section, namely an example of circular RG flow.

JHEP01(2022)175
Or in terms of g eff (θ), it starts with g = g eff (θ = 0) that has some small imaginary part, and circles around in the complex g eff (θ) plane and goes back to the zero g eff (θ → +∞) → 0 − . Therefore the careful solution we found in (B.150), especially the imaginary part of x −∞ is just to make sure we choose this circular type of RG flow, depicted as red lines in figure 10.

JHEP01(2022)175
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.