text:
# 1 Introduction

## 1 Introduction

The main subject of this paper will be the spectrum of Schrödinger equation

$$(\alpha ,l)\psi (x)=\left[\frac{d^2}{dx^2}+x^6+\alpha x^2+\frac{l(l+1)}{x^2}\right]\psi (x)=E\psi (x),$$
(1.1)
with boundary conditions $`\psi (x)0`$ as $`x\mathrm{}`$ and $`\psi (x)x^{l+1}`$ as $`x0`$.*** The other natural boundary condition at the origin is $`\psi (x)x^l`$; when we are interested in this spectral problem, we shall write the differential operator as $`(\alpha ,1l)`$.

Over the years, many interesting features of this problem have been uncovered. For $`\alpha <0`$ and $`l(l+1)=0`$, the theory corresponds to a double-well potential on the full real line, with $`l=1`$ and $`l=0`$ selecting the even and odd wavefunctions respectively. Such potentials have long been studied as toy models for instanton effects in quantum field theories . Furthermore, at $`\{\alpha =3,l=1\}`$ the ground-state energy $`E_0`$ is zero, with the remaining energy levels matching those at $`\{\alpha =3,l=0\}`$. This reflects the fact that for these values of the parameters the model provides one of the simplest examples of supersymmetric quantum mechanics . More peculiar are the properties of the spectrum at $`\alpha =(4J+2l+1)`$: a finite number of energy levels can be exactly computed  being roots of particular polynomials . At these points, the model is said to be ‘quasi-exactly solvable’. Finally, at least a couple of physically interesting spectral problems are related to the above equation via simple variable and gauge transformations. The first is a linear combination of the harmonic and Coulomb potentials:

$$\left[\frac{d^2}{dx^2}+x^2\frac{\sigma }{x}+\frac{\gamma (\gamma +1)}{x^2}\right]\phi (x)=\mathrm{\Lambda }\phi (x),$$
(1.2)
with the two sets of parameters being related by $`\gamma =(l1/2)/2,\mathrm{\Lambda }=\alpha /2,\sigma =2^{3/2}E`$. The second theory is

$$\left[\frac{d^2}{dx^2}+\frac{1}{x^{3/2}}+\frac{\beta }{x}+\frac{\delta (\delta +1)}{x^2}\right]\chi (x)=\mathrm{\Delta }\chi (x),$$
(1.3)
with $`\beta =16\alpha E^2,\delta =(l3/2)/4,\mathrm{\Delta }=4096E^4`$. Surprisingly, at $`\beta =0`$ and $`\delta (\delta +1)=0`$ the latter equation is related to the Odderon problem in QCD .

In this paper we discuss some further exact spectral relationships, relating the problems (1.1) at different values of $`\alpha `$ and $`l`$, and linking them to certain third-order differential equations. We shall find these equivalences in the framework of the ‘ODE/IM correspondence’ . Some features of this correspondence are reviewed in sections 2 and 3, and they are then used to establish five spectral equivalences in sections 4, 5, 6 and 7. Section 8 provides an alternative insight into the second and third equivalences using the so-called Bender-Dunne polynomials, and in the process uncovers differential operators which generalise the action of the supersymmetry generators to all quasi-exactly solvable points. The conclusions are in section 9, and there are two relatively self-contained appendices. The first explains a simple but efficient numerical approach to Schrödinger problems with polynomial potentials, and the second uses ideas associated with the ODE/IM correspondence to prove a generalised version of a conjecture due to Bessis, Zinn-Justin, Bender and Boettcher.

## 2 Bethe ansatz equations for the $`x^{2M}+x^{M1}`$ potential

The rôle of functional relations in the spectral theory of the Schrödinger equation has been extensively explored by Voros , but only recently has it been realised that they can lead, in certain cases, to a precise connection with the theory of integrable models . This so-called ‘ODE/IM correspondence’ has been developed in a number of directions , some of which are reviewed in . In this section, we summarise the results obtained in this context by Suzuki  concerning the Schrödinger equation with potential $`x^{2M}+\alpha x^{M1}`$, which includes the sextic potential (1.1) as a special case. Our treatment is perhaps a marginal simplification of that given in , but the only genuinely new contribution to the discussion is the inclusion of the angular momentum term. The ODE to be considered is

$$(M,\alpha ,l)\mathrm{\Phi }(x)=\left[\frac{d^2}{dx^2}+x^{2M}+\alpha x^{M1}+\frac{l(l+1)}{x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x)$$
(2.1)
with $`M`$ a positive real number, which for technical reasons will sometimes be taken to be greater than $`1`$. The goal is to find the eigenvalues $`\{E_i\}`$, those $`E`$ for which (2.1) has a solution vanishing as $`x+\mathrm{}`$, and behaving as $`x^{l+1}`$ as $`x0`$. (For $`\mathrm{}el>1/2`$, the latter condition is equivalent to the demand that the usually-dominant $`x^l`$ behaviour near the origin should be absent; outside this region, the problem is best defined by analytic continuationFor a discussion, see chapter 4 of ..) The starting-point is the uniquely-determined solution $`Y(x,E,\alpha ,l)`$ which has the following asymptotic for large, positive $`x`$ :

$$Y(x,E,\alpha ,l)\frac{x^{M/2\alpha /2}}{\sqrt{2i}}\mathrm{exp}\left(\frac{x^{M+1}}{M+1}\right),$$
(2.2)
and an associated set of functions $`Y_k`$ :

$$Y_k(x,E,\alpha ,l)=\mathrm{\Omega }^{k/2+k\alpha /2}Y(\mathrm{\Omega }^kx,\mathrm{\Omega }^{2Mk}E,\mathrm{\Omega }^{(M+1)k}\alpha ,l),\mathrm{\Omega }=\mathrm{exp}(\frac{i\pi }{M+1}).$$
(2.3)
For integer $`k`$ the $`Y_k`$’s are solutions of (2.1), and any pair $`\{Y_k,Y_{k+1}\}`$ forms a basis of solutions. In particular, $`Y_1`$ can be written as a linear combination of $`Y_0`$ and $`Y_1`$, the precise relation being

$$T(E,\alpha ,l)Y_0(x,E,\alpha ,l)=Y_1(x,E,\alpha ,l)+Y_{\mathrm{\hspace{0.33em}1}}(x,E,\alpha ,l).$$
(2.4)
The coefficient

$$T(E,\alpha ,l)=W[Y_1,Y_1]=\text{Det}\left[\begin{array}{cc}Y_1(x)& Y_1(x)\\ Y_1^{}(x)& Y_1^{}(x)\end{array}\right],$$
(2.5)
given here as a Wronskian, is called a Stokes multiplier. In terms of the original function $`Y(x,E,\alpha ,l)`$, (2.4) taken at $`\alpha `$ and at $`\alpha `$ leads to the following pair of equations:

$$T^{(+)}(E)Y^{(+)}(x,E)=\mathrm{\Omega }^{(1+\alpha )/2}Y^{()}(\mathrm{\Omega }x,\mathrm{\Omega }^{2M}E)+\mathrm{\Omega }^{(1+\alpha )/2}Y^{()}(\mathrm{\Omega }^1x,\mathrm{\Omega }^{2M}E)$$
(2.6)
$$T^{()}(E)Y^{()}(x,E)=\mathrm{\Omega }^{(1\alpha )/2}Y^{(+)}(\mathrm{\Omega }x,\mathrm{\Omega }^{2M}E)+\mathrm{\Omega }^{(1\alpha )/2}Y^{(+)}(\mathrm{\Omega }^1x,\mathrm{\Omega }^{2M}E)$$
(2.7)
where

$$T^{(\pm )}(E)=T(E,\pm \alpha ,l),Y^{(\pm )}(x,E)=Y(x,E,\pm \alpha ,l).$$
(2.8)
Keeping $`\mathrm{}el>1/2`$, the leading behaviour of $`Y`$ near the origin at generic $`E`$ is

$$Y(x,E,\alpha ,l)D(E,\alpha ,l)x^l+\mathrm{}.$$
(2.9)
At the zeroes of $`D(E)`$, the leading term is instead proportional to $`x^{l+1}`$, and $`Y(x,E,\alpha ,l)`$ is an eigenfunction of (2.1). This implies that $`D(E)`$ is proportional to the spectral determinant for the problem. (For $`M>1`$ the order of $`D(E)`$ can be shown to be less than one, so $`D(E)`$ is fixed up to a constant by the positions of its zeroes.) Setting $`D^{(\pm )}(E)=D(E,\pm \alpha ,l)`$, from (2.6) and (2.7) we have

$`T^{(+)}(E)D^{(+)}(E)`$ $`=`$ $`\mathrm{\Omega }^{(2l+1+\alpha )/2}D^{()}(\mathrm{\Omega }^{2M}E)+\mathrm{\Omega }^{(2l+1+\alpha )/2}D^{()}(\mathrm{\Omega }^{2M}E),`$ (2.10)
$`T^{()}(E)D^{()}(E)`$ $`=`$ $`\mathrm{\Omega }^{(2l+1\alpha )/2}D^{(+)}(\mathrm{\Omega }^{2M}E)+\mathrm{\Omega }^{(2l+1\alpha )/2}D^{(+)}(\mathrm{\Omega }^{2M}E).`$ (2.11)
Now let the zeroes of $`D^{(\pm )}(E)`$ be at $`\{E_k^{(\pm )}\}`$, and set $`E=E_k^{(\pm )}`$ in either (2.10) or (2.11). Both $`T^{(\pm )}(E)`$ and $`D^{(\pm )}(E)`$ are entire in $`E`$, so the LHS of the relevant equation vanishes. Factorising the functions $`D^{(\pm )}(E)`$ as products over their zeroes, for $`M>1`$ the following system of ‘Bethe ansatz’ type equations for the energy spectrum is obtained:

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{E_n^{()}\mathrm{\Omega }^{2M}E_k^{(+)}}{E_n^{()}\mathrm{\Omega }^{\mathrm{\hspace{0.33em}2}M}E_k^{(+)}}}\right)`$ $`=`$ $`\mathrm{\Omega }^{2l1\alpha };`$ (2.12)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{E_n^{(+)}\mathrm{\Omega }^{2M}E_k^{()}}{E_n^{(+)}\mathrm{\Omega }^{\mathrm{\hspace{0.33em}2}M}E_k^{()}}}\right)`$ $`=`$ $`\mathrm{\Omega }^{2l1+\alpha }.`$ (2.13)
Note that the spectra of the Hamiltonians $`(M,\alpha ,l)`$ and $`(M,\alpha ,l)`$ are completely tangled up by the Bethe ansatz constraints.

One of the products of the ODE/IM correspondence was the realisation that energy levels for Schrödinger problems can be calculated using nonlinear integral equations . For the ODE (2.1) with $`l=0`$ and $`\alpha `$ small, these were derived in , and it is straightforward to include the effect of the angular momentum term. Suppose that
(a) all the zeroes of $`D^{(\pm )}(E)`$ lie on the positive real axis of the complex $`E`$ plane, and
(b) all the zeroes of $`T^{(\pm )}(E)`$ lie away from it.
(As shown in appendix B, (a) holds if $`l>1/2`$, and (b) if $`|\alpha |<M+1|2l+1|`$. ) Setting

