A_(2|1) spectral equivalences and nonlocal integrals of motion

The existence of common symmetries underlying both ordinary differential equations (ODEs) and quantum integrable models (QIMs) was first observed for the simple anharmonic oscillator Schrödinger potential in [1], from which functional relations characteristic to QIMs were constructed. Soon after, in [2], the addition of a centrifugal term to the differential problem was shown to correspond to a twist in the boundary conditions of the lattice model. This allowed one to relate a differential eigenvalue problem to more general values of the Fock vacuum parameter in the QIM. The extension to excited states of the Virasoro module was achieved in [3] by introducing a singular part to the potential with poles located away from the origin. The overall symmetry of the differential problem in all cases, although hidden, is encoded by the A₁ = su(2) Lie algebra [4–6]. Besides the twisted XXZ spin chain, the six-vertex lattice model and the massless twisted sine-Gordon integrable field theory can also be shown to be spectrally equivalent to this Schrödinger equation. The crucial ingredients in both frameworks are the Symanzik symmetry of rotations on the complex plane in the case of the ODE and the existence of a Yang–Baxter structure for the QIM. These two elements, although not directly related, lead to equivalent functional relations—the Bethe Ansatz equation (BAE), the T–Q or the quantum Wronskian relations—which allow one to establish a parallel between the transfer matrix T(E) and Baxter's Q(E) operator to spectral determinants D(E) and Wronskians W(E) of solutions to the ODE, e.g. [7]. Consequently, the eigenvalues of Q(E) have zeros at the same points of the spectral parameter E as D(E). Thus, lattice models, interesting from the statistical viewpoint, and integrable quantum field theories (IQFTs), possibly minimal models in conformal field theories (CFTs), can be conveniently mapped onto ODEs defined on the complex plane.

In order to investigate the connection of the integrals of motion in the QIM to the differential problem, a simpler extension consists of utilizing the supersymmetric A_(2|1) = su(2|1) Lie algebra. The advantage of this relies on the existence of a well-understood quantum field theory, namely the AKNS soliton hierarchy [8, 9]. In this situation, we have close links to the integrable perturbed hairpin model and the Perk–Schultz model. Another positive aspect is that the su(2|1) algebra, being one of the simplest supersymmetric extensions of the Lie algebra, provides a good environment to test ideas in supersymmetry and superconformal field theory, which could play an important role in the fundamental descriptions of nature, according to AdS/CFT correspondence [10]. The problem in discussion is not only closely related to physical applications but also has connections to more purely mathematical questions such as number theory. In fact, the eigenvalues of this problem can be combined in spectral Zeta functions, e.g. [11], allowing one to obtain some of the nonlocal charges [12]. With this in mind, the ODEs and integrable models can therefore be used to infer information about the Zeta functions. Another important by-product of the investigation on the aforementioned correspondence came in [13], what is generally considered a reality proof for the spectra of a class of non-Hermitian $\mathcal {PT}$-symmetric Hamiltonians [14, 15].

In what follows, to fix the notation we will present the class of ODEs of interest and review how to construct spectral functional relations from its symmetries. With a few identification, we can show that this expansion corresponds to that associated with the quantum AKNS field theory and that the coefficients in this series are the eigenvalues of its local integrals of motion I_n. The nonlocal integrals of motion $\tilde{H}_n$, usually obtained more inconveniently, are constructed through a novel perturbative analysis of the ODE after a suitable change of variables and the first few are explicitly presented.

Interested in studying an ODE which could encode all the information about the quantum AKNS field theory, we start with the following differential equation:

$$\begin{eqnarray} &&\left( - \frac{{\rm d}^2}{{\rm d}y^2} +4 (P^2 + \bar{P} \kappa {\rm e}^y)+ \kappa ^2 ({\rm e}^{2y} + {\rm e}^{2 \beta ^2 y}) \right) \Phi (y) =0, \end{eqnarray} \tag{ 2.1 }$$

defined originally on the whole real line −∞ < x < ∞ for a restricted interval of the coupling $0 \le \beta ^2 \le \frac{1}{2}$. This is the equation presented in [8] if we identify the parameters there as $p = 2 \imath P , q = 2\bar{P}, n = -2\beta ^2$. For our convenience, we introduce a preliminary change of variables

$$\begin{eqnarray} &&y = \log ((-E)^{-\frac{M+1}{2M}}x^{M+1}), \qquad \Phi (y) = \frac{(M+1)^2}{\sqrt{y}} \psi (y), \end{eqnarray} \tag{ 2.2 }$$

and a consort reparameterization given by

$$\begin{eqnarray} &&P = \frac{2\ell +1}{4(M+1)}, \quad \bar{P} = \frac{\alpha }{4(M+1)}, \quad \beta ^2=\frac{1}{M+1}, \quad \kappa = \frac{(-E)^{\frac{M+1}{2M}}}{M+1}. \end{eqnarray} \tag{ 2.3 }$$

This allows us to re-express the so-called integrable perturbed hairpin model (ODE) (2.1) in a more familiar fashion, namely a generalization of the harmonic oscillator,

$$\begin{eqnarray} &&\left( - \frac{{\rm d}^2}{{\rm d}x^2} + { \frac{\ell (\ell +1)}{x^2} } + {x^{2M}} + {\alpha x^{M-1}} - E \right) \psi (x) =0, \end{eqnarray} \tag{ 2.4 }$$

defined on the half-line only, 0 < x < ∞, within the range M ⩾ 1, for which the harmonic oscillator is a limiting case. The appearance of the α-term was shown in [16] to correspond to breaking the su(2) symmetry into su(2|1). The complete knowledge about the spectrum {E_k} can give information about the ground state of the associated quantum integrable model¹.

The equation above is an extension of the x^2M potential, or of the $\mathcal {PT}$-symmetric non-Hermitian (ıx)^2M problem [14, 15], so that when x → ix, the associated eigenfunctions and eigenvalues are given by $\tilde{\psi }(x) = \psi ({\rm i} x), \lambda = -E$. The investigation of this family of potentials has motivated a number of studies on non-Hermitian quantum mechanics, and for the non-supersymmetric problem, α = 0, the reality proof for its spectra is given in [13], whereas metrics and Hermitian counterparts for the −x⁴-potential are found in [17]. More general non-Hermitian Hamiltonians have been studied outside the scope of the ODE/IM correspondence in e.g. [18–20].

The M = 1 case usually requires special attention as some of the results no longer hold and have to be modified. However, the introduction of an α ≠ 0 term has a unifying character, as the energy can be seen as rescaled by an amount α, in this case. Therefore, the usual readaptation one has to do when evaluating the extended harmonic oscillator can be easily obtained by simply taking

$$\begin{eqnarray} \psi (x; M = 1,\ell ,\alpha =-\epsilon , E= 0) &=& \psi (x; M =1,\ell ,\alpha =0 ,\epsilon ). \end{eqnarray} \tag{ 2.5 }$$

This model (M = 1, α = 0) is sometimes denoted the Hooke atom because it describes a particle in the presence of a centrifugal form alongside a linear Hooke force, and it is a benchmark for being exactly solvable. In the integrable model framework, this corresponds to a coupling of $\beta ^2=\frac{1}{2}$, known as the free fermion point because in this situation the corresponding sine-Gordon field theory can be mapped onto a non-interacting massive Thirring model. The solution of the Schrödinger problem is given in terms of Laguerre polynomials and hypergeometric functions, depending on the boundary conditions and eigenvalues, weighted by a decaying exponential factor. For general M ≠ 1, the quantization problem cannot be solved exactly in a closed form unless E = 0. In this case,