$$a^\pm (\theta )=\mathrm{\Omega }^{2l+1\pm \alpha }\frac{D^{()}(\mathrm{\Omega }^{2M}E)}{D^{()}(\mathrm{\Omega }^{2M}E)}\text{with}E=\mathrm{exp}\left(\frac{2M}{M+1}\theta \right),$$
(2.14)
the equations are then

$`\mathrm{ln}a^\pm (\theta )`$ $`=`$ $`i{\displaystyle \frac{\pi }{2}}(2l+1\pm {\displaystyle \frac{\alpha }{M}})ib_0e^\theta `$ (2.15)
$`+{\displaystyle _{𝒞_1}}𝑑\theta ^{}K_1(\theta \theta ^{})\mathrm{ln}(1+a^\pm (\theta ^{})){\displaystyle _{𝒞_2}}𝑑\theta ^{}K_1(\theta \theta ^{})\mathrm{ln}(1+{\displaystyle \frac{1}{a^\pm (\theta ^{})}})`$
$`+{\displaystyle _{𝒞_1}}𝑑\theta ^{}K_2(\theta \theta ^{})\mathrm{ln}(1+a^{}(\theta ^{})){\displaystyle _{𝒞_2}}𝑑\theta ^{}K_2(\theta \theta ^{})\mathrm{ln}(1+{\displaystyle \frac{1}{a^{}(\theta ^{})}})`$
The integration contours $`𝒞_1`$ and $`𝒞_2`$ run just below and just above the real axis respectively, and the constant $`b_0=\pi ^{1/2}\mathrm{\Gamma }(\frac{1}{2M})/(2M\mathrm{\Gamma }(\frac{3}{2}+\frac{1}{2M}))`$ is fixed via a consideration of the WKB asymptotics of $`D^{(\pm )}(E)`$ for $`|E|\mathrm{}`$, $`\mathrm{arg}(E)0`$. An integral expression for the kernels $`K_1`$ and $`K_2`$ at general values of $`M`$ was found in . In this paper we are particularly interested in the sextic potential, and here we note that the special nature of this case is reflected in the fact that at $`M=3`$ the kernel functions can be explicitly integrated, and have the simple forms

$$K_1(\theta )=\frac{\sqrt{3}}{2\pi (2\mathrm{cosh}\theta +1)},K_2(\theta )=\frac{\sqrt{3}}{2\pi (2\mathrm{cosh}\theta 1)}.$$
(2.16)
A search along the real $`\theta `$ axis for the zeroes of the functions $`1+a^\pm (\theta )`$ provides the energy levels of the Hamiltonians $`(M,\pm \alpha ,l)`$. While we do not have a rigorous proof, we expect that the solution to (2.15), which can readily be obtained numerically by iteration, is unique for $`\alpha `$ and $`l`$ in the stated range. This is one way to justify the claim that, for $`l>1/2`$ and $`|\alpha |<M+1|2l1|`$, the Bethe Ansatz equations (2.12) and (2.13), together with the ‘analytic properties’ (a) and (b) and the WKB asymptotic which determined $`b_0`$, characterise the set of numbers $`\{E_k^{(\pm )}\}`$ uniquely. (An alternative approach to this question might be to generalise the analysis of  and appendix A of , based on the so-called quantum Wronskian relations.)

To treat more general values of $`\alpha `$ and $`l`$, we will find it most convenient to work directly with the Bethe ansatz equations. In contrast to the integral equation, these do not have a unique solution. As is standard in studies of the Bethe ansatz, it is useful to take logarithms. For $`l>1/2`$ and $`|\alpha |<M+2l+2`$ (see appendix B), all energies $`E_k^{(\pm )}=E_0^{(\pm )},E_1^{(\pm )},E_2^{(\pm )}\mathrm{}`$ are positive. For this situation we assume that there is a one-to-one correspondence between these energies and the integers $`k=0,1,2,\mathrm{}`$ :

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{()}\mathrm{\Omega }^{2M}E_k^{(+)}}{E_n^{()}\mathrm{\Omega }^{2M}E_k^{(+)}}}\right)`$ $`=`$ $`i\pi \left[{\displaystyle \frac{2l+1+\alpha }{M+1}}+2k+1\right];`$ (2.17)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{(+)}\mathrm{\Omega }^{2M}E_k^{()}}{E_n^{(+)}\mathrm{\Omega }^{2M}E_k^{()}}}\right)`$ $`=`$ $`i\pi \left[{\displaystyle \frac{2l+1\alpha }{M+1}}+2k+1\right],`$ (2.18)
where the logarithms are all on the principal branch: $`\pi i\mathrm{ln}<\pi `$. At larger values of $`|\alpha |`$ and $`|l|`$ some of the low-lying energies might become negative, and in such cases care must be taken to keep track of the nontrivial monodromy of the log function.

## 3 The Bethe ansatz approach to a third-order equation

This section summarises the derivation of equations of Bethe ansatz type for a third-order ODE with ‘potential’ $`x^{3N}`$, following . We start with the equation

$$\left[\frac{d^3}{dx^3}+x^{3N}+\frac{L}{x^3}G\left(\frac{1}{x^2}\frac{d}{dx}\frac{1}{x^3}\right)\right]\varphi (x)=E\varphi (x),$$
(3.1)
and, as in , rewrite it as

$$\left[\mathrm{D}(𝐠)+x^{3N}\right]\varphi (x)=E\varphi (x)$$
(3.2)
where $`𝐠=\{g_0,g_1,g_2\}`$ with $`g_0+g_1+g_2=3`$, and

$$\mathrm{D}(𝐠)=\mathrm{D}(g_22)\mathrm{D}(g_11)\mathrm{D}(g_0),\mathrm{D}(g)=\left(\frac{d}{dx}\frac{g}{x}\right).$$
(3.3)
The relationship between $`𝐠`$ and $`\{G,L\}`$ is

$$G=g_0g_1+g_0g_2+g_1g_22,L=2g_0g_1g_2(g_0g_1+g_0g_2+g_1g_2).$$
(3.4)
Again we introduce a uniquely-defined function $`y(x,E,𝐠)`$, which solves (3.1) and tends to zero as $`x\mathrm{}`$ along the positive real axis as

$$y(x,E,𝐠)\frac{x^N}{i\sqrt{3}}\mathrm{exp}\left(\frac{x^{N+1}}{N+1}\right).$$
(3.5)
Given $`y(x)`$, bases of solutions are constructed just as in the second-order case. Set

$$y_k(x,E,𝐠)=\omega ^ky(\omega ^kx,\omega ^{3Nk}E,𝐠),\omega =\mathrm{exp}\left(\frac{2\pi i}{3N+3}\right).$$
(3.6)
For integer $`k`$, $`y_k`$ solves (3.1), and $`\{y_k,y_{k+1},y_{k+2}\}`$ form a basis. We can therefore expand $`y_0`$ as

$$y_0C^{(1)}y_1+C^{(2)}y_2y_3=0$$
(3.7)
with coefficients $`C^{(1)}`$ and $`C^{(2)}`$ – Stokes multipliers – which are independent of $`x`$. Eliminating $`C^{(2)}`$, we have

$$W_{02}C^{(1)}W_{12}+W_{23}=0,$$
(3.8)
where the Wronskians of pairs of solutions,

$$W_{k_1k_2}=W[y_{k_1}y_{k_2}]=\text{Det}\left[\begin{array}{cc}y_{k_1}(x)& y_{k_2}(x)\\ y_{k_1}^{}(x)& y_{k_2}^{}(x)\end{array}\right],$$
(3.9)
were used. (Note, since $`y_{k_1}`$ and $`y_{k_2}`$ solve a third-order equation, these are nontrivial functions of $`x`$.) Now multiply by $`y_1`$ and use the relation $`y_1W_{02}=y_0W_{12}+y_2W_{01}`$ to find

$$C^{(1)}y_1W_{12}=y_0W_{12}+y_2W_{01}+y_1W_{23}.$$
(3.10)
In this form the relation can be rewritten in terms of just two functions, $`y(x,E,𝐠)`$ and $`W(x,E,𝐠)=W_{01}(x,E,𝐠)`$ :

$`C^{(1)}(E)y(\omega ^1x,\omega ^{3N}E)W(\omega ^1x,\omega ^{3N}E)=\omega ^1y(x,E)W(\omega ^1x,\omega ^{3N}E)`$ $`+`$
$`y(\omega ^2x,\omega ^{6N}E)W(x,E)+\omega y(\omega ^1x,\omega ^{3N}E)W(\omega ^2x,\omega ^{6N}E).`$ (3.11)
We initially suppose that $`\mathrm{}e(g_0)<\mathrm{}e(g_1)<\mathrm{}e(g_2)`$. Then as $`x0`$ the leading behaviours of $`y`$ and $`W`$ are

$$y(x)D^{(1)}(E,𝐠)x^{g_0},W(x)D^{(2)}(E,𝐠)x^{g_0+g_11}.$$
(3.12)
Using equations (3.12) the relation (3.11) becomes

$`C^{(1)}(E)D^{(1)}(\omega ^{3N}E)D^{(2)}(\omega ^{3N}E)=\omega ^{g_01}D^{(1)}(E)D^{(2)}(\omega ^{3N}E)+`$
$`\omega ^{g_11}D^{(1)}(\omega ^{6N}E)D^{(2)}(E)+\omega ^{2g_0g_1}D^{(1)}(\omega ^{3N}E)D^{(2)}(\omega ^{6N}E).`$ (3.13)
Again we shall consider this expression at the zeroes of $`D^{(1)}`$ and $`D^{(2)}`$. It is convenient to write these as $`E_k^{(1)}`$ and $`\omega ^{3N/2}E_k^{(2)}`$ , so that

$$D^{(1)}(E_k^{(1)},𝐠)=0,D^{(2)}(\omega ^{3N/2}E_k^{(2)},𝐠)=0.$$
(3.14)
As in the second-order case, these functions have a spectral interpretation. In particular, the vanishing of $`D^{(1)}(E)`$ signals the existence of a solution, at that value of $`E`$, to (3.2), decaying as $`x\mathrm{}`$, and having a faster-than-usual decay at the origin:

$$y(x)x^{\mathrm{min}(g_1,g_2)},x0.$$
(3.15)
(Since $`g_1`$ and $`g_2`$ can be complex, by $`\mathrm{min}(g_1,g_2)`$ we mean whichever of $`g_1`$ and $`g_2`$ has the smallest real part.) Evaluating (3.13) at $`E\{E_k^{(1)}\}`$ and $`E\{\omega ^{3N/2}E_k^{(2)}\}`$ and imposing the entirety of $`C^{(1)}(E)`$ leads to the following set of $`SU(3)`$-related BA equations, with $`k=0,1,2,\mathrm{}`$ :

$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{E_j^{(1)}\omega ^{3N}E_k^{(1)}}{E_j^{(1)}\omega ^{3N}E_k^{(1)}}}\right)\left({\displaystyle \frac{E_j^{(2)}\omega ^{\frac{3N}{2}}E_k^{(1)}}{E_j^{(2)}\omega ^{\frac{3N}{2}}E_k^{(1)}}}\right)`$ $`=`$ $`\omega ^{g_0g_1};`$ (3.16)
$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{E_j^{(2)}\omega ^{3N}E_k^{(2)}}{E_j^{(2)}\omega ^{3N}E_k^{(2)}}}\right)\left({\displaystyle \frac{E_j^{(1)}\omega ^{\frac{3N}{2}}E_k^{(2)}}{E_j^{(1)}\omega ^{\frac{3N}{2}}E_k^{(2)}}}\right)`$ $`=`$ $`\omega ^{2g_1+g_03}.`$ (3.17)
Since $`g_0+g_1+g_2=3`$, the right-hand sides of equations (3.16) and (3.17) can be given a more symmetrical appearance by rewriting them as $`\omega ^{2g_0+g_23}`$ and $`\omega ^{2g_2g_0+3}`$ respectively. These equations, together with WKB-like asymptotics for $`D^{(1)}(E)`$ and $`D^{(2)}(E)`$, fix the numbers $`E_k^{(1)}`$ and $`E_k^{(2)}`$ up to discrete ambiguities, which for the problems in hand can be eliminated given some facts about the approximate positions of the zeroes of the functions $`D^{(1,2)}(E)`$ and some associated functions $`a^{(1,2)}(E)`$. These are analogous to the analyticity conditions $`(𝐚)`$ and $`(𝐛)`$ of the previous section, and are described in more detail in . It is also possible to solve the system via a nonlinear integral equation, but this will not be needed here.

## 4 The first spectral equivalence

The first spectral equivalence follows from observing that at $`N=1`$, $`\omega ^{3N}=1`$ and $`\omega ^{\frac{3N}{2}}=i`$. The $`SU(3)`$ BA equations therefore simplify to

$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{E_j^{(2)}iE_k^{(1)}}{E_j^{(2)}+iE_k^{(1)}}}\right)`$ $`=`$ $`\omega ^{2g_0+g_23};`$ (4.1)
$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{E_j^{(1)}iE_k^{(2)}}{E_j^{(1)}+iE_k^{(2)}}}\right)`$ $`=`$ $`\omega ^{2g_2g_0+3}.`$ (4.2)
These equations coincide with the system (2.12), (2.13) at $`M=3`$ provided the right-hand sides of the two BA sets are equated:

$$(2g_0+g_23)/3=(2l1\alpha )/4,(2g_2+g_03)/3=(2l1+\alpha )/4.$$
(4.3)
Combined with a matching of the analytic properties (a) and (b), this suggests the following relationship between quantities in the two problems:

$`D^{(1)}(\kappa ^1E,𝐠)`$ $`=`$ $`f(\alpha ,l)D^{(+)}(E,\alpha ,l),`$ (4.4)
$`D^{(2)}(i\kappa ^1E,𝐠)`$ $`=`$ $`f(\alpha ,l)D^{()}(E,\alpha ,l).`$ (4.5)
The proportionality factors $`f(\alpha ,l)`$ and $`\kappa `$ cannot be determined by a comparison of the Bethe ansatz equations alone. However, as in , $`\kappa `$ can be calculated by comparing the large negative $`E`$ asymptotics of $`D^{(1)}`$ and $`D^{(2)}`$. The result, independent of $`\alpha `$ and $`l`$, is $`\kappa =4/(3\sqrt{3})`$. Solving (4.3), the parameters $`(\alpha ,l)`$ and $`𝐠`$ are related as

$$\alpha \alpha (g_0,g_2)=2(2g_0g_2),ll(g_0,g_2)=(2g_232g_0)/6,$$
(4.6)
and

$$g_0=(1\alpha 6l)/4,g_1=(1+\alpha /2),g_2=(7\alpha +6l)/4.$$
(4.7)
Thus we have a spectral equivalence between the following eigenvalue problems