$$\begin{eqnarray} \qquad \nonumber \psi (M,\ell , \alpha ,E = 0,x) &\propto & x^{\frac{1}{2} (1\pm (2 \ell +1))} {\rm e}^{-\frac{x^{M+1}}{M+1}} L_{-\frac{M+1\pm (2 \ell +1)+\alpha }{2 (M+1)}}^{\frac{2 \ell +1}{M+1}}\left(\frac{2 x^{M+1}}{M+1}\right),\\ \qquad \psi (M,\ell , \alpha ,E = 0,x) &\propto & x^{\frac{1}{2} (1\pm (2 \ell +1))} {\rm e}^{-\frac{x^{M+1}}{M+1}} U_{\frac{M+1\pm (2 \ell +1)+\alpha }{2 (M+1)}}^{\frac{M+1\pm (2 \ell +1)}{M+1}}\left(\frac{2 x^{M+1}}{M+1}\right), \end{eqnarray} \tag{ 2.6 }$$

where L and U denote the Laguerre and hypergeometric functions.

It is important to note that equation (2.4) is invariant under discrete rotations by multiples of a certain angle on the complex plane, alongside a change in the constants α and E. The so-called Symanzik transformations for this system read [21]

$$\begin{eqnarray} &&x \; \rightarrow \; \omega ^{-k} x, \quad \alpha \rightarrow (-1)^{k}\alpha , \quad E \; \rightarrow \; \omega ^{2k}E, \quad \omega = \exp \frac{\imath \pi }{M+1} . \end{eqnarray} \tag{ 2.7 }$$

Since symmetries in general help one to solve differential equations, we can use this invariance to generate a new solution, suitably normalized,

$$\begin{eqnarray} &&\psi _k (x; M,\ell ,E,\alpha ) \equiv (-1)^{k+1}\imath \, \omega ^{\frac{k}{2}(1-{\rm e}^{\imath k \pi }\alpha )} \psi (\omega ^{-k}x;M,\ell ,\omega ^{2k}E, {\rm e}^{\imath k \pi }\alpha ), \end{eqnarray} \tag{ 2.8 }$$

which solves the ODE for any $k \in \mathbb{Z}$. It is the unique subdominant solution, specified by the condition that it decays exponentially at infinity, in the region $S_k= \lbrace | \textrm {arg}(x) - \frac{k\pi }{M+1} | < \frac{\pi }{2(M+1)} \rbrace$. Once one knows a solution of the second-order differential equation, the symmetry on the complex plane can be used to generate a new linearly independent solution, and one can use these two to write any solution of the problem. However, allowing (2.4) in the complex plane requires a consistent choice of boundary conditions, as square-integrable wavefunctions can only exist in appropriate regions. The Stokes sectors in the complex plane define wedges where eigenfunctions exponentially decay and blow up in an alternating pattern. The problem we are mainly interested in throughout this paper is the radial problem, defined on the positive real axis and imposed to vanish at the origin and infinity, within the sector S₀.

The kth rotated solution, ψ_k, is subdominant in S_k and dominant in S_{k + 1}. Moreover, {ψ_k, ψ_{k + 1}} forms a basis of solution and W[ψ_{k − 1}, ψ_{k + 1}] will vanish if and only if there exists a solution decaying in both sectors. Thus, the solution ψ₀ = ψ, one may find for the problem, is, for instance, linearly independent of ψ_±1, solutions on the complex plane rotated by an angle ω^±1 with respect to ψ₀. In our problem, defined on the positive real line, with vanishing boundary conditions at the origin and at infinity, ψ₀ is defined on the real axis, whereas ψ_±1 satisfy boundary conditions on certain regions of the complex plane. To be more precise, the appropriate boundary condition at infinity for our radial problem implies that it decays exponentially according to

$$\begin{eqnarray} &&\psi (x; M,\ell ,E,\alpha ) \sim \Psi (x; M,\ell ,E,\alpha ) \equiv x^{-\frac{M+\alpha }{2}}\exp \left( -\frac{x^{M+1}}{M+1} \right). \end{eqnarray} \tag{ 2.9 }$$

Note that when M is not an integer, the potential V(x) is not single valued, so a branch cut, chosen along the negative real axis, is necessary.

Another symmetry of the general ODE (2.4) is that it is invariant under ℓ → −(ℓ + 1). Therefore, there are two linearly independent solutions, $\psi ^\pm (x; M,\ell ,E,\alpha ) \equiv \psi (x; M,\frac{1}{2} \pm \frac{2\ell +1}{2},E,\alpha )$. The parameter ℓ dictates the behavior near the origin and the expressions above indicate that, although (2.6) has a singular behavior near the origin of the type x^−ℓ, it can be easily analytically continued to a regular solution. Thus, by simply using this transformation, one recovers square integrability, an important property when considering the hermiticity of the differential problem. Near the origin, the solutions behave like

$$\begin{eqnarray} &&\psi ^{\pm }(x; M,\ell ,E,\alpha ) \sim \varphi ^{\pm } (x;M,\ell ,E,\alpha ) \equiv x^{\frac{1}{2} \pm \frac{2\ell +1}{2}}, \quad \textrm {as} \quad x \rightarrow 0, \end{eqnarray} \tag{ 2.10 }$$

and one can use them as a basis of expansion for a solution in this region. Because the problem is not solved exactly in general for M ≠ 1, one must instead resort to constructing solutions near the origin by the Frobenius method and using the WKB method to find asymptotic solutions which can furnish us with non-perturbative information about its analytical structure.

The basis of solutions φ^± around the origin (2.10) forms a linearly independent set satisfying the following Wronskian identities:

$$\begin{eqnarray} &&\fl W_{[-,+]} = W[\varphi ^-(x; M,\ell ,E,\alpha ),\varphi ^+(x; M,\ell ,E,\alpha )] = 2\ell +1, \end{eqnarray} \tag{ 2.11 }$$

$$\begin{eqnarray} &&\fl (2\ell +1)\psi (x; M,\ell ,E,\alpha ) = W_{[0,+]}(M,\ell ,E,\alpha ) {\varphi ^- (x; M,\ell ,E,\alpha )} \nonumber\\ &&+ W_{[0,-]}(M,\ell ,E,\alpha ) {\varphi ^+ (x; M,\ell ,E,\alpha )}, \end{eqnarray} \tag{ 2.12 }$$

with x-independent coefficients given by the Wronskians W_{[0, ±]} = W[ψ, φ^±] of a solution ψ₀(x; M, ℓ, E, α) with this pair, so that one can use them to expand a general solution in this region. The exact relation between them is given in terms of Wronskians. Considering that W_{[0, 1]} = W_{[ − 1, 0]} can be imposed to be 1 with a suitable normalization condition for ψ_k. In order to remove completely the x-dependence of this relation, one takes the Wronskian of this equation with the pair φ^±, giving

$$\begin{eqnarray} &&\fl W_{[-1,1]}(M,\ell ,E,\alpha ) W_{[0,\pm ]}(M,\ell ,E,\alpha ) = \omega ^{-\frac{\alpha \pm (2\ell +1)}{2}}W_{[0,\pm ]} (M, \ell , \omega ^{-2}E,-\alpha ) \nonumber\\ &&+\, \omega ^{+\frac{\alpha \pm (2\ell +1)}{2}}W_{[0,\pm ]} (M, \ell , \omega ^{+2}E,-\alpha ). \end{eqnarray} \tag{ 2.13 }$$

One of the quantities on the left-hand side involves the asymptotic behavior in the sectors S_±1 of the complex plane and is therefore associated with the lateral problem and the $\mathcal {PT}$-symmetric differential equation; namely W_{[ − 1, 1]} vanishes if and only if, for a certain value E_n, the wavefunctions are such that ψ_±1 vanish asymptotically. Thus, this Wronskian encloses information regarding the quantization condition of this problem. On the other hand, the remaining object W_{[0, ±]} depends on the asymptotic behavior on the real line as well as near the origin, being associated with the radial problem we are interested in.

The objects

$$\begin{eqnarray} &&C(E,\alpha ) = W_{[-1,1]}(M,\ell ,E,\alpha ),\quad D^{(\mp )}(E,\alpha ) = W_{[0,\pm ]}(M,\ell ,E,\alpha ) \end{eqnarray} \tag{ 2.14 }$$

are entire functions and we can use the Weierstrass factorization theorem to express them as an infinite product. A zero of D(E, α) indicates that the eigenproblem has a solution for which the eigenfunction corresponding to the eigenvalue E satisfies the correct boundary conditions. Therefore, we can conclude that this coincides with the spectral determinant, up to a normalizing factor, vanishing at the possible eigenvalues of the differential equation. What is more, the definitions above furnish us with the analogue of Baxter T–Q relations (also in [12] and reference [1] therein)

$$\begin{eqnarray} &&\fl C(E,\alpha ) D^{(\pm )}(E,\alpha ) = \omega ^{+\frac{2\ell +1+\alpha }{2}} D^{(\pm )}(\omega ^{+2}E, -\alpha ) + \omega ^{-\frac{2\ell +1+\alpha }{2}} D^{(\pm )}(\omega ^{-2}E, -\alpha ), \end{eqnarray} \tag{ 2.15 }$$

whose zeros $E=E_j^{(\alpha )}$ are characterized by

$$\begin{eqnarray} &&\frac{D^{(\pm )}\big(\omega ^{-2}E_j^{(\alpha )}, -\alpha \big)}{ D^{(\pm )}\big(\omega ^{+2}E_j^{(\alpha )}, -\alpha \big)} = - \omega ^{(2\ell +1+\alpha )}. \end{eqnarray} \tag{ 2.16 }$$

Part of the solutions for this equation correspond to C(E, α) and the remaining ones are associated with zeros of D(E, α). Because the problem defined on the real line is Hermitian, the zeros of D(E, α) should be real. Therefore, any complex eigenvalue should be a zero of C(E, α), but there is no guarantee that the zeros of C(E, α) cannot be real.

In order to use the Weierstrass factorization, we must take into account that a simple WKB approximation tells us that the zeros accumulate toward infinity along the positive real axis according to $E_n \sim n^{\frac{2M}{M+1}}$, so if M > 1 we can establish a convergent product D^{( − )} ≡ D in the form $D(E,\alpha )= D(0,\alpha )\prod _{k=0}^\infty (1-\frac{E}{E_k^{(\alpha )}})$, which leads immediately to an expression of the form of the BAEs,

$$\begin{equation} \prod _{k=0}^\infty \left( \frac{E_k^{({\alpha })} - \omega ^{-2}E_j^{({-\alpha })}}{E_k^{({\alpha })} - \omega ^{2}E_j^{({-\alpha })}} \right) = - \omega ^{(2\ell +1 + {\alpha })}. \end{equation} \tag{ 2.17 }$$

In this case, they correspond to two sets of coupled equations, both for the values of +α and −α. The problems defined for ±α have different properties, reflecting on the eigenvalues, or Bethe roots. In the differential equation framework, for one case we have a convex potential, whereas for the other a double-well potential [22]. For M = 1, we cannot guarantee the convergence of the infinite product above and changes are expected to be necessary. Relations (2.5), though, make it simple to obtain the correct equations for the T–Q and quantum Wronskian relations. One should note that, rather than the Q(E, α)-operator itself, the spectral determinant D(E, α) should be identified with the quantity A(E, α), defined by

$$\begin{equation} Q ( \omega ^{\pm 2} E,-\alpha )= \omega ^{\pm \frac{2\ell +1+ \alpha }{2}} A (\omega ^{\pm 2}E,-\alpha ). \end{equation} \tag{ 2.18 }$$

To obtain the correct normalization factor to the Hadamard factorization, we observe that the solution of the differential equation at vanishing value of energy can be found exactly in (2.6). By normalizing (2.9) and (2.10) appropriately, the expansion of the spectral determinant above has the expected leading behavior

$$\begin{eqnarray} &&\log D(E,\alpha ) \sim \frac{\pi }{\cos \left(\frac{\pi }{2 M}\right)} \left(-\frac{E }{\nu }\right)^{\frac{M+1}{2 M}}, \quad \nu \equiv 2\sqrt{\pi } \frac{\Gamma \left(\frac{3M+1}{2M}\right)}{\Gamma \left(\frac{2M+1}{2M}\right)} ^\frac{2M}{M+1}, \end{eqnarray} \tag{ 2.19 }$$

well known in the literature. An important result of the spectral determinant expansion in the spectral parameter above is the appearance of quantities which correspond to the eigenvalues of the commuting set of local conserved operators in the quantum AKNS model. The reformulation here presented is well known now and allows one to transform our original differential problems into a set of coupled nonlinear equations. In the following section, we will present the exact identification of such coefficients with the integrals of motion corresponding to the quantum integrable model and giving explicit expressions for them in terms of integral equations.

It is an important result that one can transform our initial eigenvalue problem into algebraic equations, but the simplification is only relative as the number of equations to solve becomes infinite. In any case, this alternative method which is based not on the differential equation directly but rather on the functional relations one can construct from it allows one to extend the toolbox to solve the original eigenvalue problem as well as gives access to important information in an exact analytic form. In the thermodynamic limit, this approach in essence transforms the differential problem into nonlinear integral equations, usually referred to as NLIEs. This has a central role in our strategy to compute nonlocal conserved quantities for the QIM from an analysis of the ODE.

The starting point is to consider the C–D relation (2.15), which is equivalent to the T–Q relation,

$$\begin{eqnarray} T(E ,\alpha )Q(E ,\alpha ) &=& Q (\omega ^{2} E,-\alpha )+ Q (\omega ^{-2} E ,-\alpha ), \end{eqnarray} \tag{ 3.20 }$$

if we identify the spectral determinant D(E, α) with the quantity A(E, α) in (2.18). One then defines a new object,

$$\begin{eqnarray} &&a(E,\alpha ) = \omega ^{2\ell +1+\alpha }\frac{D(\omega ^{+2} E(\alpha ),-\alpha )}{D(\omega ^{-2}E(\alpha ),-\alpha )}, \end{eqnarray} \tag{ 3.21 }$$

so that for the zeros E_k(α) of A(E(α), −α), the associated Bethe equation may be recast in the form of a(E_k, α) + 1 = 0. This is actually the same system of two coupled equations: one for +α and one for −α. We make the assumption that T(E, α) and A(E, α), or equivalently C(E, α) and D(E, α), respectively, are entire functions so that we can factorize them in terms of their zeros in a convergent product through a Weierstrass decomposition.

One can write the quantization condition as

$$\begin{eqnarray} &&\fl \log a(E,\alpha ) = \log \left( \omega ^{2\ell +1+\alpha } \right) + \sum _{k=0}^\infty \log \left( \frac{1- \omega ^2 \frac{E(\alpha )}{E_k(-\alpha )}}{1- \omega ^{-2} \frac{E(\alpha )}{E_k(-\alpha )}} \right) =2 {\rm i} \pi \left(m_k+\frac{1}{2}\right), \end{eqnarray} \tag{ 3.22 }$$

so that for each eigenvalue E_k, there are infinite possibilities for the integer parameter m_k. When we take the logarithm of the expression above, we introduce an arbitrariness in the choice of the branch, or the Riemann sheet, denoted by an integer m_k, in such a way that an eigenvalue E_k is associated with one of the integers m_k. A certain eigenstate of the Q operator is characterized by a unique set of integers {m_k}. However, not every set of integers will correspond to an eigenstate of Q.

We can replace the infinite sum over zeros by a contour integral. The total contour can be split into two parts: one part containing the positive real line and one contour for the remaining complex plane. Using the reparameterization $E ={\rm e}^{\frac{2M}{M+1}\theta }$, the logarithm expansion of the quantity a(θ, α) can be expanded as an integral which can be written as a sum of infinitely many residues plus a diverging contribution at the large circle [23, 24], after the inclusion of a term responsible for the leading behavior following (2.19):