$$\left[\frac{d^2}{dx^2}+x^6+\alpha x^2+\frac{l(l+1)}{x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x),\mathrm{\Phi }|_{x0}x^{l+1};$$
(4.8)
$$\kappa \left[\frac{d^3}{dx^3}+x^3+\frac{L}{x^3}G\left(\frac{1}{x^2}\frac{d}{dx}\frac{1}{x^3}\right)\right]\varphi (x)=E\varphi (x),\varphi |_{x0}x^{\mathrm{min}(g_1,g_2)},$$
(4.9)
where

$$G=g_0g_1+g_0g_2+g_1g_22,L=2g_0g_1g_2(g_0g_1+g_0g_2+g_1g_2),$$
(4.10)
and the parameters in the two models are related as in (4.6) and (4.7). Note that the general $`SU(3)`$-related equation at $`N=1`$ is mapped onto the general sextic-potential problem. The number of parameters matches up, because the third-order equation allows for two linearly-independent angular momentum type terms. The different rôles that these parameters play in the second-order problem will be important for the next spectral equivalence that we discuss.

## 5 The second spectral equivalence

Our second spectral equivalence is related to the enhanced symmetries of the third-order equation. To be more precise, the differential equation (3.2) is unchanged under permutations of $`\{g_0,g_1,g_2\}`$, while the values of $`\alpha `$ and $`l`$ which appear in the corresponding second-order equation, given by (4.6), are not. If we make a continuation in $`\{g_0,g_1,g_2\}`$ which swaps $`g_1`$ and $`g_2`$ while leaving $`g_0`$ unchanged, then both the third-order equation itself, and the specification (3.12) of $`D^{(1)}`$, are unchanged, and so

$$D^{(1)}(\kappa ^1E,𝐠)D^{(1)}(\kappa ^1E,𝐠).$$
(5.1)
Using (4.4), this means that

$$D^{(+)}(E,\stackrel{~}{\alpha },\stackrel{~}{l})=\frac{f(\alpha ,l)}{f(\stackrel{~}{\alpha },\stackrel{~}{l})}D^{(+)}(E,\alpha ,l),$$
(5.2)
with

$$\stackrel{~}{\alpha }\alpha (g_0,g_1)=(3\alpha +6l)/2,\stackrel{~}{l}l(g_0,g_1)=(\alpha +2l1)/4.$$
(5.3)
It will sometimes be convenient to put this in matrix form. If $`\alpha \alpha =(\alpha ,l,1)^T`$, then

$$\stackrel{~}{\alpha }\alpha =𝕋\alpha \alpha ,\text{with}𝕋=\left(\begin{array}{ccc}1/2& 3& 3/2\\ 1/4& 1/2& 1/4\\ 0& 0& 1\end{array}\right).$$
(5.4)
At this stage we do not know how to calculate $`f(\alpha ,l)`$ and $`f(\stackrel{~}{\alpha },\stackrel{~}{l})`$ exactly, but (5.2) can be combined with (B.15) from appendix B to give their ratio:

$$\frac{f(\alpha ,l)}{f(\stackrel{~}{\alpha },\stackrel{~}{l})}=\frac{\mathrm{\Gamma }\left(\stackrel{~}{l}+\frac{1}{2}\right)}{\mathrm{\Gamma }\left(l+\frac{1}{2}\right)}.$$
(5.5)
Note that this is singular or zero at negative-half-integer values of $`\stackrel{~}{l}`$ or $`l`$, at which a ‘resonance’ is expected in one or other of the spectral problems . Away from these points, we have a spectral equivalence between

$$\left[\frac{d^2}{dx^2}+x^6+\alpha x^2+\frac{l(l+1)}{x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x),\mathrm{\Phi }|_{x0}x^{l+1},$$
(5.6)
and

$$\left[\frac{d^2}{dx^2}+x^6+\frac{(3\alpha +6l)}{2}x^2+\frac{(\alpha +2l1)(\alpha +2l+3)}{16x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x),$$
(5.7)
with $`\mathrm{\Phi }|_{x0}x^{(\alpha +2l1)/4+1}`$. An alternative viewpoint on this equivalence in terms of intertwining operators will be given in §8 below, while some direct numerical checks are reported in appendix A.

## 6 The third spectral equivalence

As was mentioned in the introduction, at the special values $`\alpha =\alpha _J(l)=(4J+2l+1)`$, with $`J`$ a positive integer, the model (1.1) is ‘quasi-exactly solvable’ (QES), and the first $`J`$ energy levels can be computed exactly. For $`J=1`$, the single exactly-solvable energy is the ground state and the model is an example of supersymmetric quantum mechanics. This is signalled by the fact that the Hamiltonian at $`\alpha =\alpha _1(l)=(2l+5)`$ can factorised in terms of first-order operators as

$$(\alpha _1,l)\left[\frac{d^2}{dx^2}+x^6(2l+5)x^2+\frac{l(l+1)}{x^2}\right]=𝒬^{}𝒬^+,$$
(6.1)
where

$$𝒬^{}=\left[\frac{d}{dx}+x^3\frac{l+1}{x}\right],𝒬^+=\left[\frac{d}{dx}+x^3\frac{l+1}{x}\right].$$
(6.2)
The SUSY ‘partner’ Hamiltonian $`\widehat{}=(\widehat{\alpha }_1,\widehat{l}_1)`$ is obtained through the intertwining relation $`(\widehat{\alpha }_1,\widehat{l}_1)𝒬^+=𝒬^+(\alpha _1,l)`$ with $`\widehat{\alpha }_1(l)=12l`$ and $`\widehat{l}_1(l)=l+1`$ :

$$\widehat{}=𝒬^+𝒬^{}=\left[\frac{d^2}{dx^2}+x^6+(12l)x^2+\frac{(l+1)(l+2)}{x^2}\right].$$
(6.3)
The wavefunctions of the two models are simply related by

$$\widehat{\psi }_i(x)=𝒬^+\psi _i(x).$$
(6.4)
However, the ground-state wavefunction of $`(\alpha ,l)`$ is $`\psi _0=x^{l+1}\mathrm{exp}(\frac{x^4}{4}),`$ and this is annihilated by $`𝒬^+`$. As a result, $``$ and $`\widehat{}`$ are spectrally equivalent save for the extra level at $`E=0`$ only present in $`\text{Spec}()`$. This (very standard) result makes it natural to ask whether similar ‘partner potentials’ might exist at higher values of $`J`$, sharing the same spectra as the QES Hamiltonians $`(\alpha _J(l),l)`$ apart from the first $`J`$ levels. We shall find that this question has a surprisingly simple answer by using the Bethe ansatz approach to the spectral problem.

Setting $`\alpha =\alpha _J(l)=(4J+2l+1)`$, the $`J`$ exactly-solvable levels $`E_0^{(+)}`$, $`E_1^{(+)}`$, …$`E_{J1}^{(+)}`$ lie in the sector ‘$`(+)`$’. The BAE (2.17), (2.18) for $`M=3`$ are then

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{()}iE_k^{(+)}}{E_n^{()}+iE_k^{(+)}}}\right)`$ $`=`$ $`i\pi [J+2k+1];`$ (6.5)
$`{\displaystyle \underset{n=0}{\overset{J1}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{(+)}iE_k^{()}}{E_n^{(+)}+iE_k^{()}}}\right)+{\displaystyle \underset{n=J}{\overset{\mathrm{}}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{(+)}iE_k^{()}}{E_n^{(+)}+iE_k^{()}}}\right)`$ $`=`$ $`i\pi [l+3/2+J+2k],`$ (6.6)
where the integer $`k`$ runs from 0 to $`\mathrm{}`$. Next, we will use the fact that the exactly-solvable energy levels appear symmetrically, as $`E_i^{(+)}=E_{Ji1}^{(+)}`$, to simplify the first sum on the LHS of (6.6). Recalling that, since $`\alpha `$ is negative, the $`E^{()}`$ are positive for $`l>1/2`$, and keeping track of the monodromy of the logarithms by using the reflection formula

$$\mathrm{ln}\left(\frac{xi}{x+i}\right)=2\pi i\mathrm{ln}\left(\frac{xi}{x+i}\right),(x0),$$
(6.7)
we obtain

$$\underset{n=0}{\overset{J1}{}}\mathrm{ln}\left(\frac{E_n^{(+)}iE_k^{()}}{E_n^{(+)}+iE_k^{()}}\right)=i\pi J.$$
(6.8)
Finally, ignoring the first $`k=0,1,\mathrm{}J1`$ instances of (6.5) and relabelling $`E_{k+J}^{(+)}E_k^{(+)}`$ we end up with

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{()}iE_k^{(+)}}{E_n^{()}+iE_k^{(+)}}}\right)`$ $`=`$ $`i\pi [J+2k+1];`$ (6.9)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left({\displaystyle \frac{E_n^{(+)}iE_k^{()}}{E_n^{(+)}+iE_k^{()}}}\right)`$ $`=`$ $`i\pi [l+2k+3/2],`$ (6.10)
where $`k`$ again runs from $`0`$ to $`\mathrm{}`$. Comparing with (2.17), (2.18) we now reinterpret the left-hand sides of equations (6.9) and (6.10) as the quantisation conditions for the energy levels of a new potential, with parameters $`\widehat{\alpha }_J`$ and $`\widehat{l}_J`$ :

$`i\pi \left[J+2k+1\right]`$ $`=`$ $`i\pi \left[{\displaystyle \frac{2\widehat{l}_J+1+\widehat{\alpha }_J}{4}}+2k+1\right],`$ (6.11)
$`i\pi \left[l+2k+3/2\right]`$ $`=`$ $`i\pi \left[{\displaystyle \frac{2\widehat{l}_J+1\widehat{\alpha }_J}{4}}+2k+1\right].`$ (6.12)
Solving,

$$\widehat{\alpha }_J=2J2l1=(3+\alpha +6l)/2,\widehat{l}_J=J+l=(\alpha +2l1)/4,$$
(6.13)
and this can again be put in matrix form, for $`\alpha \alpha =(\alpha ,l,1)^T`$, as

$$\widehat{\alpha }\alpha =\alpha \alpha ,\text{with}=\left(\begin{array}{ccc}1/2& 3& 3/2\\ 1/4& 1/2& 1/4\\ 0& 0& 1\end{array}\right).$$
(6.14)
Then, for $`J`$,

$$E_{k+J}^{(+)}(\alpha _J,l)=E_k^{(+)}(\widehat{\alpha }_J,\widehat{l}_J),(k=0,1,2,\mathrm{})$$
(6.15)
and, modulo the exactly-solvable levels, there is a spectral equivalence between

$$\left[\frac{d^2}{dx^2}+x^6(4J+2l+1)x^2+\frac{l(l+1)}{x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x),\mathrm{\Phi }|_{x0}x^{l+1}$$
(6.16)
and

$$\left[\frac{d^2}{dx^2}+x^6+(2J2l1)x^2+\frac{(J+l)(J+l+1)}{x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x),\mathrm{\Phi }|_{x0}x^{J+l+1}.$$
(6.17)
If $`J`$ is a negative integer, the mapping still makes sense, but it acts in the opposite sense: shifting $`llJ`$, (6.17) becomes the QES problem for $`\alpha _{|J|}`$, and (6.16) its partner with the QES levels removed.

Finally, we still have the freedom to apply the ‘tilde-duality’ $`𝕋`$ discussed in §5, so (6.17) is in turn isospectral to

$$\left[\frac{d^2}{dx^2}+x^6+(2J+4l+2)x^2+\frac{(J\frac{1}{2})(J+\frac{1}{2})}{x^2}\right]\mathrm{\Phi }(x)=E\mathrm{\Phi }(x),\mathrm{\Phi }|_{x0}x^{J1/2+1}.$$
(6.18)
This chain of equivalences will be discussed further in the conclusions.

Since the quasi-exactly solvable energies and the associated wavefunctions are in principle exactly known, one could eliminate them one by one using the Darboux transformation, though this would be a lengthy business for large values of $`J`$. What seems surprising about the results (6.17) and (6.18) is that the potential can have such a simple form once all of these levels have been subtracted.

## 7 The fourth and fifth spectral equivalences

The equivalence of the second-order equation with an $`SU(3)`$-related third-order equation suggests two further spectral equivalences. As explained in , the $`_2`$ symmetry of the $`SU(3)`$ Dynkin diagram is reflected in a relation between the functions $`y`$ and $`W`$ that were introduced in §3. Explicitly,

$$y(x,E,𝐠^{})=W[y_{1/2},y_{1/2}](x,E,𝐠),$$
(7.1)
where $`𝐠^{}=\{g_0^{},g_1^{},g_2^{})`$ and $`g_i^{}=2g_{2i}^{}`$ . On the second-order side of the story, a similar Wronskian appears, but this time in the formula (2.5) for $`T(E,\alpha ,l)`$:

$$T(E,\alpha ,l)=W[Y_1,Y_1](E,\alpha ,l).$$
(7.2)
Before the two equations can be compared, the $`x`$-dependence must be eliminated from (7.1), and as usual this is done by considering the behaviour as $`x0`$. Extending the definition (3.12), we define three functions $`D_{[i]}^{(1)}`$, with $`D_{}^{(1)}D_{[0]}^{(1)}`$, by

$$y(x,E,𝐠)=\underset{i=0}{\overset{2}{}}D_{[i]}^{(1)}(E,𝐠)\chi _i(x,E,𝐠),$$
(7.3)
where the solutions $`\chi _i(x,E,𝐠)`$ to (3.2) are defined by $`\chi _ix^{g_i}+O(x^{g_i+3})`$ as $`x0`$. To match (7.2), we expand out the RHS of (7.1) and then project onto the component behaving as $`x^{g_0+g_21}=x^{g_1^{}}`$ to find

$`D_{[1]}^{(1)}(E,𝐠^{})`$ $`=`$ $`(g_2g_0)[\omega ^{(g_0g_2)/2}D_{[0]}^{(1)}(\omega ^{3N/2}E,𝐠)D_{[2]}^{(1)}(\omega ^{3N/2}E,𝐠)`$ (7.4)
$`\omega ^{(g_2g_0)/2}D_{[0]}^{(1)}(\omega ^{3N/2}E,𝐠)D_{[2]}^{(1)}(\omega ^{3N/2}E,𝐠)].`$
On the other hand, the function $`Y(x,E,\alpha ,l)`$ can be expanded as

$$Y(x,E,\alpha ,l)=D(E,\alpha ,l)X(x,E,\alpha ,l)+D(E,\alpha ,1l)X(x,E,\alpha ,1l)$$
(7.5)
with $`X(x,E,\alpha ,l)x^l`$ as $`x0`$. (Cf. eq. (5.2) of , but note that the definition of $`D(E,l)`$ used in  differs from that used here by a factor of $`(2l+1)^1`$. ) Now substitute into (7.2) taken at $`(E,\alpha ,l)`$ :

$`T(E,\alpha ,l)`$ $`=`$ $`(2l+1)[\mathrm{\Omega }^{2l1}D(\mathrm{\Omega }^{2M}E,\alpha ,l)D(\mathrm{\Omega }^{2M}E,\alpha ,1l)`$ (7.6)
$`\mathrm{\Omega }^{2l+1}D(\mathrm{\Omega }^{2M}E,\alpha ,l)D(\mathrm{\Omega }^{2M}E,\alpha ,1l)].`$
If $`M=3`$ and $`N=1`$, then $`\mathrm{\Omega }^{2M}=\omega ^{3N/2}=e^{\pi i/2}`$, and, if $`𝐠`$ and $`l`$ are related by (4.6), $`\omega ^{(g_0g_2)/2}=\mathrm{\Omega }^{2l1}=e^{\pi i(2l+1)/4}`$. Furthermore, from (4.4),

$$D_{[0]}^{(1)}(E,𝐠)=f(\alpha ,l)D(\kappa E,\alpha ,l),D_{[2]}^{(1)}(E,𝐠)=f(\alpha ,1l)D(\kappa E,\alpha ,1l).$$
(7.7)
(The second relation is obtained by a continuation which swaps $`g_0`$ and $`g_2`$.) Using these identifications and comparing (7.4) and (7.6) gives our fourth spectral equivalence:

$$D_{[1]}^{(1)}(\kappa ^1E,𝐠^{})=\frac{3}{2}f(\alpha ,l)f(\alpha ,1l)T(E,\alpha ,l),$$
(7.8)
where $`g_i^{}=2g_{2i}^{}`$ and $`\alpha =2(2g_0g_2)`$, $`l=(2g_232g_0)/6`$ . As with the first equivalence, this relates spectral data for differential equations of different orders. The spectral interpretation of functions such as $`T`$ in the ODE/IM correspondence was discussed in , and is reviewed and extended to the current context in appendix B below.

We can obtain a relation between objects in the second-order equation by using the first spectral equivalence to rewrite the LHS of (7.8), and this constitutes our fifth and final spectral equivalence. The only subtlety is that (4.4) involves $`D_{[0]}^{(1)}`$, not $`D_{[1]}^{(1)}`$, and this can be overcome by a continuation in the $`g_i`$ . Swapping $`g_0`$ and $`g_1`$ and tracing back,

$$T(E,\alpha ,l)=\frac{2f(\stackrel{~}{\alpha },\stackrel{~}{l})}{3f(\alpha ,l)f(\alpha ,1l)}D(E,\stackrel{~}{\alpha }(\alpha ,l),\stackrel{~}{l}(\alpha ,l)).$$
(7.9)
The proportionality factor can also be found explicitly, by considering (7.9) at $`E=0`$ and using formulae (B.15) and (B.16). The result:

$$T(E,\alpha ,l)=\frac{2\sqrt{i\pi }}{\mathrm{\Gamma }\left(\stackrel{~}{l}(\alpha ,l)+\frac{1}{2}\right)}D(E,\stackrel{~}{\alpha }(\alpha ,l),\stackrel{~}{l}(\alpha ,l)).$$
(7.10)
Via the second equivalence, this can be rewritten as

$$T(E,\alpha ,l)=\frac{2\sqrt{i\pi }}{\mathrm{\Gamma }\left(\stackrel{~}{l}(\alpha ,l)\frac{1}{2}\right)}D(E,\stackrel{~}{\alpha }(\alpha ,l),1\stackrel{~}{l}(\alpha ,l)).$$
(7.11)
As explained in appendix B, $`T`$ is the spectral determinant for a ‘lateral connection’ problem, with the wavefunction lying on a contour in the complex plane joining a pair of Stokes sectors at infinity. In contrast, $`D`$ is the spectral determinant for a ‘radial connection’ problem, with the wavefunction living on a half-line. Similar equivalences, albeit for slightly different potentials, have been found in , and it would be interesting to see whether similar methods could be applied in this case.

The mappings of parameters involved in these relations can be streamlined by introducing two further matrices, $`𝔸`$ and $`𝕃`$, again acting on the vectors $`\alpha \alpha =(\alpha ,l,1)^T`$:

$$𝔸=\left(\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1\end{array}\right);𝕃=\left(\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 1\\ \hfill 0& \hfill 0& \hfill 1\end{array}\right),$$
(7.12)
and setting

$$\gamma (\alpha \alpha )\gamma ((\alpha ,l,1)^T)=2\sqrt{i\pi }/\mathrm{\Gamma }(l+\frac{1}{2}).$$
(7.13)
The tilde-duality (5.2) is then

$$\gamma (\alpha \alpha )D(E,\alpha \alpha )=\gamma (𝕋\alpha \alpha )D(E,𝕋\alpha \alpha ),$$
(7.14)
while (7.10) and (7.11) are, respectively,

$$T(E,𝔸\alpha \alpha )=\gamma (𝔸𝕋\alpha \alpha )D(E,𝔸𝕋\alpha \alpha ),T(E,\alpha \alpha )=\gamma (𝕃𝕋\alpha \alpha )D(E,𝕃𝕋\alpha \alpha ).$$
(7.15)
Note also that $`=𝔸𝕋𝔸`$, so the first relation of (7.15) can be rewritten as $`T(E,\alpha \alpha )=\gamma (\alpha \alpha )D(E,\alpha \alpha )`$. At the QES points, $`\alpha \alpha =\alpha \alpha _J=(\alpha _J(l),l,1)^T`$, the spectrum encoded by $`D(E,\alpha \alpha _J)`$ is equal to that of $`D(E,\alpha \alpha _J)`$, apart from the QES levels. As will be explained in the next section, these levels are in fact the zeroes of the Bender-Dunne polynomial $`P_J(E)`$, so at the QES points we have

$$P_J(E)T(E,\alpha \alpha _J)D(E,\alpha \alpha _J),$$
(7.16)
which is a relation between spectral data for lateral and radial connection problems with the same Hamiltonian.

The algebra of the matrices we have introduced is best described by first defining $`𝕄=𝔸𝕃`$. Then a set of defining relations for $`𝕃`$, $`𝕄`$ and $`𝕋`$ is

$$𝕃^2=𝕄^2=𝕋^2=(𝕃𝕄)^2=(𝕄𝕋)^2=(𝕃𝕋)^3=𝕀.$$
(7.17)
Thus $`𝕄`$ commutes with $`𝕃`$ and $`𝕋`$, while $`<𝕃,𝕋>`$ forms the Weyl group of $`SU(3)`$. However, in general only $`𝕋`$ yields a spectral equivalence of the $`D(E,\alpha \alpha )`$. We will return to this point in the conclusions.

We end this section with two further remarks about the fifth set of equivalences. First, they can also be obtained entirely in the context of the second-order differential equation. For $`M=3`$, manipulating equation (2.11) and using the fact that $`\mathrm{\Omega }^{2M}=i`$ leads to the following functional relation, special to this particular value of $`M`$:

$`2\mathrm{sin}(\frac{\pi }{4}(2l+1\alpha ))D^{(+)}(E)`$ $`=`$ (7.18)
$`\mathrm{\Omega }^{(2l+1\alpha )/2}T^{()}(iE)D^{()}(iE)\mathrm{\Omega }^{(2l+1\alpha )/2}T^{()}(iE)D^{()}(iE).`$
Taking (7.18) at $`E=E_k^{(+)}`$ , combining it with (2.11), also at $`E=E_k^{(+)}`$, and finally expressing the result in a factorised form over the zeroes of $`T^{()}(E)`$ (which we denote by $`\{\lambda _k^{()}\}`$ ) yields the following set of constraints:

$$\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{\lambda _n^{()}\mathrm{\Omega }^{2M}E_k^{(+)}}{\lambda _n^{()}\mathrm{\Omega }^{\mathrm{\hspace{0.33em}2}M}E_k^{(+)}}\right)=\mathrm{\Omega }^{4l2}.$$
(7.19)
A complementary set arises from (2.10), taken at $`E=\lambda _k^{()}`$ :

$$\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{E_n^{(+)}\mathrm{\Omega }^{2M}\lambda _k^{()}}{E_n^{(+)}\mathrm{\Omega }^{\mathrm{\hspace{0.33em}2}M}\lambda _k^{()}}\right)=\mathrm{\Omega }^{2l+1\alpha }.$$
(7.20)
Together, (7.19) and (7.20) form a set of Bethe ansatz equations of exactly the same form as (2.12) and (2.13), save for the replacement of $`\alpha `$ and $`l`$ on the right-hand sides of (2.12) and (2.13) by $`\stackrel{~}{\alpha }(\alpha ,l)`$ and $`\stackrel{~}{l}(\alpha ,l)`$, respectively. By comparing the left-hand sides and exploiting the analytic properties derived in appendix B, one obtains by another route the second and fifth equivalences (5.2), (7.9):

$$E_k^{(+)}(\alpha ,l)=E_k^{(+)}(\stackrel{~}{\alpha },\stackrel{~}{l})\text{and}\lambda _k^{()}(\alpha ,l)=E_k^{()}(\stackrel{~}{\alpha },\stackrel{~}{l}).$$
(7.21)
The previous approach, which proceeded via the symmetries of the third-order equation, was perhaps more elegant. The advantage of this alternative method is that the only analytic properties used are those of spectral determinants of the second-order equation, and these are known rigorously from the results in appendix B.

The second remark relates to the fact that the lateral connection problem solved by $`T`$ is closely related to $`𝒫𝒯`$-symmetric quantum mechanics. These problems are not, in any obvious sense, self-adjoint, and the reality properties of their spectra have long been of interest . Since the spectral problems on the right-hand sides of (7.10) or (7.11) are self-adjoint for $`\stackrel{~}{l}>1/2`$ (respectively $`1\stackrel{~}{l}>1/2`$), these two identities give us a simple way to understand the reality of the spectra encoded by $`T`$ in these particular cases. However, in appendix B below we will give a proof of reality which is both more general and more direct, and so we will not pursue this any further.

## 8 A relation with Bender-Dunne polynomials

The dualities that we have been discussing have an interesting relationship with the so-called Bender-Dunne polynomials. These were introduced in  as a way of understanding quasi-exact solvability, but here we will also be interested in their properties at general values of the parameters. Briefly, one searches for a solution to (1.1) of the form

$$\psi (x,E,\alpha ,l)=e^{x^4/4}x^{l+1}\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{1}{4}\right)^n\frac{P_n(E,\alpha ,l)}{n!\mathrm{\Gamma }(n+l+3/2)}x^{2n}.$$
(8.1)
For this to solve the differential equation (1.1), the coefficients $`P_n`$ must satisfy the following recursion relation:

$$P_n(E)=EP_{n1}(E)+16(n1)(nJ1)(n+l1/2)P_{n2}(E),(n1)$$
(8.2)
where, as before, $`J=J(\alpha ,l)=(\alpha +2l+1)/4`$. The value of $`P_0(E)`$, which determines the normalisation of $`\psi (x)`$, is conventionally taken to be $`1`$; from (8.2), $`P_1=E`$, and $`P_n`$ is a polynomial of degree $`n`$ in $`E`$, known as a Bender-Dunne polynomial. So long as $`ln3/2`$ for any $`n^+`$, (8.1) will yield an everywhere-convergent series solution to (1.1). Furthermore, this solution automatically satisfies the boundary condition $`\psi x^{l+1}`$ at $`x=0`$ ; but at general values of $`E`$, it will grow exponentially as $`x\mathrm{}`$. We now ask whether there are transformations of the parameters $`\alpha `$ and $`l`$ which leave the Bender-Dunne polynomials invariant. It is easily seen that if $`J`$ and $`l`$ are replaced by $`\stackrel{~}{J}=l1/2`$ and $`\stackrel{~}{l}=J1/2`$, then the recursion relation is unchanged. Translated back to the parameters $`\alpha `$ and $`l`$, this implies that

$$P_n(E,\alpha ,l)=P_n(E,\stackrel{~}{\alpha },\stackrel{~}{l})$$
(8.3)
where

$$\stackrel{~}{\alpha }=2J+4l+2=3/2\alpha /2+3l,\stackrel{~}{l}=J1/2=\alpha /4+l/21/4.$$
(8.4)
This matches the ‘second spectral equivalence’ found earlier.

We will return to this case later, but first we discuss the special points where the model is quasi-exactly solvable, for which a similar game can be played. If $`\alpha `$ and $`l`$ are such that $`J(\alpha ,l)`$ is a positive integer, the second term on the RHS of (8.2) vanishes at $`n=J+1`$, and all the subsequent $`P_n`$ therefore factorise:

$$P_{n+J}(E,\alpha _J,l)=P_J(E,\alpha _J,l)Q_n(E,\alpha _J,l),(n>0,J=(\alpha _J+2l+1)/4).$$
(8.5)
Hence, if $`P_J(E)`$ vanishes then so do all $`P_{nJ}(E)`$ and the series (8.1) terminates, automatically yielding a normalisable solution to (1.1). The $`J`$ zeroes of $`P_J(E)`$ are the $`J`$ exactly-solvable levels for the model and, as observed by Bender and Dunne, this provides a simple way to understand the quasi-exact solvability of the model. Now we would like to go further and discuss the remaining levels. The polynomials $`Q_n`$ satisfy the recursion

$$Q_n(E)=EQ_{n1}(E)+16(n+J1)(n1)(n+J+l1/2)Q_{n2}(E),(n1)$$
(8.6)
with initial conditions $`Q_0=1`$, $`Q_1=E`$. This matches the recursion relation for $`P_n(E)`$, so long as $`J`$ and $`l`$ in (8.2) are replaced by $`\widehat{J}=J`$ and $`\widehat{l}=J+l`$. Hence, if

$$\widehat{\alpha }_J=2J2l1=\alpha _J/23l3/2,\widehat{l}_J=J+l=\alpha _J/4+l/21/4,$$
(8.7)
then

$$Q_n(E,\alpha _J,l)=P_n(E,\widehat{\alpha }_J,\widehat{l}_J).$$
(8.8)
This corresponds to the ‘third spectral equivalence’, and it has an interesting consequence for the series expansion (8.1), which we rewrite as

$`\psi (x,E,\alpha _J,l)`$ $`=`$ $`e^{x^4/4}x^{l+1}\left[\mathrm{}+{\displaystyle \underset{n=J}{\overset{\mathrm{}}{}}}\left(\frac{1}{4}\right)^n{\displaystyle \frac{P_n(E,\alpha _J,l)}{n!\mathrm{\Gamma }(n+l+3/2)}}x^{2n}\right]`$
$`=`$ $`e^{x^4/4}x^{l+1}\left[\mathrm{}+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(\frac{1}{4}\right)^{n+J}{\displaystyle \frac{P_J(E,\alpha _J,l)Q_n(E,\alpha _J,l)}{(n+J)!\mathrm{\Gamma }(n+J+l+3/2)}}x^{2(n+J)}\right]`$ (8.9)
the dots standing for lower-order terms. This can be compared with the expansion of the wavefunction $`\psi (x,E,\widehat{\alpha }_J,\widehat{l}_J)`$. Using $`\widehat{l}_J=J+l`$ and the equality (8.8), this is

$$\psi (x,E,\widehat{\alpha }_J,\widehat{l}_J)=e^{x^4/4}x^{l+J+1}\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{1}{4}\right)^n\frac{Q_n(E,\alpha _J,l)}{n!\mathrm{\Gamma }(n+J+l+3/2)}x^{2n}.$$
(8.10)
It is now easy to see that $`\psi (x,E,\alpha _J,l)`$ is mapped onto a function proportional to $`\psi (x,E,\widehat{\alpha }_J,\widehat{l}_J)`$ by the differential operator

$$𝒬_J^{}(l)=e^{x^4/4}x^{l+J+1}\left(\frac{1}{x}\frac{d}{dx}\right)^Je^{x^4/4}x^{l1}=x^J\left[\frac{1}{x}\frac{d}{dx}+x^2\frac{l+1}{x^2}\right]^J.$$
(8.11)
This is enough to see that the following intertwining relation between differential operators must hold:

$$𝒬_J^{}(l)(\alpha _J(l),l)=(\widehat{\alpha }_J(l),\widehat{l}_J(l))𝒬_J^{}(l).$$
(8.12)
(ConsiderWe would like to thank Peter Bowcock for a discussion of this point. the difference between the LHS and RHS. This is a linear $`(J+2)^{\mathrm{th}}`$-order differential operator, independent of $`E`$, and it is easily seen that it annihilates the functions $`\psi (x,E,\alpha _J(l),l)`$. These functions are linearly independent for different values of $`E`$, while a $`(J+2)^{\mathrm{th}}`$-order operator can annihilate at most $`(J+2)`$ independent functions, unless it is identically zero. This establishes the equality.) We also used Maple to verify (8.12) directly. Finally, $`𝒬_J^{}(l)`$ respects the boundary conditions: if $`\psi (x)`$ decays as $`x\mathrm{}`$ then so does $`𝒬_J^{}(l)\psi (x)`$ , and if $`\psi (x)`$, given as a series by (8.1), has leading behaviour $`x^{l+1}`$ at the origin, then $`𝒬_J^{}(l)\psi (x)`$ has leading behaviour $`x^{l+J+1}`$. Thus $`𝒬_J^{}(l)`$ maps eigenfunctions of the problem $`(\alpha _J,l)`$ to those of $`(\widehat{\alpha }_J(l),\widehat{l}_J(l))`$, or to zero. The eigenfunctions mapped to zero are those for which $`P_J(E,\alpha _J,l)`$ vanishes (the lower-order terms ‘…’ in (8.1) are clearly annihilated by $`𝒬_J^{}(l)`$), and these are precisely the exactly-solvable levels. This provides an alternative derivation of the duality found with the aid of the Bethe ansatz equations in §6, and shows that in $`𝒬_J^{}`$ we have found the generalisation of the supersymmetry operator $`𝒬^+𝒬_1^{}`$ to the QES problems (6.16) with $`J>1`$.

So far we have discussed the action of $`𝒬_J^{}(l)`$ on solutions to the spectral problem with $`x^{l+1}`$ boundary conditions. But given the intertwining relation (8.12) it is natural to look for an action on solutions satisfying the other, $`x^l`$, boundary condition at the origin. As a relation between differential operators, the fact that $`(a,b)=(a,1b)`$ means, trivially, that

$$𝒬_J^{}(l)(\alpha _J(l),1l)=(\widehat{\alpha }_J(l),1\widehat{l}_J(l))𝒬_J^{}(l).$$
(8.13)
It can also be checked that, in general, the relevant boundary conditions are respected, so that (8.13) holds as an intertwining relation between eigenvalue problems. Substituting $`1l`$ for $`l`$ throughout, the relation is

$$𝒬_J^{}(1l)(\alpha _J(1l),l)=(\widehat{\alpha }_J(1l),1\widehat{l}_J(1l))𝒬_J^{}(1l).$$
(8.14)
Thus $`𝒬_J^{}(1l)`$ and its adjoint intertwine between the spectral problems

$$(4J+2l+1,l)\text{ and }(2J+2l+1,lJ),$$
(8.15)
and, in general, no eigenfunctions are annihilated. Furthermore, the mapping (8.15) is exactly the ‘second spectral equivalence’ of §5 above, specialised to cases where the initial pair of parameters $`(\alpha ,l)`$ satisfies $`\alpha =4J+2l+1`$. This hints at an alternative way to obtain (8.14): just as was done at the QES points using (8.8), one can compare the series expansions using (8.3). At a formal level, for any value of $`\alpha `$ and $`l`$ the series for $`\psi (x,E,\alpha ,l)`$ is mapped onto that for $`\psi (x,E,\stackrel{~}{\alpha },\stackrel{~}{l})`$ by

$$𝒫_\mu (l)=e^{x^4/4}x^{\mu l}\left(\frac{1}{x}\frac{d}{dx}\right)^\mu e^{x^4/4}x^l,$$
(8.16)
where $`\mu =\alpha /4+l/2+1/4`$, and the action of a fractional power of the derivative is, again formally, defined by

$$\left(\frac{1}{x}\frac{d}{dx}\right)^\mu x^{2n}=2_{}^\mu \frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(n+1\mu )}x^{2(n\mu )}.$$
(8.17)
In cases when $`\mu `$ is a positive integer, $`𝒫_\mu (l)`$ becomes an ordinary differential operator, equal to $`𝒬_\mu (1l)`$, and (8.14) is recovered. It is interesting to speculate about the existence of some kind of spontaneously-broken ‘fractional’ generalised supersymmetry lying behind the tilde-duality at arbitrary values of the parameters, but we leave this for future work.

In concluding this section, we would like to mention the recent article , which we noticed as we were finishing the writing of this paper. By a completely different route, involving a study of a concept called ‘$`𝒩`$-fold supersymmetry’ §§§generalised supersymmetries of this sort are also called ‘higher-derivative’, or ‘non-linear’ . the authors of this work have also introduced higher-order analogues of the supersymmetry generators. Although the connection with quasi-exact solvability is not mentioned, one can check that the ‘cubic’ case of the type A $`𝒩`$-fold supersymmetry of  reproduces the result (8.12) above, albeit with a slightly different presentation of the operators. (In fact, using the most general form of their operators, one can also obtain an intertwining relation for the more general QES sextic potentials involving an additional $`x^4`$ term.) We should also mention that a connection between certain other forms of non-linear supersymmetry and quasi-exact solvability has recently been pointed out in . However the forms of the supersymmetry generators explicitly treated in that paper do not cover the case of the sextic potential discussed above.

## 9 Conclusions

The main purpose of this paper has been to illustrate how spectral properties and symmetries of interesting differential operators can be handled using tools originally developed in the context of integrable models. These lead to some novel spectral equivalences, and, as will be shown in appendix B below, they also allow for an elementary proof of a reality property which has been surprisingly elusive when studied by more conventional methods. Some of the equivalences we have subsequently been able to re-derive by other means, and in this respect, the rôle of higher-order generalisations of the supersymmetry operators is particularly intriguing, especially in the light of their independent appearance in . It would be very interesting to find out whether the connection between such operators and quasi-exact solvability that we have observed is more general. In this paper, we obtained the operators $`𝒬_J^{}`$ through a direct examination of power series solutions; how they fit into the more algebraic schemes for understanding quasi-exact solvability, as developed in, for example, , is another question that deserves further study.

In many ways, the first and fourth spectral equivalences, between second- and third-order equations, are the most unexpected of our results. They can be traced back to the collapse of the $`SU(3)`$ Bethe ansatz equations at $`N=1`$. A similar phenomenon occurs in $`SU(n)`$-related BA systems for $`n>3`$, which are related to higher-order differential equations via the ODE/IM correspondence . However, in these cases the resulting ‘reduced’ systems are not so readily identifiable, and so at this stage we lack an interpretation of the phenomenon in terms of the properties of differential equations. On the other hand, via the first equivalence we do at least see that quasi-exact solvability is not restricted to second-order spectral problems.

One can also ask whether quasi-exact solvability might have a rôle on the ‘integrable models’ side of the ODE/IM correspondence. Bethe roots correspond to zeroes of $`D(E)`$, and at the QES points the locations of a finite subset of these can be found exactly. From the identity (7.16), the locations of the remaining roots coincide with the zeroes of $`T(E)`$, while, as follows from the form of the Bender-Dunne wavefunctions, the QES roots themselves correspond to coincidences in the locations of zeroes of $`D^{(\pm )}(E)`$ and $`T^{()}(E)`$. However, we do not know of any special significance of these facts.

Finally, it is interesting to draw the full set of spectral problems that can be reached from a QES starting-point using the dualities that we have been discussing:

$`(4J2l1,l)`$ $``$ $`(2J2l1,J+l)`$
$``$ $``$
$`(2J+4l+2,J\frac{1}{2})`$ $``$ $`(2J+4l+2,J\frac{1}{2})`$
Vertical arrows correspond to the second spectral equivalence $`𝕋`$, while the upper horizontal arrow is the level-eliminating third equivalence, $``$. The two problems on the bottom row are related by the transformation $`𝕃`$: $`l1l`$. They correspond to the same Schrödinger equation, and differ only in the boundary condition imposed at the origin. It follows from the diagram that the ‘regular’ eigenvalue problem for this equation, that with the $`x^{J1/2}`$ boundary condition, has exactly the same spectrum as the irregular problem, with the exception of the first $`J`$ levels. This can be understood by noticing that, at $`l=J1/2`$, there is a ‘resonance’ between the regular and irregular solutions of the Schrödinger equation (see, for example, ). The Bender-Dunne series expansions for the first $`J`$ levels of the irregular problem truncate before the resonance is reached, whilst the wavefunctions for the remaining levels are completely dominated by the effect of the resonance, and hence match those of the regular problem. In fact, such a square of spectral problems can be drawn starting from any values of the parameters $`\alpha `$ and $`l`$, on account of the identity $`=𝕋𝕃𝕋`$. But it is only at the QES points that the horizontal directions correspond to (partial) spectral equivalences, since the resonance between regular and irregular solutions just described only occurs when $`J`$ is an integer. Thus we have some novel points at which the sextic potential can be considered to be quasi-exactly solvable, a ‘dual’ interpretation of quasi-exact solubility for this model in terms of the resonance of irregular solutions of the Schrödinger equation, and an alternative interpretation of the level-elimination at work in the ‘supersymmetric’ third spectral equivalence.

Acknowledgements —  We would like to thank Sharry Borgan, Peter Bowcock, Francesco Cannata, Davide Fioravanti, Brent Everitt, Giuseppe Mussardo, John Parker, Junji Suzuki, Alexander Turbiner, André Voros and Jean Zinn-Justin for useful conversations and help. RT also thanks Ferdinando Gliozzi and the Department of Theoretical Physics of Torino University for kind hospitality during the final stages of this project, and we all thank Hideaki Aoyama, Carl Bender, Andrei Mezincescu and Miloslav Znojil for helpful correspondence. The work was supported in part by a TMR grant of the European Commission, reference ERBFMRXCT960012. RT was also supported by an EPSRC VF grant, number GR/N27330.

## Appendix A Solving radial Schrödinger equations using Maple

In this appendix we discuss the numerical treatment of the radial Schrödinger equation. The standard method is to integrate the ODE at varying values of energy, imposing boundary conditions either at the origin or infinity, and searching for the values of energy at which the other boundary condition is satisfied (see for example ). This approach runs into problems in the region $`\mathrm{}el<1/2`$, where the eigenvalue problem is better defined via analytic continuation. The series solution alternative that we describe here avoids this difficulty, can be implemented in only a few lines of Maple and seems to be a rather efficient method for finding the lower-lying levels, at least for a polynomial potential. We chose not to use a series of the sort described in §8 above, as the factor of $`\mathrm{exp}(x^4/4)`$ means that any finite truncation of the series always decays at $`x\mathrm{}`$; it is only in the infinite sum that the exponential growth of a solution at generic $`E`$ is recovered. This makes it hard to detect eigenvalues reliably. Instead, we generated a pure power series directly in Maple, using an algorithm based on the method of Cheng . The power series y produced by the program depends both on x and on E, and by construction it satisfies the boundary condition at x=0. At an eigenvalue, the solution must decay at large x, and by choosing a suitably-large value x0 and searching for values of E at which y(x0,E)=0, the eigenvalues can be located with high accuracy. The value of x0 must be large enough that the asymptotic behaviour of the true solution has set in, but small enough that the approximated power series can be relied on. (The level of the approximation is controlled by the variable iterations in the program below.) This can be checked by examining plots of the candidate wavefunctions. We first give the code that we used, the particular example producing figure A below. The values of $`\alpha `$ and $`l`$ are specified in the second line.

```
> Digits:=20: iterations:=40:
> alpha:=-4: l:=0:
> V:=x^6+alpha*x^2;
> L:=poly->sum(coeff(poly,x,n)*x^(n+2)/(n+2)/(n+2*l+3),n=0..degree(poly,x)):
> P:=1:for i from 1 to iterations do P:= simplify(1+L((V-E)*P)) end do:
> y:=x^(l+1)*P:
> spectrum:=fsolve(eval(y,x=3.2)=0);