$$\begin{eqnarray} &&\fl \log a(\theta , \alpha )= \log (\omega ^{2\ell +1+ \alpha }) -\frac{2 \pi \imath }{\nu ^{\frac{M+1}{2 M}} }{\rm e}^{\theta } +\int _{c_1} K_+(\theta -w) \log (a(w, \alpha )+1) \, {\rm d}w \nonumber\\ &&-\int _{c_2} K_+(\theta -w) \log \left(\frac{1}{a(w, \alpha )}+1\right) \, {\rm d}w \nonumber\\ &&+\int _{c_1} K_-(\theta -w) \log (a(w,-\alpha )+1) \, {\rm d}w \nonumber\\ &&-\int _{c_2} K_-(\theta -w) \log \left(\frac{1}{a(w,-\alpha )}+1\right) \, {\rm d}w , \end{eqnarray} \tag{ 3.23 }$$

where the contours c₁ and c₂ are placed just above and below the real axis and the exact forms of the kernels K_±(θ) were given in an integral form in [25, 26]. These integrals can be written as infinite sums, so the kernel may be expressed as

$$\begin{eqnarray} &&K_{\pm } (\theta ) = \case{1}{2} (\varphi (\theta )\pm \chi (\theta )), \end{eqnarray} \tag{ 3.24 }$$

with

$$\begin{eqnarray} \varphi (\theta ) &=& \sum _{n=1}^{\infty } \frac{M}{\pi }{\rm e}^{-2 \theta M n} \tan (\pi M n)- \sum _{n=1}^{\infty } \frac{1}{\pi }{\rm e}^{\theta (1-2 n)} \cot \left(\frac{\pi (2 n-1)}{2 M}\right), \end{eqnarray} \tag{ 3.25 }$$

$$\begin{eqnarray} \chi (\theta ) &=& \sum _{n=1}^{\infty } \frac{M}{2 \pi } {\rm e}^{\theta M (1-2 n)} \cot \left(\frac{\pi }{2} M (2 n-1)\right)- \sum _{n=1}^{\infty } \frac{1}{2 \pi } {\rm e}^{-2 \theta n} \tan \left(\frac{\pi n}{M}\right). \end{eqnarray} \tag{ 3.26 }$$

Therefore, the sums over the residues only give an asymptotic series expansion and from them we can construct the conserved charges. The NLIE for the KdV su(2) problem, with vanishing α, involves a function which coincides with the soliton–soliton scattering amplitude for the sine-Gordon model. The same NLIE was earlier derived for the BAE associated with the XXZ model and by taking the appropriate continuous limit leads to the sine-Gordon field theory. Thus, the integrals of motion of the aforementioned models are related to each other. In this case, the scattering matrices are the ones that should appear in the quantum AKNS field theory and the conserved charges are equivalent in both models.

After some identifications, we can relate the integrals of motion to the integrals appearing in the expansion. In the case of the nonlocal ones, we have

$$\begin{eqnarray} &&\fl \tilde{I}_{2 n} = -\frac{1}{4 \pi }\left(\frac{M+1}{2^{M}}\right)^{2 n}\frac{M}{\cos (\pi M n)} \left( \int _{c_1} {\rm e}^{2 n M {\rm i} w} \log \left({a(\alpha ,w)}+1\right) \, {\rm d}w \right.\nonumber\\ &&\left. - \int _{c_2} {\rm e}^{2 n M {\rm i} w} \log \left(\frac{1}{a(\alpha ,w)}+1\right) \, {\rm d}w +\int _{c_1} {\rm e}^{2 n M {\rm i} w} \log \left({a(-\alpha ,w)}+1\right) \, {\rm d}w \right.\nonumber\\ &&\left. - \int _{c_2} {\rm e}^{2 n M {\rm i} w} \log \left(\frac{1}{a(-\alpha ,w)}+1\right) \, {\rm d}w \right), \end{eqnarray} \tag{ 3.27 }$$

$$\begin{eqnarray} &&\fl \tilde{I}_{2 n-1}=-\frac{1}{4 \pi }\left(\frac{M+1}{2^{M} }\right)^{2 n-1}M \cos \left(\frac{\pi M (2n-1)}{2}\right) \left( \int _{c_1} {\rm e}^{(2 n-1) M {\rm i} w} \log ({a(\alpha ,w)}+1) \, {\rm d}w \right.\nonumber\\ &&\left. - \int _{c_2} {\rm e}^{(2 n-1) M {\rm i} w} \log \left(\frac{1}{a(\alpha ,w)}+1\right) \, {\rm d}w \right.\nonumber\\ &&\left. -\int _{c_1} {\rm e}^{(2 n-1) M {\rm i} w} \log \left({a(-\alpha ,w)}+1\right) \, {\rm d}w \right.\nonumber\\ &&\left. + \int _{c_2} {\rm e}^{(2 n-1) M {\rm i} w} \log \left(\frac{1}{a(-\alpha ,w)}+1\right) \, {\rm d}w \right) \end{eqnarray} \tag{ 3.28 }$$

and similar expressions for the local ones [27, 28]. We note that in order to have a consistent notation with respect to the literature, it is convenient to make the following identifications: $\tilde{H}_{n}=\tilde{I}_{2 n}$ and $\tilde{H}_{\frac{2 n-1}{2}}=\tilde{I}_{2 n-1}$ for when α ≠ 0. For vanishing α, the contributions of the integrals coincide, so the last two terms in the series above cease to exist and one recovers the results of [4–6, 27]. A WKB analysis of (2.4) leads to an expansion of the spectral determinant in terms of local conserved charges which can be found in e.g. [3, 27] when α = 0. The local charges for the supersymmetric case are investigated in detail in [28] but also have both odd and even contributions, with the latter vanishing with α, as expected.

The expressions above though provide a recursive scheme which is efficient to make numerical comparisons, they are also an important tool from a practical viewpoint, due to the difficulty of computing the objects I_n and—especially—$\tilde{H}_n$, the local and nonlocal conserved charges of the model. Having in mind the equivalence between Baxter's Q(E) operator and the spectral determinant D(E), we will show in the following section how the nonlocal integrals of motion can be determined in a reasonable way.

We have seen that the supersymmetric model studied here has an important feature, namely it interpolates the properties of the quantum KdV field theory or, equivalently, that of the sine-Gordon model and the XXZ spin lattice, which are all well understood. Besides being characterized by local integrals of motion of the kind previously presented for the AKNS analogue, there is another set of important commuting operators, the denoted nonlocal integrals of motion. Originally defined to be the coefficients in a Taylor expansion of the transfer matrix operator, e.g. [5, 29], they can be expressed as

$$\begin{eqnarray} &&T(\lambda )=2 \cos (2 \pi p)+\sum _{n=1}^{\infty } G_n \lambda ^{2 n}, \end{eqnarray} \tag{ 4.29 }$$

with λ being the spectral parameter related to E according to E = λ²ρ, where $\rho =2^{\frac{2 M}{M+1}} (M+1)^{\frac{2 M}{M+1}} \Gamma (\frac{M}{M+1})^2$.

An alternative way of expressing the nonlocal conserved charges, in the case of an underlying su(2) symmetry, consists of defining them as H_n in terms of the coefficients of an expansion of the Q-operator in powers of the variable $\lambda \frac{\Gamma (1-\beta ^2)}{\beta ^2}$. The dual nonlocal integrals of motion, specified by $\tilde{H}_n(\beta ^2)=H_n(\frac{1}{\beta ^2})$, give rise to an expansion of the operator A(E), which corresponds to the spectral determinant D(E) of the associated Schrödinger problem. This means that we can write down an asymptotic expansion, possibly with a limited radius of convergence,

$$\begin{eqnarray} &&\log D(E) \sim \sum _{n=1}^{\infty } \left( \frac{2^M}{M+1} \right)^{2n} \tilde{H}_n (-E)^{-n (M+1)}, \end{eqnarray} \tag{ 4.30 }$$