```

A specific level, determined by the integer plotlevel, is then plotted as follows:

```
> with(plots): plotlevel:=0:Ee:=spectrum[plotlevel+1];
> xmax:=2.5:ymin:=-35:ymax:=100:
> display([
  plot(eval(80*y,E=Ee),x=0..xmax,ymin..ymax,color=blue,linestyle=2,thickness=2),
  plot((V+l*(l+1)/x^2),x=0..xmax,ymin..ymax,color=red,thickness=2),
  seq(plot([[0.05,spectrum[lev]],[0.12,spectrum[lev]]],
             color=black,linestyle=1,thickness=3), lev=1..7),
  plot([[0.01,spectrum[plotlevel+1]],[0.17,spectrum[plotlevel+1]]],
        color=black,linestyle=1,thickness=1)
  ]);

```

The levels are contained in the list spectrum, and in table 1 below we compare semiclassical, nonlinear integral equation (NLIE) and Maple results for two sides of the $`SU(3)`$-inspired duality of §5. To compile the table, we increased Digits to $`40`$ and iterations to $`50`$; nevertheless, each column of data still took less than 6 minutes of CPU time on a 650 MHz Pentium III machine, running under Linux.

From the table, it is clear that the NLIE (2.15) is able to find the energy levels with high accuracy, provided $`|\alpha |<M+1|2l+1|`$, and for such values of the parameters it seems to be the most reliable method to find the full set of energy levels. The power series approach also works extremely well for the low-lying energy levels, but loses accuracy for the higher levels, at least if the value of iterations is kept reasonably small. Nevertheless, combining power series and semiclassical methods allows us to obtain good results, over the full spectrum, for any values of $`\alpha `$ and $`l`$.

The figures which follow illustrate some typical wavefunctions. In each case the low-lying energy levels are shown as bars near the left-hand side of the plot, with the level under scrutiny having double length. The dashed line is the (un-normalised) wavefunction associated with this energy level, and the solid line the potential itself.

The first set of four figures illustrates the second spectral equivalence, connecting the potentials with $`\alpha ,l=4,0`$ and $`7/2,5/4`$ respectively. Note that the dual problem is an example where the eigenfunction is not square-integrable; nevertheless, the numerical solution converges and reproduces exactly the required spectrum. Figures A and A show the ground state, and figures A and A the first excited state.

Figure A:
=α-4, =l0, =E01.0057683Figure Aα=7/2l=5/4E0=1.0057683
Figure A:
=α-4, =l0, =E01.0057683Figure Aα=7/2l=5/4E0=1.0057683\begin{array}[]{ll}\includegraphics[width=199.4681pt,height=199.4681pt]{pa1.ps}{}{}{}{}&\includegraphics[width=199.4681pt,height=199.4681pt]{pa2.ps}\\
\parbox{195.12767pt}{
\raggedright{\small Figure \ref{pa1}:
$\alpha=-4$, $l=0$, $E_{0}=1.0057683$
}\@add@raggedright}&\parbox{216.81pt}{
\raggedright{\small Figure \ref{pa2}:
$\alpha=7/2$, $l=-5/4$, $E_{0}=1.0057683$
}\@add@raggedright}\end{array}
Figure A:
=α-4, =l0, =E110.572585Figure Aα=7/2l=5/4E1=10.572585
Figure A:
=α-4, =l0, =E110.572585Figure Aα=7/2l=5/4E1=10.572585\begin{array}[]{ll}\includegraphics[width=199.4681pt,height=199.4681pt]{ppa1.ps}{}{}{}{}&\includegraphics[width=199.4681pt,height=199.4681pt]{ppa2.ps}\\
\parbox{195.12767pt}{
\raggedright{\small Figure \ref{ppa1}:
$\alpha=-4$, $l=0$, $E_{1}=10.572585$
}\@add@raggedright}&\parbox{216.81pt}{
\raggedright{\small Figure \ref{ppa2}:
$\alpha=7/2$, $l=-5/4$, $E_{1}=10.572585$
}\@add@raggedright}\end{array}
Figure A depicts the quasi-exactly solvable case for the third energy level, at $`J=2`$ and $`l=0`$. This corresponds to the ground state of the SUSY ($``$) dual potential, shown in figure A. The fourth level of the QES problem has the same energy as the first excited state of the SUSY dual, and these two are shown in figures A and A. Figures A and A then show the first two states of the potential related to the SUSY dual by the second spectral equivalence, the tilde-duality $`𝕋`$.

Figure Aα=9l=0E2=16.919850
Figure A:
=α3, =l2, =E016.919850Figure Aα=9l=0E2=16.919850
Figure A:
=α3, =l2, =E016.919850\begin{array}[]{ll}\includegraphics[width=199.4681pt,height=199.4681pt]{p1.ps}{}{}{}{}&\includegraphics[width=199.4681pt,height=199.4681pt]{p2.ps}\\[3.0pt]
\parbox{195.12767pt}{
\raggedright{\small Figure \ref{p1}:
$\alpha=-9$, $l=0$, $E_{2}=16.919850$
}\@add@raggedright}&\parbox{216.81pt}{
\raggedright{\small Figure \ref{p2}:
$\alpha=3$, $l=2$, $E_{0}=16.919850$
}\@add@raggedright}\end{array}
$$\begin{array}{cc}\text{}\hfill & \text{}\hfill \\ \text{Figure }\text{A}\text{}\alpha =9\text{}l=0\text{}E_3=32.240265\hfill & \text{Figure }\text{A}\text{}\alpha =3\text{}l=2\text{}E_1=32.240265\hfill \end{array}$$

Figure A:
=α6, =l/32, =E016.919850Figure Aα=6l=3/2E1=32.240265
Figure A:
=α6, =l/32, =E016.919850Figure Aα=6l=3/2E1=32.240265\begin{array}[]{ll}\includegraphics[width=199.4681pt,height=199.4681pt]{p3.ps}{}{}{}{}&\includegraphics[width=199.4681pt,height=199.4681pt]{pp3.ps}\\[3.0pt]
\parbox{195.12767pt}{
\raggedright{\small Figure \ref{p3}:
$\alpha=6$, $l=3/2$, $E_{0}=16.919850$
}\@add@raggedright}&\parbox{195.12767pt}{
\raggedright{\small Figure \ref{pp3}:
$\alpha=6$, $l=3/2$, $E_{1}=32.240265$
}\@add@raggedright}\end{array}
Finally, in figures A and A we illustrate how ‘extra’ energy levels can appear with the irregular boundary condition in resonance situations, a phenomenon that was mentioned at the end of the conclusions above. Again, we take $`J=2`$. To avoid numerical difficulties, in figure A we shifted the angular momentum slightly away from the exactly-resonant value.

$$\begin{array}{cc}\text{}\hfill & \text{}\hfill \\ \text{Figure }\text{A}\text{}\alpha =9\text{}l=0\text{}E_0=4.898979\hfill & \text{Figure }\text{A}\text{}\alpha =6\text{}l=5/2+10^6\text{}E_0=4.898974\hfill \end{array}$$

## Appendix B An elementary proof of the Bessis, Zinn-Justin, Bender and Boettcher conjecture

A conjecture of Bessis and Zinn-Justin , generalised by Bender and Boettcher , states that the eigenvalues $`\lambda _k`$ of the $`𝒫𝒯`$-symmetric Schrödinger equation

$$\left[\frac{d^2}{dx^2}(ix)^{2M}\right]\psi _k(x)=\lambda _k\psi _k(x),\psi _k(x)L^2(𝒞)$$
(B.1)
are real and positive for $`M1`$. The contour $`𝒞`$ on which the wavefunction is defined can be taken to be the real axis for $`M<2`$ ; beyond this point, the contour should be deformed down into the complex plane so as to remain in the same pair of Stokes sectors . (For an informal review in the context of the ODE/IM correspondence, see also .) This conjecture has provoked a fair amount of work in recent years, a sample being refs. . In this appendix we consider a slightly more general class of $`𝒫𝒯`$-symmetric spectral problems, namely

$$\left[\frac{d^2}{dx^2}(ix)^{2M}\alpha (ix)^{M1}+\frac{l(l+1)}{x^2}\right]\psi _k(x)=\lambda _k\psi _k(x),\psi _k(x)L^2(𝒞)$$
(B.2)
with $`M`$, $`\alpha `$ and $`l`$ real. Again, for $`M<2`$ the contour $`𝒞`$ can be taken to be the real axis, though if $`l(l+1)0`$ it should be distorted so as to pass below the origin. We shall prove reality of the spectrum for $`M>1`$, $`\alpha <M+1+|2l+1|`$, and positivity for $`M>1`$, $`\alpha <M+1|2l+1|`$. The spectrum might be real for a greater range of $`\alpha `$, but strict positivity certainly fails on the lines $`\alpha =M+1|2l+1|`$. Even with the restrictions on $`\alpha `$, our result includes the previously-considered cases: for $`\alpha =l(l+1)=0`$ and $`M=3/2`$, a version of the original Bessis – Zinn-Justin conjecture is recovered; allowing $`M`$ to vary then gives the generalisation discussed by Bender and Boettcher, while the conjecture for $`\alpha =0`$ and $`l`$ small was proposed in . (Strictly speaking the original BZ-J conjecture concerned the potential $`x^2+igx^3`$ with $`g`$ real; our discussion applies to the strong-coupling limit of this problem.)

Setting $`\mathrm{\Phi }(x)=\psi (x/i)`$, (B.2) becomes

$$\left[\frac{d^2}{dx^2}+x^{2M}+\alpha x^{M1}+\frac{l(l+1)}{x^2}\right]\mathrm{\Phi }_k(x)=\lambda _k\mathrm{\Phi }_k(x),\mathrm{\Phi }_k(x)L^2(i𝒞),$$
(B.3)
and has the same form as (2.1) with $`E=\lambda _k`$, though with different boundary conditions: to qualify as an eigenfunction, $`\mathrm{\Phi }`$ must decay as $`|x|\mathrm{}`$ along the contour $`i𝒞`$. However, it is an easy generalisation of the discussion in §7 of  that the function $`T(\lambda ,\alpha ,l)`$ defined in (2.5) is the spectral determinant associated to the spectral problem (B.3). This identification allows us to study the generalised BZ-JBB conjecture (B.2) using techniques inspired by the Bethe ansatz.

We start from equation (2.10):

$$T^{(+)}(E)D^{(+)}(E)=\mathrm{\Omega }^{(2l+1+\alpha )/2}D^{()}(\mathrm{\Omega }^{2M}E)+\mathrm{\Omega }^{(2l+1+\alpha )/2}D^{()}(\mathrm{\Omega }^{2M}E),$$
(B.4)
and define the zeroes of $`T^{(+)}(E)=T(E,\alpha ,l)`$ to be the set $`\{\lambda _k\}`$. (Note that for $`\alpha =0`$, (B.4) reduces to the T-Q system obtained in .) Putting $`E=\lambda _k`$ in (B.4) and using, for $`M>1`$, the factorised form for $`D^{()}(E)`$ gives the following constraints on the $`\lambda _k`$’s:

$$\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{E_n^{()}+\mathrm{\Omega }^{2M}\lambda _k}{E_n^{()}+\mathrm{\Omega }^{2M}\lambda _k}\right)=\mathrm{\Omega }^{2l1\alpha },k=0,1,\mathrm{}$$
(B.5)
Since the original eigenproblem (B.2) is invariant under $`l1l`$, we can assume $`l1/2`$ without any loss of generality. Then each $`E_n^{()}`$ is an eigenvalue of an Hermitian operator $`(M,\alpha ,l)`$, and hence is real. Furthermore a Langer transformation  (see also ) shows that the $`E_n^{()}`$ solve a generalised eigenproblem with an everywhere-positive ‘potential’, and so are all positive, for $`\alpha <1+2l`$. This can be sharpened by considering the value of $`D^{()}(E)|_{E=0}^{}`$ . From (B.15) below, this first vanishes when $`\alpha =M+2l+2`$. Until this point is reached, no eigenvalue $`E_n^{()}`$ can have passed the origin, and all must be positive. (It might be worried that negative eigenvalues could appear from $`E=\mathrm{}`$ , but this possibility can be ruled out by a consideration of the Langer-transformed version of the equation.)

Taking the modulus2 of (B.5), using the reality of the $`E_k^{()}`$, and writing the eigenvalues of (B.2) as $`\lambda _k=|\lambda _k|\mathrm{exp}(i\delta _k)`$, we have

$$\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{(E_n^{()})^2+|\lambda _k|^2+2E_n^{()}|\lambda _k|\mathrm{cos}(\frac{2\pi }{M+1}+\delta _k)}{(E_n^{()})^2+|\lambda _k|^2+2E_n^{()}|\lambda _k|\mathrm{cos}(\frac{2\pi }{M+1}\delta _k)}\right)=1.$$
(B.6)
For $`\alpha <M+2l+2`$ , all the $`E_n^{()}`$ are positive, and each single term in the product on the LHS of (B.6) is either greater than, smaller than, or equal to one depending only on the relative values of the cosine terms in the numerator and denominator. These are independent of the index $`n`$. Therefore the only possibility to match the RHS is for each term in the product to be individually equal to one, which for $`\lambda _k^{}0`$ requires

$$\mathrm{cos}(\frac{2\pi }{M+1}+\delta _k)=\mathrm{cos}(\frac{2\pi }{M+1}\delta _k),\text{or}\mathrm{sin}(\frac{2\pi }{M+1})\mathrm{sin}(\delta _k)=0.$$
(B.7)
Since $`M>1`$, this latter condition implies

$$\delta _k=n\pi ,n$$
(B.8)
and this establishes the reality of the eigenvalues of (B.2) for $`M>1`$ and $`\alpha <M+2l+2`$ or, relaxing the condition on $`l`$, $`\alpha <M+1+|2l+1|`$.

One might ask what goes wrong for $`M<1`$, since from  (and, for the case $`l0`$, ) it is known that most of the $`\lambda _k`$ become complex as $`M`$ falls below $`1`$, at least for $`\alpha =0`$. The answer is that if $`M<1`$, the order of $`D^{()}(E)`$ is greater than $`1`$, the factorised form of $`D^{()}(E)`$ provided by Hadamard’s theorem no longer has such a simple form, and the proof breaks down.

The borderline case $`M=1`$ is the simple harmonic oscillator, exactly solvable for all $`l`$ and $`\alpha `$ . Starting from the discussion in §3 of , it is easily seen that

$$T(E,\alpha ,l)|_{M=1}^{}=\frac{2\pi }{\mathrm{\Gamma }\left(\frac{1}{2}+\frac{2l+1+E\alpha }{4}\right)\mathrm{\Gamma }\left(\frac{1}{2}\frac{2l+1E+\alpha }{4}\right)}$$
(B.9)
and so the eigenvalues of (B.2) are at $`\lambda =4n+2\alpha \pm (2l+1)`$, $`n=0,1,\mathrm{}`$ . All are real for all real values of $`\alpha `$ and $`l`$, and all are positive for $`\alpha <2|2l+1|`$.

To discuss positivity at general values of $`M>1`$, we can continue in $`M`$, $`\alpha `$ and $`l`$ away from a point in this latter region, $`\{M=1,\alpha <2|2l+1|\}`$. So long as $`\alpha `$ remains less than $`M+1+|2l+1|`$, all eigenvalues will be confined to the real axis during this process, and the first passage of an eigenvalue from positive to negative values will be signalled by the presence of a zero in $`T(\lambda ,\alpha ,l)`$ at $`\lambda =0`$. Fortunately, $`T(\lambda ,\alpha ,l)|_{\lambda =0}^{}`$ can be calculated exactly, extending an argument given for $`\alpha =0`$ in . First, one notices that the function

$$\phi (x)=\left(\frac{M+1}{2}\right)^{\frac{M+\alpha }{2M+2}}x^{\frac{M1}{2M+2}}Y(\left(\frac{M+1}{2}\right)^{\frac{1}{M+1}}x^{\frac{2}{M+1}},E,\alpha ,l)$$
(B.10)
solves the Schrödinger equation

$$\left[\frac{d^2}{dx^2}+x^2\sigma x^{\frac{22M}{M+1}}+\frac{\gamma (\gamma +1)}{x^2}\right]\phi (x)=\mathrm{\Lambda }\phi (x),$$
(B.11)
where

$$\sigma =\left(\frac{2}{M+1}\right)^{\frac{2M}{M+1}}E,\gamma =\frac{2l+1}{M+1}\frac{1}{2},\mathrm{\Lambda }=\frac{2\alpha }{M+1}.$$
(B.12)
(This transformation, which can be found via a pair of Langer transformations, leads to equation (1.2) in the case $`M=3`$.) Further, $`\phi (x)`$ has the large-$`x`$ asymptotic

$$\phi (x)\frac{1}{\sqrt{2i}}x^{\frac{1}{2}\frac{\alpha }{M+1}}\mathrm{exp}\left(\frac{1}{2}x^2\right).$$
(B.13)
At $`E=0`$, $`\sigma =0`$ and (B.11) is the simple harmonic oscillator, which can be solved exactly in terms of the confluent hypergeometric function $`U(a,b,z)`$. Matching asymptotics at large $`x`$,

$$\phi (x)|_{E=0}^{}=\frac{1}{\sqrt{2i}}x^{\gamma +1}e^{x^2/2}U(\frac{1}{2}(\gamma +\frac{3}{2})\frac{1}{4}\mathrm{\Lambda },\gamma +\frac{3}{2},x^2).$$
(B.14)
Reversing the variable changes, extracting the leading behaviour as $`x0`$ and comparing with (2.9), we find

$$D(E,\alpha ,l)|_{E=0}^{}=D^{(+)}(E)|_{E=0}^{}=\frac{1}{\sqrt{2i}}\left(\frac{M+1}{2}\right)^{\frac{2l+1\alpha }{2M+2}\frac{1}{2}}\frac{\mathrm{\Gamma }\left(\frac{2l+1}{M+1}\right)}{\mathrm{\Gamma }\left(\frac{2l+1+\alpha }{2M+2}+\frac{1}{2}\right)}.$$
(B.15)
Now $`T(E,\alpha ,l)|_{E=0}^{}`$ follows from (B.4), remembering that $`D^{()}(E)=D(E,\alpha ,l)`$ :

$$T(E,\alpha ,l)|_{E=0}^{}=T^{(+)}(E)|_{E=0}^{}=\left(\frac{M+1}{2}\right)^{\frac{\alpha }{M+1}}\frac{2\pi }{\mathrm{\Gamma }\left(\frac{1}{2}+\frac{2l+1\alpha }{2M+2}\right)\mathrm{\Gamma }\left(\frac{1}{2}\frac{2l+1+\alpha }{2M+2}\right)}.$$
(B.16)
The first zero arrives at $`E=\lambda =0`$ when $`\alpha =M+1|2l+1|`$, and so for all $`\alpha <M+1|2l+1|`$ , the spectrum is entirely positive, as claimed.

In finishing, we return to the reality of the spectrum encoded by $`T(\lambda )`$. We have proved that, if $`M>1`$, the eigenvalues $`\lambda _k`$ are real for all real $`\alpha <M+1+|2l+1|`$. One might conjecture that this reality should hold for all real $`\alpha `$ and $`l`$. However, this is definitely not the case: for $`M=3`$, at the QES points and with $`l`$ sufficiently negative, an examination of the Bender-Dunne polynomials shows that the exactly-calculable part of the spectrum of $`D(E,\alpha \alpha _J)`$ has at least one pair of complex-conjugate eigenvalues. The identity $`T(E,\alpha \alpha _J)=\gamma (\alpha \alpha _J)D(E,\alpha \alpha _J)`$, which follows from the results obtained in §7 above, shows that $`T(E,\alpha \alpha _J)`$ must share these complex zeroes. If $`(\alpha _J,l)`$ are real then so are $`(\widehat{\alpha }_J,\widehat{l}_J)`$, and so such examples demonstrate that $`T`$ can have complex zeroes even while $`M>1`$ and $`\alpha `$ and $`l`$ are real. It would be worthwhile to map out the full extent of the region where the spectrum is entirely real, but we have not yet done this, beyond a quick check that it appears to extend at least some way beyond the domain $`\alpha <M+1+|2l+1|`$ covered by the proof given in this appendix.

Note added in proof.
(1) We have recently obtained some further results on the region within which the spectrum of the $`𝒫𝒯`$-symmetric problem discussed in appendix B becomes complex. These can be found in .
(2) An alternative treatment of a class of $`𝒫𝒯`$-symmetric quantum mechanical problems similar to those discussed in appendix B can be found in .