for large values of the spectral parameter E.

However, we have seen that the addition of the α-term brought another set of local integrals of motion, characterized by even indices, to the pre-existing ones, specified by odd indices only. Therefore, we would expect the same to happen with the nonlocal charges. This implies having an expansion

$$\begin{eqnarray} &&T(\lambda ,\alpha )=T(0,\alpha ) + \sum _{n=1}^{\infty } G_\frac{n}{2} \lambda ^{n}, \end{eqnarray} \tag{ 4.31 }$$

which reduces to (4.29) when all even powers of λ vanish with α. So the first nonlocal conserved charge $G_\frac{1}{2}$ appears in the power of λ, rather than G₁ at order λ², and has no equivalence when α = 0, unlike G₁. Moreover, the dual nonlocal charges also appear as twice as many,

$$\begin{eqnarray} &&\log D(E,\alpha ) \sim \sum _{n=1}^{\infty } \left( \frac{2^M}{M+1} \right)^{n} \tilde{H}_{\frac{n}{2}} (-E)^{-n \frac{(M+1)}{2}}, \end{eqnarray} \tag{ 4.32 }$$

and the first integral of motion in the expansion would have no counterpart in the su(2) scenario and $\tilde{H}_1$ would be the first nonvanishing contribution, since no $\tilde{H}_{\frac{1}{2}}$ would be present in that case.

The nonlocal integral equations allow us to compute the charges numerically but analytic expressions are difficult to obtain, particularly in the case of the nonlocal integrals of motion. In [4–6], the G_n are formally computed in generality for α = 0 and exactly evaluated for first two contributions in the expansion. It would be in principle advantageous if one could determine these important objects by investigating the associated ODE with nonvanishing α by the existing spectral correspondence. In fact, here we show that these properties of the operators A(E, α) and Q(E, α) can be obtained in a systematic fashion by investigating the same original differential equation (2.4). Although the construction of the local integrals of motion can be done directly from the latter, in the case of nonlocal charges some reparameterizations prove to be very convenient. For that, we introduce the new variables,

$$\begin{equation} s= \log \big( \pm \varepsilon ^{-\frac{M}{M+1}} x \big), \quad \Sigma (s)= \frac{1}{\sqrt{x}}\psi (x), \quad \mathcal {F} = \varepsilon ^M = (\pm E)^{-\frac{M+1}{2}}. \end{equation} \tag{ 4.33 }$$

This transforms the problem into an equivalent one, but free from singularities at the origin again,

$$\begin{eqnarray} &&\fl \left[ -\frac{{\rm d}^2}{{\rm d}s^2}+ \left( \frac{2\ell +1}{2}\right)^2 \mp \,{\rm e}^{2s}+ (\pm )^{M+1} \alpha {\rm e}^{(M+1)s}\mathcal {F} + {\rm e}^{2(M+1)s}\mathcal {F}^2 \right]\Sigma (s)=0, \end{eqnarray} \tag{ 4.34 }$$

which corresponds to a potential expanded in terms of powers of $\mathcal {F}$. This can be associated with the AKNS problem by taking $M=\frac{1-\beta ^2}{\beta ^2},\; \ell = 2\frac{P}{\beta ^2} -\frac{1}{2},\;\ \alpha = 4\frac{\bar{P}}{\beta ^2}$.

For simplicity, we can work with the choice $\varepsilon = E^{-\frac{M+1}{2M}}$. Instead, one could also choose $\varepsilon = \left(-E\right)^{-\frac{M+1}{2M}}$, which is equivalent to replacing s by $s + \imath \frac{\pi }{2}$ in the former case. The two possibilities for the choice of the sign are closely related to the different boundary conditions associated with the radial and lateral $\mathcal {PT}$-symmetric problem. The series structure on the parameter ε tells us that for small values of $\mathcal {F}$, i.e. large values of E, this equation can be solved perturbatively, order by order. The unperturbed potential with $\mathcal {F}=0$, which also corresponds to s → −∞, has two linearly independent solutions given by Bessel functions:

$$\begin{eqnarray} &&\Sigma _\pm (s)= \mathcal {N}_\pm \Gamma \left(1 \pm \frac{2\ell +1}{2}\right) J_{\pm \frac{2\ell +1}{2}}\left({\rm e}^{s} \right), \end{eqnarray} \tag{ 4.35 }$$

satisfying required boundary conditions of vanishing at s = ±∞, respectively, as well as the familiar relations, with $\mathcal {N}_\pm =\pi ^{\pm \frac{1}{2}} \Gamma (\frac{2 \ell +1}{2}+\frac{1\pm 1}{2})^{\mp 1}$,

$$\begin{eqnarray} W[\Sigma _-(s),\Sigma _+(s)] &=& 2,\nonumber\\ \qquad\qquad\Sigma (s) &=& \case{1}{2}D^+(E,\alpha ) \Sigma _+(s) + \case{1}{2}D^-(E,\alpha ) \Sigma _-(s). \end{eqnarray} \tag{ 4.36 }$$

In fact, the solutions Σ_±(s) have the first-order contribution of the type φ^±(x) in (2.10), after the appropriate change of variables.

As s → −∞, we have that Σ₊(s) → 0, so we can establish the proportionality of the spectral determinant D(E, α) = D^{( − )}(E, α) = W_{[0, +]}(E, α) with the wavefunction in this limit. However, expression (4.36) is generally true in this regime and does not rely on any interpretation of the D^±(E, α) as spectral determinants which vanish for appropriate boundary conditions. In fact, in the present context the transformed differential equation does not exactly constitute an eigenvalue problem, albeit closely related to one. Therefore, these Wronskians are simply the projection of a solution in the aforementioned basis and the problem of determining them relies on a solution of the differential equation.

Because the latter can be seen as a perturbed equation in powers of the parameter $\mathcal {F}$, we can assume a general solution of the form

$$\begin{eqnarray} \Sigma (s) &=& \Sigma ^{(+)}(s)+\Sigma ^{(-)}(s), \end{eqnarray} \tag{ 4.37 }$$

$$\begin{eqnarray} \Sigma ^{(\pm )}(s) &=& \Delta _0^{(\pm )}(s)+\sum _{j=1}^{\infty } \mathcal {F}^j \Delta _j^{(\pm )}(s), \end{eqnarray} \tag{ 4.38 }$$

where in the lowest order $\Delta _0^{(\pm )}(s) = \Sigma _{\pm }(s)$ is the solution of the unperturbed equation ($\mathcal {F}=0$) satisfying the necessary boundary conditions. Besides, we have added all the remaining perturbative contributions, so that the result is expressible as a series in integer powers of $\mathcal {F}$. Restricting for now our analysis to the ground state of the Q-operator for the perturbed Hamiltonian (4.34), we can expand the wavefunction in powers of $\mathcal {F}$ and solve the associated Schrödinger equation perturbatively at each order. First,

$$\begin{eqnarray} &&\fl \Delta _1^{(\pm )}(s)= \mp \alpha \frac{ \Gamma \left(1\pm \frac{2 \ell +1}{2}\right)^2}{2 \ell +1} \Sigma _\mp (s) \int _{s_{11}}^s {\rm e}^{(M+1) y} J_{\pm \frac{2 \ell +1}{2}}\left({\rm e}^y\right){}^2 \, {\rm d}y + \nonumber\\ &&\hphantom{\fl \Delta _1^{(\pm )}(s)}\quad\pm \frac{\alpha }{2} \frac{\pi }{ \cos (\pi \ell )} \Sigma _\pm (s) \int _{s_{12}}^s {\rm e}^{(M+1) y} J_{-\frac{2 \ell +1}{2}}\left({\rm e}^y\right) J_{\frac{2 \ell +1}{2}}\left({\rm e}^y\right) \, {\rm d}y; \end{eqnarray} \tag{ 4.39 }$$

then, by solving the resulting equation as the coefficient of $\mathcal {F}^2$,

$$\begin{eqnarray} &&\fl \Delta _2^{(\pm )}(s)= \mp \frac{\Gamma \left(1\pm \frac{2 \ell +1}{2} \right)^2 }{2 \ell +1} \Sigma _\mp (s) \int _{s_{21}}^s {\rm e}^{2 (M+1) y} J_{\pm \frac{2 \ell +1}{2}}\left({\rm e}^y\right){}^2 \, {\rm d}y+ \nonumber\\ &&\hphantom{\fl \Delta _2^{(\pm )}(s)=}- \frac{\alpha }{2} \Gamma \left(\pm \frac{2 \ell +1}{2}\right) \Sigma _\mp (s) \int _{s_{21}}^s {\rm e}^{(M+1) y} \Delta _1(y) J_{\pm \frac{2 \ell +1}{2}}\left({\rm e}^y\right) \, {\rm d}y + \nonumber \\ &&\hphantom{\fl \Delta _2^{(\pm )}(s)=}\pm \frac{1}{2} \frac{\pi }{\cos (\pi \ell )} \Sigma _\pm (s) \int _{s_{22}}^s {\rm e}^{2 (M+1) y} J_{-\frac{2 \ell +1}{2}}\left({\rm e}^y\right) J_{\frac{2 \ell +1}{2}}\left({\rm e}^y\right) \, {\rm d}y + \nonumber\\ &&\hphantom{\fl \Delta _2^{(\pm )}(s)=} - \frac{\alpha }{2} \Gamma \left(\mp \frac{2 \ell +1}{2}\right) \Sigma _\pm (s) \int _{s_{22}}^s {\rm e}^{(M+1) y} \Delta _1(y) J_{\mp \frac{2 \ell +1}{2}}\left({\rm e}^y\right) \, {\rm d}y \end{eqnarray} \tag{ 4.40 }$$

and so on, successively, for higher orders. We see that Δ₂(s) depends explicitly on Δ₁(s) and the same pattern occurs at all orders: to determine the solution at each order, it is necessary to know the solutions at previous ones. In fact, the solution for a generic order can be found in a closed formula described by the expression

$$\begin{eqnarray*} &&\fl \nonumber \Delta _{k}^{(\pm )}(s) = - \frac{\Gamma \left(\frac{2 \ell +1}{2} \right)}{2} \Sigma _-(s) \int _{s_{k1}}^s {\rm e}^{(M+1) y} J_{\frac{2 \ell +1}{2}}\left({\rm e}^y\right) \left( {\rm e}^{(M+1) y} \Delta _{k-2}^{(\pm )}(y)+ \alpha \Delta _{k-1}^{(\pm )}(y) \right) \, {\rm d}y \\ &&\nonumber \hphantom{\fl \nonumber \Delta _{k}^{(\pm )}(s) =} - \frac{\Gamma \left(-\frac{2 \ell +1}{2} \right)}{2} \Sigma _+(s) \int _{s_{k2}}^s {\rm e}^{(M+1) y} J_{-\frac{2 \ell +1}{2}}\left({\rm e}^y\right) \left( {\rm e}^{(M+1) y} \Delta _{k-2}^{(\pm )}(y)+ \alpha \Delta _{k-1}^{(\pm )}(y) \right) \, {\rm d}y . \end{eqnarray*}$$

The constants of integration are appropriately fixed to satisfy the required boundary conditions so that contributions which do not behave appropriately are excluded form Δ_k(s).

When α vanishes, the contribution of Δ₁(s) is reabsorbed in the unperturbed solution. Curiously, the contribution Δ₂(s; α = 0) is almost the same as Δ₁(s; α = 1) and that is because $V \sim \alpha {\rm e}^{(M+1) s} \mathcal {F} + {\rm e}^{2 (M+1) s} \mathcal {F}^2$, i.e. when the supersymmetric parameter vanishes, the first contribution, coming at the order $\mathcal {F}^2$, is the square of the first contribution when such parameter is nonzero, at first order in $\mathcal {F}$.

The expressions above seem cumbersome but from them we can see that it is convenient to introduce the following two quantities:

$$\begin{eqnarray*} \fl \mathcal {J}_n (y) &=& \frac{ _2F_3\left(\ell +1,\frac{n}{2}(M+1)+\frac{2\ell +1}{2};\frac{2\ell +1}{2}+1,\frac{n}{2}(M+1)+\frac{2\ell +1}{2}+1,2 (\ell +1);-{\rm e}^{2y}\right)}{2^{2 \ell +1} {\rm e}^{-2((M+1) n+2 \ell +1)y}((M+1) n+2 \ell +1)\Gamma \left(\frac{2\ell +1}{2}+1\right)^2 }, \end{eqnarray*}$$

$$\begin{eqnarray} \fl \mathcal {K}_n (y) &=& \frac{ _2F_3\left(\frac{1}{2},\frac{n}{2}(M+1);\frac{n}{2}(M+1)+1,1-\frac{2\ell +1}{2},1+\frac{2\ell +1}{2};-{\rm e}^{2y}\right)}{(M+1) n \Gamma \left(1-\frac{2\ell +1}{2}\right) \Gamma \left( 1+ \frac{2\ell +1}{2}\right)} \, {\rm e}^{2(M+1) n y} , \end{eqnarray} \tag{ 4.41 }$$

which arise as the indefinite integral forms involving $J_{\frac{2 \ell +1}{2}}\left({\rm e}^y\right)^2$ and $J_{\frac{2 \ell +1}{2}}\left({\rm e}^y\right)J_{-\frac{2 \ell +1}{2}}\left({\rm e}^y\right)$ in the expressions for Δ_n(s) above. The equivalent result for $J_{-\frac{2 \ell +1}{2}}\left({\rm e}^y\right)^2$ can be trivially obtained by just replacing ℓ → −(ℓ + 1) in the former. Taking the lower integration limit to be x_i = ∞, or s_i = ∞, and evaluating it near the origin, x = 0 or s = −∞, as this is the region we are interested in the wavefunction, we obtain the quantities

$$\begin{eqnarray} &&J_n = \frac{\Gamma \left(\frac{1}{2}-\frac{n}{2} (M+1)\right) \Gamma \left(\frac{2\ell +1}{2}+\frac{n}{2} (M+1)\right)}{2 \sqrt{\pi } \Gamma \left(1-\frac{n}{2} (M+1)\right) \Gamma \left(\frac{2\ell +1}{2}-\frac{n}{2} (M+1)+1\right)}, \end{eqnarray} \tag{ 4.42 }$$

$$\begin{eqnarray} &&K_n = \frac{\Gamma \left(\frac{n}{2} (M+1)\right) \Gamma \left(\frac{1}{2}-\frac{n}{2} (M+1)\right)}{2 \sqrt{\pi } \Gamma \left(1-\frac{2\ell +1}{2}-\frac{n}{2} (M+1)\right) \Gamma \left(1+\frac{2\ell +1}{2}-\frac{n}{2} (M+1)\right)}. \end{eqnarray} \tag{ 4.43 }$$

These are related to each other in a very simple way,

$$\begin{equation*} \frac{J_n}{K_n}= \frac{ \cos \frac{\pi }{2} \left((M+1) n-1\right) }{ \cos \frac{\pi }{2} \left((M+1) n-1+(2 \ell +1)\right) }, \end{equation*}$$

and can be easily related to the quantity $\tilde{I}_2=\tilde{H}_1$ calculated in [5],

$$\begin{eqnarray} &&J_2 = \frac{2}{\pi } \sin {M\pi } \left( \frac{2^M}{M+1} \right)^2 \tilde{H}_1, \end{eqnarray} \tag{ 4.44 }$$

with the sin  contribution arriving because we expand $\mathcal {F} = E^{-\frac{M+1}{2}}$. These constants, K_n and J_n, play an important role in the theory as they will then appear when one is evaluating the solutions at each order.

We have seen that in every order there are contributions in the direction of both Σ_±(s) and one can collect them and determine an expansion for the then called spectral determinant. This implies that, up to the Weierstrass normalization factor which fixes D(0, α), an expansion of the spectral determinant in the form

$$\begin{eqnarray} &&D(E,\alpha ) \sim 1+ \sum _{n=1}^\infty F_n \, \mathcal {F}^n, \quad \textrm {where} \quad F_n = \lim _{s \rightarrow - \infty } \frac{\Delta _n^{(+)}(s)}{\Delta _0^{(-)}(s)}. \end{eqnarray} \tag{ 4.45 }$$

In the limit we are interested in, the above expressions allow us to compute the quantity F₁, appearing as the coefficient in the first-order expansion in terms of $E^{-\frac{M+1}{2}}$ and it turns out to be given by a particularly simple expression which is proportional to the constant J₁ we have just defined,

$$\begin{eqnarray} &&F_1 = \alpha \frac{\pi }{2} J_1. \end{eqnarray} \tag{ 4.46 }$$

This object has a simple form and does not contribute to the su(2) limit because it is proportional to the supersymmetric parameter α.

Solving the second-order term is trickier as it depends on the solution Δ₁(s), which is to be integrated in the object

$$\begin{eqnarray} &&B_1(s) = \int _{s_0}^s {\rm e}^{(M+1) y} \Delta _1(y) J_{\pm \frac{2 \ell +1}{2}}\left({\rm e}^y\right) \, {\rm d}y. \end{eqnarray} \tag{ 4.47 }$$

The evaluation in this case, when lim_{s → −∞}, shows that, although having a more complicated form, again the solution can be written down as a function of the fundamental quantities J_n, or K_n. It yields that

$$\begin{eqnarray} F_2 &=& \frac{\pi }{2} J_2 + \alpha \frac{ \Gamma \left(\frac{1}{2}-\ell \right) \Gamma \left(\ell +\frac{1}{2}\right) \Gamma \left(\ell +\frac{3}{2}\right)^2 }{4 \pi ^2 (2 \ell +1)} \left(b_1+b_2\right), \end{eqnarray} \tag{ 4.48 }$$

where the remaining constants b_i, which arise from the evaluation of B₁(s) at −∞, involve functions of the kind ₅F₄,

$$\begin{eqnarray} b_1 &=& \frac{2 \Gamma (-M) \Gamma (M+\ell +1)}{M \Gamma \left(\frac{1}{2}-M\right) \Gamma (-M+\ell +1)} \, _5F_{4,1}, \\ \nonumber b_2 &=& \frac{\pi \Gamma (M+1) \Gamma (-M-\ell -1)}{4^{\ell } (M+2 \ell +2)^2 \Gamma \left(\ell +\frac{3}{2}\right)^2 \Gamma (-M-2 \ell -2) \Gamma \left(-M-\ell -\frac{1}{2}\right)} \, _5F_{4,2}, \end{eqnarray} \tag{ 4.49 }$$

and we have briefly used the quantities ₅F_{4, 1} and ₅F_{4, 2}, respectively, as

$$\begin{eqnarray*} &&\fl _5F_4\left(\frac{1}{2},-M,-\frac{M}{2},-\ell ,\ell +1;\frac{1}{2}-M,1-\frac{M}{2},-M-\ell ,-M+\ell +1;1\right), \end{eqnarray*}$$

$$\begin{eqnarray*} &&\fl _5F_4\left(M+1,\ell +1,\frac{M}{2}+\ell +1,M+\ell +\frac{3}{2},M+2 \ell +2;\ell +\frac{3}{2},\right. \\ &&\left. \frac{M}{2}+\ell +2,M+\ell +2,2 \ell +2;1\right). \end{eqnarray*}$$

The appearance of these hypergeometric functions already in $\tilde{H}_1$ makes the complications presented by the introduction of the supersymmetric term evident, since for α = 0 such contributions would appear at higher orders, starting at $\tilde{H}_2$ or G₂ as in [12]. The present results, though, do not rely on any knowledge concerning the properties of Zeta functions.

Now we are finally in a position to evaluate the first few terms in the expansion coming from (4.45), so that the logarithm expansion of the spectral determinant takes the form

$$\begin{eqnarray} &&\fl \log D(E,\alpha ) \sim \mathcal {F} {F_1} + \mathcal {F}^2 \left( {F_2}- \frac{F_1^2}{2} \right) \nonumber\\ &&+\,\mathcal {F}^3\left(F_3 + \frac{F_1^3}{3}-F_2 F_1\right) +\mathcal {F}^4\left(F_4-\frac{F_2^2}{2}-F_3F_1+F_2 F_1^2-\frac{F_1^4}{4} \right) + \cdots .\nonumber\\ \end{eqnarray} \tag{ 4.50 }$$

In the limit of vanishing α, the odd (or half integer, according to other notations) components vanish, as expected from the su(2)-symmetric theory, for which only even (or integers, alternatively) terms are present. For this situation, the first contribution is given simply by

$$\begin{eqnarray} &&F_2 = (-1)^{-(M+1)}\left( \frac{2^M}{M+1} \right)^2 \tilde{I}_2, \end{eqnarray} \tag{ 4.51 }$$

with the standard identification

$$\begin{eqnarray} &&\tilde{I}_2 = \tilde{H}_1= (M+1)^2\Gamma (-(2 M+1)) \frac{\Gamma (M+1)}{\Gamma (-M)} \frac{\Gamma \left(\frac{2 \ell +1}{2} +M+1\right)}{\Gamma \left(\frac{2 \ell +1}{2}-M\right)}, \end{eqnarray} \tag{ 4.52 }$$

in agreement with (4.32). In possession of expressions (4.45), we can keep on with (4.41) and determine higher contributions to the solution expansion and therefore extract higher nonlocal integrals of motion. Those terms however become cumbersome with simplifications when α = 0, as odd F_n disappear from (4.50).

4.1. Towards a generalization for excited states

It is possible to extend (2.4) once more, so that the excited states of the Q operator can also be accessed by the introduction of singularities in the complex plane to the potential [3, 28]. Nevertheless, it would be interesting to verify if the method currently under investigation provides any information about the integrals of motion evaluated away from the vacuum of the Q operator. It was shown in [3], for α = 0, that higher level eigenvalues of the Baxter operator can be obtained by adding the potential poles of the form $-2\frac{{\rm d}^2}{{\rm d}x^2} \sum _{k=1}^L \log \left( x^{2(M+1)} - z_k \right)$. For the first excited state, for instance,

$$\begin{equation} z_1 = \frac{(2 \ell +1)^2-4 M^2}{M}. \end{equation} \tag{ 4.53 }$$

The corresponding part introduced in [3] that needs to be added to (4.34) is also expressed as a power series in the parameter $\mathcal {F}$,

$$\begin{eqnarray} &&U_L(s) = 4(M+1)\sum _{k=1}^L \sum _{j=1}^\infty \frac{(2 j (M+1)-1)}{z_k^{j}}\, {\rm e}^{2 s (j M+j-1)} \mathcal {F}^{\frac{2 (j M+j-1)}{M+1}}, \end{eqnarray} \tag{ 4.54 }$$

but the expansion is in rational, not integer, powers—unless M = 1. Therefore, one needs to accommodate this feature in the solution series according to

$$\begin{eqnarray} &&\Sigma (s)=\Delta _0(s)+\sum _{j=1}^{\infty } \mathcal {F}^j \Delta _j(s)+\sum _{j=1}^{\infty } \mathcal {F}^{\frac{2 (j M+j-1)}{M+1}} \Upsilon _j(s), \end{eqnarray} \tag{ 4.55 }$$

where ϒ_j are the contributions due to the presence of the higher level potential (L ≠ 0), involving non-integer powers of the perturbation parameter. Here, we simplify the notation by eliminating the upper indices ^{( ± )} and present only the positive ones, recalling the fact that the negative ones are related to these by a simple change in the sign of the combination $\frac{2\ell +1}{2}$.

In principle, there are infinitely many ϒ_j because this is the expansion of the U_L potential, and investigating the behavior of this introduction, in a similar way as we have just done for L = 0, consists of finding solutions to the ϒ_j(s) in expansion (4.55). It turns out that these terms can be solved in a closed form for a generic jth coefficient,

$$\begin{eqnarray} &&\fl \Upsilon _j(s) = \Sigma _-(s) {(2 j (M+1)-1)}(M+1)(2 \ell +1)\Gamma \left(\frac{2 \ell +1}{2} \right)^2 \nonumber\\ &&\hphantom{\fl \Upsilon _j(s) =} \times\, \sum _{k=1}^{L} \frac{1}{z_k^{j}} \; \int _{s_{j1}}^s {\rm e}^{2 y (j (M+1)-1)} J_{\frac{2 \ell +1}{2} }\left({\rm e}^y\right){}^2 {\rm d}y \nonumber\\ &&\hphantom{\fl \Upsilon _j(s) =}+\, \Sigma _+(s) {(2 j (M+1)-1)}(M+1)\frac{2\pi }{\cos \left( \pi \ell \right)} \nonumber\\ &&\hphantom{\fl \Upsilon _j(s) =}\times\, \sum _{k=1}^{L} \frac{1}{z_k^{j}} \; \int _{s_{j2}}^s {\rm e}^{2 y (j (M+1)-1)} J_{\frac{2 \ell +1}{2} }({\rm e}^y )J_{-\frac{2 \ell +1}{2} }({\rm e}^y){\rm d}y . \end{eqnarray} \tag{ 4.56 }$$

Once more, in order to determine the contribution of the excited states of level L to spectral determinant expansion, we must evaluate the wavefunction in limit where we take s → −∞. In this case, the integrals, which have a similar form as those in (4.41), can be evaluated. With this expression at hand, we have determined how the introduction of the poles z_k in the ODE affects the spectral determinant expansion,

$$\begin{eqnarray} &&\fl \Upsilon _j(s) = \Gamma \left(1-\frac{2 \ell +1}{2} \right) J_{-\frac{2 \ell +1}{2}}\left({\rm e}^s\right) \frac{(2 j (M+1)-1)}{z_k^{j}}(M+1)(2 \ell +1)\Gamma \left(\frac{2 \ell +1}{2} \right)^2 \nonumber \\ &&\times \, \frac{1}{2 \sqrt{\pi }} \frac{\Gamma \left(\frac{1}{2} (3-2 j (M+1))\right)}{\Gamma \left(\frac{1}{2} (3-2 j(M+1))+\frac{1}{2}\right)} \frac{\Gamma \left(\frac{1}{2} (2 j (M+1)-3)+\ell +1\right)}{\Gamma \left(\frac{1}{2} (3-2 j (M+1))+\ell +1\right)} \sum _{k=1}^{L} \frac{1}{z_k^{j}}\nonumber\\ &&+\, \Gamma \left(1+\frac{2 \ell +1}{2} \right) J_{\frac{2 \ell +1}{2}}\left({\rm e}^s\right) \frac{(2 j (M+1)-1)}{z_k^{j}}(M+1)\frac{2\pi }{\cos \left( \pi \ell \right)} \nonumber\\ &&\times \, \frac{1}{2 \sqrt{\pi }} \frac{\Gamma \left(\frac{1}{2} (3-2 j (M+1))\right)}{\Gamma \left(\frac{1}{2} (3-2 j (M+1))-\ell \right)} \frac{\Gamma \left(\frac{1}{2} (2 j (M+1)-2)\right)}{\Gamma \left(\frac{1}{2} (3-2 j (M+1))+\ell +1\right)} \sum _{k=1}^{L} \frac{1}{z_k^{j}} ,\nonumber\\ \end{eqnarray} \tag{ 4.57 }$$

evaluated at s = −∞. Therefore, the expansion suggested allows us to know the exact form of each of the higher level coefficients in the expansion. The complete evaluation of the function above requires, though, the knowledge about the exact positions of the poles, which is known exactly for L = 1. However, even if we include the expression for z₁, the overall contribution of the ϒ_j does not appear in D(E, α) in the precise form one is interested, namely the integer powers of $\mathcal {F}$. This prevents us from determining the effect of excited states on the nonlocal integrals of motion, indicating that the above parameterization is not good to investigate the effects of higher levels L ⩾ 0.

The situation is different when M = 1, with the introduction of poles contributing to lower orders according to

$$\begin{eqnarray} &&U_L(s) = \sum _{k=1}^L \frac{24}{z_k} \, {\rm e}^{2s} \, \mathcal {F}+ \sum _{k=1}^L \frac{56}{z_k^2} \, {{\rm e}}^{6s} \, \mathcal {F}^3+ \cdots , \end{eqnarray} \tag{ 4.58 }$$

so that the lowest contribution comes simply as a shift in the term of order $\mathcal {F}$ in (4.34). In this situation, it is enough to make the following replacement in F₁ and F₂, (4.46) and (4.48), respectively:

$$\begin{eqnarray} &&\alpha \longrightarrow \sum _{k=1}^L \frac{24}{z_k}, \end{eqnarray} \tag{ 4.59 }$$

which for the first excited state can be determined analytically with (4.53). If the location of the poles z_k(α) is known for when α ≠ 0, then in order to evaluate the first contribution of the excited states, one needs to take $\alpha \longrightarrow \alpha + \sum _{k=1}^L \frac{24}{z_k(\alpha )}$ instead.

The fact that the poles are located away from the origin guarantees a well-behaved modification to the supersymmetric parameter. Higher contributions can also be calculated, but to the best of our knowledge these objects have not been explicitly presented for a direct comparison. In fact, the study of excited states deserves separate attention. The appearance of the new terms to the original potential at orders $\mathcal {F}^3, \mathcal {F}^5, \ldots$ does not allow us to extract the higher level contributions from the vacuum in a simple way. The solution for these cases is obtained by identifying ϒ_j(s) with Δ_j(s) and following the same procedure for the modified equation at each order.

The spectral equivalence between su(2|1) supersymmetric spectral determinants D(λ) associated with a class of ODEs and the eigenvalues of the Baxter $\hat{Q}(\lambda )$-operator of the quantum AKNS field theory has been investigated. We have presented a framework capable of proving us with the nonlocal conserved charges associated with the quantum model, $\tilde{H}_n$ in order to determine the whole NLIE expansion together with the local conserved quantities I_n. This constitutes an alternative and efficient method which allows one to determine both odd and even charges, which are missed out from the spectral Zeta function approach of [12], and extend the results in [27, 8] for a supersymmetric potential in an elegant way which can be reduced to the previous case when the deformation parameter, responsible for taking an su(2)-symmetric problem into an su(2|1) analogue, is set to vanish.

Many questions remain open, especially regarding the understanding of the higher level ODE/IM correspondence for algebras beyond the local integrals of motion associated with su(2) and su(2|1) [3, 28]. Here, the first step is taken in order to also accommodate the nonlocal conserved quantities. The introduction of poles in the complex plane may in some cases destroy the reality of the spectrum for a certain ODE. Therefore, the excitation of higher $\hat{Q}(\lambda )$-states may be seen as a mechanism of breaking of the $\mathcal {PT}$-symmetry for a non-Hermitian Hamiltonian, a phenomenon which deserves attention on its own. Once we have established this equivalence between these systems also for the excited states, information extracted on one hand may be valuable on the other.

I would like to thank Clare Dunning, Joe Watkins and Giota Adamopoulou for useful discussions as well as EPSRC (EP/G039526/1) and SISSA/INFN-Trieste for financial support.