# Painlevé equation, anharmonic oscillators and orthogonality

M. Bertola, E. Chavez, T. Grava

Exactly solvable anharmonic oscillator,

degenerate orthogonal polynomials and Painlevé II

M. Bertola†‡⋆<EMAIL_ADDRESS>sissa.it} E. Chavez-Heredia‡♢
<EMAIL_ADDRESS>T. Grava ‡♢<EMAIL_ADDRESS>

  1. ${\dagger}$

Department of Mathematics and Statistics, Concordia University
1455 de Maisonneuve W., Montréal, Québec, Canada H3G 1M8

  2. ${\ddagger}$

SISSA, International School for Advanced Studies, via Bonomea 265, Trieste,
Italy and INFN sezione di Trieste

  3. ${\star}$

Centre de recherches mathématiques, Université de Montréal
C. P. 6128, succ. centre ville, Montréal, Québec, Canada H3C 3J7

  4. ${\diamondsuit}$

School of Mathematics, University of Bristol, Fry Building, Bristol, BS8 1UG,
UK

Abstract.

The paper addresses a conjecture of Shapiro and Tater on the similarity
between two sets of points in the complex plane; on one side is the values of
$t\in\mathbb{C}$ for which the spectrum of the quartic anharmonic oscillator
in the complex plane

$\derivative[2]{y}{x}-\left(x^{4}+tx^{2}+2Jx\right)y=\Lambda y,$

with certain boundary conditions, has repeated eigenvalues. On the other side
is the set of zeroes of the Vorob’ev-Yablonskii polynomials, i.e. the poles of
rational solutions of the second Painlevé equation. Along the way, we indicate
a surprising and deep connection between the anharmonic oscillator problem and
certain degenerate orthogonal polynomials.

###### Contents

  1. 1 Introduction and results
  2. 2 Link to Painlevé
    1. 2.1 Gauge transformation
  3. 3 A study of quasi-polynomials
    1. 3.1 Stokes phenomenon
    2. 3.2 Connection to degenerate orthogonal polynomials
    3. 3.3 Repeated eigenvalues
  4. 4 Exact WKB analysis
    1. 4.1 Schrödinger, Riccati and WKB
    2. 4.2 Stokes graphs and connection formulæ.
  5. 5 Quantization conditions
    1. 5.1 The Jimbo-Miwa case
    2. 5.2 The Shapiro-Tater case
      1. 5.2.1 Repeated eigenvalue condition
    3. 5.3 Comparison of quantization conditions
    4. 5.4 The elliptic region
    5. 5.5 Analysis of the two lattices
  6. A Contribution of the integrals near the turning points
    1. A.1 Proof of Theorem 5.4
  7. B Numerical verifications
    1. B.1 Overlapping patterns

##  1 Introduction and results

The second Painlevé equation is an ODE in the complex domain given by

$\derivative[2]{u}{t}=2u^{3}+tu+\alpha,$ (PII)

with $\alpha\in\mathbb{C}$ and $t\in\mathbb{C}$. It admits rational solutions
when $\alpha=n\in\mathbb{Z}$ and this was recognized by Vorob’ev and
Yablonskii in two separate papers [Vor65, Yab59]. For $\alpha=n$ there are
$n^{2}$ poles of the rational solution, of which $n(n-1)/2$ correspond to
poles with residue $+1$ and the remaining $n(n+1)/2$ correspond to poles with
residue $-1$. Both sets of poles are zeros of certain polynomials defined
recursively (see below), and that are referred to as Vorob’ev-Yablonsky (1)
polynomials.

The study of the asymptotic behaviour of the poles of rational solutions of
Painlevé equations has received significant attention both in the community of
researchers interested in Painlevé theory [CM03, Mas18] and also due to their
occurrence in the description of asymptotic behaviour in the semiclassical
limit of the Sine-Gordon equation [BM12, BM14, BM15, BB15]. The poles form a
very regular “triangular” pattern as it can be seen in Fig. 2.

A seemingly disconnected problem consists in the study of the spectrum of the
following boundary value problem for the anharmonic oscillator

$\displaystyle\derivative[2]{y(x)}{x}-\left(x^{4}+tx^{2}+2Jx\right)y(x)=\Lambda
y(x)$ (1.1) $\displaystyle y(x)\to 0\mbox{ as $x\to\infty$ and
$\arg(x)=\pi,\pm\pi/3$,}$ (1.2)

where $t$ and $J$ are in general complex parameters, and $x$, $y(x)$ are the
complex independent and dependent variable, respectively. In other words (1.1)
is an ODE in the complex domain.

The equation is equivalent to a tri-confluent Heun equation by a gauge
transformation $y(x)=p(x){\rm e}^{x^{3}/3+t/2x}$. Following the definition in
[BB98] a quantum mechanical potential is Quasi Exactly Solvable (QES) if a
finite portion of the energy spectrum and associated eigenfunctions can be
found exactly and in closed form, while the potential is Exactly Solvable (ES)
if all the spectrum and associated eigenfunctions can be found exactly and in
closed form. In particular it was shown in [ST22] that the eigenvalue problem
(1.1) with only two out of the three boundary conditions at infinity is Quasi
Exactly Solvable. Here we show that the quartic anharmonic oscillator in
(1.1), with the boundary conditions (1.2) is Exactly Solvable and admits a
solution if and only if $J=n+1\in\mathbb{N}$ is a positive integer. The
corresponding eigenfunctions are quasi–polynomials (i.e. a polynomial times a
fixed exponential factor).

These polynomials have been studied in the works of Eremenko and Gabrielov
[EG11, EG12] and also Mukhin and Tarasov [MT13]. We actually find several
additional structures relating these polynomials to a special type of
orthogonality (see below). For every given $t\in\mathbb{C}$, the corresponding
eigenvalues of these quasi–polynomial eigenfunctions are the zeros of a
secular equation of degree $n+1$ in the shifted variable
$\lambda=\Lambda-\frac{t^{2}}{4}$. There are particular discrete values of
$t\in\mathbb{C}$ for which the spectrum of the boundary problem (1.1)-(1.2)
becomes degenerate, namely, one of the eigenvalues has higher algebraic
multiplicity.

For a fixed $J\in\mathbb{N}$, these particular values of the parameter
$t\in\mathbb{C}$ form a pattern shown in Fig. 2 that is ostensibly very
similar to the pattern of the poles of the rational solutions of the Painlevé
II equation (i.e. the zeroes of the Vorob’ev–Yablonskii polynomials).

Figure 1: Scaled roots of the Vorob’ev-Yablonsky polynomials $Y_{n}(n^{2/3}s)$
in red, and roots of the discriminant $D_{n}(n^{2/3}s)$ in black, for $n=30$.
This particular scaling was conjectured by Shapiro and Tater in [ST22].

This similarity is the object of a conjecture that B. Shapiro and M. Tater
floated several years ago, but only recently formalized in [ST22] (Conjecture
2 ibidem).

The present paper addresses this Shapiro-Tater (ST) conjecture. Further, along
the way, we investigate the rather surprisingly deep connections of the
boundary value problem (1.1)-(1.2) with “degenerate” orthogonal polynomials:
in fact we show that the polynomials arising from the solutions of (1.1)-(1.2)
satisfy an excess of orthogonality conditions. We also quantify the conjecture
by estimating the discrepancy of the patterns of the two lattices coming from
the VY polynomials on one side and the ST boundary problem on the other,
explaining also why they have the regular pattern that appears numerically.

For the zeros of the VY polynomials, the analysis of these patterns was
extensively already explained in [BM14, BM15, BB15]; the rather surprising
result is the explanation of the connection between the two seemingly distant
problems that we have briefly outlined above.

##### Detailed results.

We now explain the conjecture in more detail. Suppose we are looking for
quasi-polynomials solutions of (1.1), that is, solutions of the form

$y(x)=p(x)e^{\theta(x)},\quad\text{ where
}\theta(x)=\frac{x^{3}}{3}+\frac{tx}{2}$ (1.3)

and $p(x)$ is a polynomial. We may sometimes use the notation $\theta(z;t)$ to
emphasize the dependence of $\theta$ on the parameter t. A substitution into
(1.1) leads to

$\derivative[2]{p}{x}+(2x^{2}+t)\derivative{p}{x}-2(J-1)xp=\lambda\,p,\;\;\text{with
}\lambda=\Lambda-\frac{t^{2}}{4}.$

Upon setting $J=n+1$, one notices that the ODE

$\derivative[2]{p}{x}+(2x^{2}+t)\derivative{p}{x}-2nxp=\lambda p,$ (1.4)

preserves the finite dimensional linear space of polynomials of degree at most
$n$. The left hand side of (1.4) is a linear operator acting on the space of
polynomials of degree at most $n$ that can be written as a $(n+1)\times(n+1)$
matrix $M_{n}(t)$ with respect to the usual basis
$(1,x,x^{2},\dots,x^{n})^{t}$, namely

$M_{n}(t):=\begin{bmatrix}0&-2n&&&&&\\\ t&0&-2(n-1)&&&&\\\ 2&2t&0&&&&\\\
&6&3t&0&\ddots&&\\\ &&12&\ddots&\ddots&-4&\\\ &&&\ddots&&0&-2\\\
&&&&n(n-1)&nt&0\end{bmatrix}\,.$ (1.5)

We observe that the matrix $M_{n}(t)$ is a $4$-diagonal matrix. Eigenvalues of
this matrix correspond to the eigenvalues $\lambda$ in (1.4), and the
eigenvectors correspond to quasi-polynomial eigenfunctions of the differential
operator (1.1) via the identification (1.3). Therefore we can compute the
spectrum of (1.1) from the characteristic polynomial

$C_{n}(t,\lambda):=\det(M_{n}(t)-\lambda I).$ (1.6)

Now consider the discriminant of $C_{n}(t,\lambda)$ with respect to $\lambda$,
i.e.

$D_{n}(t):=\text{Disc}_{\lambda}(C_{n})(t).$ (1.7)

The roots of $D_{n}(t)=0$ are precisely the values of $t$ such that the matrix
$M_{n}(t)$ has repeated eigenvalues (algebraic multiplicity greater than 1).

The Shapiro-Tater conjecture [ST22] relates the zeros of $D_{n}(t)$ to the
poles of rational solutions of PII. The Painlevé II equation has rational
solutions $u(t)$ if and only if $\alpha=n\in\mathbb{Z}$, and we denote them by
$u_{n}(t)$. Furthermore $u_{n}(t)$ is given explicitly in terms of the
Vorob’ev-Yablonsky polynomials $Y_{n}(t)$ in the form:

$u_{n}(t)=\derivative{t}\log\frac{Y_{n-1}(t)}{Y_{n}(t)}$ (1.8)

The VY polynomials are constructed recursively [Vor65, Yab59] as follows:

${Y}_{n+1}(t){Y}_{n-1}(t)=t{Y}_{n}^{2}(t)-4\Big{[}{Y}_{n}^{\prime\prime}(t){Y}_{n}(t)-\big{(}{Y}_{n}^{\prime}(t)\big{)}^{2}\Big{]},\hskip
14.22636ptn\geq 1,\ t\in\mathbb{C}$

with ${Y}_{0}(t)=1,{Y}_{1}(t)=t$. They can also be represented in terms of
Schur function indexed by the “staircase partition” [KO96]

$Y_{n}(t)=\left(-\frac{4}{3}\right)^{n(n+1)/6}\left(\prod_{k=1}^{n}(2k-1)!!\right)S_{(n,n-1,\dots,1)}\left(\left(-\frac{3}{4}\right)^{\frac{1}{3}}t,0,1,0,0,\dots\right).$
(1.9)

The conjecture of Shapiro and Tater [ST22] is the numerical observation that
can be summarized in the conjecture loosely formulated below.

######  Conjecture 1.1.

The roots of the rescaled discriminant
$D_{n}(s)=\text{Disc}_{\lambda}(C_{n})(n^{\frac{2}{3}}s)$ in (1.7) and the
roots of the rescaled Vorob’ev-Yablonsky polynomials $Y_{n}(n^{\frac{2}{3}}s)$
form two coinciding lattices as $n\to\infty$.

The evidence leading Shapiro and Tater to make their conjecture was purely
numerical, as seen in Fig. 2. The numerical picture seems so precise that one
may at first be tempted to compare the polynomials $Y_{n}$ with $D_{n}$;
however a simple inspection shows that their coefficients are not close to
each other, as seen from Table 1.

$n$ |  $D_{n}(t)$
---|---
1 |  $t$
2 |  $\displaystyle t^{3}+\frac{27}{8}$
3 |  $\displaystyle t^{6}+\frac{35}{2}t^{3}-\frac{243}{4}$
4 |  $\displaystyle t^{10}+\frac{215}{4}t^{7}+\frac{89}{8}t^{4}+\frac{4084101}{512}t$
5 |  $\displaystyle t^{15}+\frac{255}{2}t^{12}+\frac{76211}{32}t^{9}+\frac{3730405}{64}t^{6}-\frac{8700637815}{4096}t^{3}-\frac{125005275}{32}$

$n$ |  $Y_{n}(t)$
---|---
1 |  $t$
2 |  $t^{3}+4$
3 |  $t^{6}+20t^{3}-80$
4 |  $t^{10}+60t^{7}+11200t$
5 |  $t^{15}+140t^{12}+2800t^{9}+78400t^{6}-3136000t^{3}-6272000$

Table 1: The first five monic Vorob’ev–Yablonskii polynomials $Y_{n}(t)$ and
discriminant polynomials $D_{n}(t)$ .

######  Remark 1.2.

When $\alpha=n$ the rational solutions $u_{n}(t)$ of (PII) have two types of
poles: those of residue $+1$ and those of residue $-1$. This follows from the
fact that all the zeroes of $Y_{n}(t)$ are simple, $Y_{n}(t),Y_{n+1}(t)$ do
not share any roots [Tan00]and that rational solutions have the shape

$u_{n}(t)=\derivative{t}\log\frac{Y_{n-1}(t)}{Y_{n}(t)}=\frac{Y^{\prime}_{n-1}(t)}{Y_{n-1}(t)}-\frac{Y^{\prime}_{n}(t)}{Y_{n}(t)}.$
(1.10)

Therefore for fixed $n\in\mathbb{N}$ the poles of residue $+1$ correspond to
zeroes of $Y_{n-1}(t)$ and poles of residue $-1$ correspond to zeroes of
$Y_{n}(t)$

The results of the paper can be grouped into two categories: the first
consists in “structural” analysis of the boundary value problem (1.1)-(1.2),
in Section 3. The second category involves the asymptotic study for large $n$
and the description of the ST lattice, Sections 4, 5.1, 5.2. The first set of
results can be summarized in the following points:

  1. 1.

The boundary problem (1.1)-(1.2) has solutions if and only if
$J=n+1\in\mathbb{N}$; in this case $\Lambda$ takes at most $n+1$ values for
each $t\in\mathbb{C}$. See Proposition 3.2. For fixed $J\in\mathbb{N}$ these
values $(t,\Lambda)$ are called the Exactly-Solvable (ES) spectrum.

  2. 2.

For given $J=n+1\in\mathbb{N}$ and $(t,\Lambda)$ in the ES spectrum, the
corresponding solution of (1.1) is of the form $y(x)=p_{n}(x){\rm
e}^{2\theta(x;t)}$ with $p_{n}(x)$ a polynomial of degree $n$ and
$\theta(x;t)=\frac{x^{3}}{3}+\frac{tx}{2}$. This polynomial is degenerate
orthogonal in the sense that

$\displaystyle\left(\varkappa\int_{\infty_{1}}^{\infty_{3}}+\widetilde{\varkappa}\int_{\infty_{3}}^{\infty_{5}}\right)p_{n}(x)x^{j}{\rm
e}^{2\theta(x;t)}{\mathrm{d}}x=0,\ \ j=0,\dots,n,$ (1.11)

where $\varkappa$ and $\tilde{\varkappa}$ are some constants. Here
$\infty_{k}$ denotes a contour that tends to infinity with asymptotic
direction of argument $\frac{i\pi k}{3}$ with $k=1,3,5$. Observe that the
orthogonality involves all powers including the $n$-th power, and this means
that the $(n+1)$-st Hankel determinant of the moments of the above pairing is
zero. This result is contained in Section 3.2 and in particular Theorem 3.4.
These results extend the several properties established for these polynomials
in [EG11, EG12, MT13].

  3. 3.

In fact more is true: given $t\in\mathbb{C}$ suppose that
$\varkappa,\widetilde{\varkappa}$ are such that (1.11) admits a nontrivial
polynomial $p_{n}$ for solution. Then $y(x)=p_{n}(x){\rm
e}^{\frac{x^{3}}{3}+tx}$ is a solution of the boundary problem (1.1)–(1.2).
This is proved in Theorem 3.6.

  4. 4.

The two above results, together, yield a complete characterization of the ES
spectrum of the quartic anharmonic oscillator with boundary conditions (1.2)
in terms of degenerate orthogonality. See Corollary 3.8

  5. 5.

We then investigate the consequences of the requirement that an eigenvalue is
repeated; we prove that this is equivalent to the additional condition that
both integrals in (1.11) for the degenerate orthogonal polynomials vanish
independently, that is equivalent to:

$\displaystyle\int_{\infty_{2j+1}}^{\infty_{2j+3}}p_{n}(x)^{2}{\rm
e}^{2\theta(x;t)}{\mathrm{d}}x=0,\quad j=0,1.$ (1.12)

Note that their linear combination in (1.11) vanishes due to the degeneracy of
the orthogonality pairing, so that this condition implies only one additional
constraint on the value of $t$. This is proven in Theorem 3.9. Another amusing
consequence is that in these cases the antiderivative of $p_{n}(x)^{2}{\rm
e}^{\frac{2x^{3}}{3}+tx}$ can be shown to be also a quasipolynomial. See
Corollary 3.11.

The second set of results involves the asymptotic analysis for large $n$ and
contains the actual proof of the Shapiro-Tate conjecture. This begins in
Section 4 where we introduce the following rescaled variables for the Shapiro-
Tater and Jimbo-Miwa cases

$\displaystyle\mbox{Shapiro-Tater},\quad z$
$\displaystyle=\hbar^{\frac{1}{3}}x,\qquad{s}=\hbar^{\frac{2}{3}}t,\qquad
E=\hbar^{\frac{4}{3}}\Lambda,\quad\hbar^{-1}=n+1$ $\displaystyle\mbox{Jimbo-
Miwa},\quad z$
$\displaystyle=\hbar^{\frac{1}{3}}x,\qquad{s}=\hbar^{\frac{2}{3}}a,\qquad\hat{b}=\hbar^{\frac{4}{3}}b,\;\;E=\frac{7s^{2}}{36}+10\hat{b},\quad\hbar^{-1}=n+\frac{1}{2}\,.$

The reason for this scaling is that it yelds the same WKB-type equation with a
small parameter $\hbar$ and a $n$-independent quartic potential

$\displaystyle\hbar^{2}\derivative[2]{y}{z}-Q(z;{s},E)y=0,$ $\displaystyle
Q(z;{s},E)=z^{4}+{s}z^{2}+2z+E\,.$

This puts both the ST and JM anharmonic oscillators on the same footing and
allows us to use the exact WKB method to compute asymptotic expressions for
the Stokes phenomenon of both systems simultaneously.

  1. 1.

In Section 4 we recall the “exact WKB analysis” following [KT05] and set up
notation used in later sections. This allows us to express the Stokes’
parameters for the quartic anharmonic oscillator in terms of the so–called
exact Vorös symbols, namely, integrals of formal series in the small parameter
$\hbar$. This section does not contain new results and is mostly a preparation
for the two subsequent ones.

  2. 2.

In Section 5.1 we use exact WKB analysis to re-derive the asymptotic
conditions for the zeros of the Vorob’ev Yablonskii polynomials already
appeared in [BM14, BM15, BB15] in a different way, based on the matching of
the Stokes’ phenomenon of the associated linear problem for the second
Painlevé equation (PII), transformed into a problem for the anharmonic
oscillator using an idea originally used by Masoero in [Mas10a, Mas10b]. This
analysis produces the asymptotic conditions in Theorem 5.2 that implicitly
describe the location of the zeros of the 1 polynomials in terms of certain
contour integrals.

  3. 3.

In Section 5.2 we derive a similar asymptotic description of the points of the
ES spectrum that correspond to multiple eigenvalues by combining our previous
results in Section 3 and the exact WKB analysis. The key result in this
section is Theorem 5.4 which leads to the quantization conditions (5.23).

  4. 4.

Finally in Section 5.3 we proceed with the comparison of the two sets of
quantization conditions that describe the two lattices of points. These are
the equations (5.23) and (5.24), that we report hereafter (to leading order):

$\displaystyle
2(n+1)\int_{\tau_{1}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=\ln\left(\frac{-1}{1+\boldsymbol{\tau}(s,E)}\right)-2i\pi(m_{1}+1)$
$\displaystyle
2(n+1)\int_{\tau_{2}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=\ln\left(-1-\frac{1}{\boldsymbol{\tau}(s,E)}\right)-2i\pi(m_{2}+1)$
$\displaystyle
2(n+1)\int_{\tau_{3}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=\ln\left({\boldsymbol{\tau}(s,E)}\right)-2i\pi(m_{3}+1)$
$\displaystyle\boldsymbol{\tau}(s,E)=\frac{\displaystyle\int_{\tau_{1}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}}{\displaystyle\int_{\tau_{2}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}},\quad\imaginary(\boldsymbol{\tau}(s,E))>0.$

where $\tau_{j}$, $j=0,1,2,3$ are the zeros of $Q(z_{+};s,E)$. The three
integers satisfy $m_{1}+m_{2}+m_{3}=n-1$ due to the fact that the sum of the
three integrals on the left is $-2i\pi(n+1)$ while the sum of the three
logarithms is $2i\pi$ (principal determination) due to the definition of
$\boldsymbol{\tau}(s,E)$. On the other hand, the quantization conditions for
the Vorob’ev-Yablonskii zeroes, to the same order of approximation, read

$\displaystyle\left(2n+1\right)\int_{\tau_{j}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=-{i\pi}-2i\pi
k_{j}$ $\displaystyle k_{1}+k_{2}+k_{3}=n-1.$

In Proposition 5.6 we show that both lattices form a regular pattern where the
local lattice generators are slowly modulated vectors. In the scaled
$s=t\hbar^{\frac{2}{3}}$ plane (where $\hbar^{-1}=n+\frac{1}{2}$ for the VY
polynomials and $\hbar^{-1}=n+1$ for the ST case) the two lattices fill a
region of triangular shape that was analyzed in [BM14, BM15, BB15, ST22] with
corners at $s_{j}={\rm e}^{\frac{i\pi}{3}j}\frac{3}{2^{1/3}}$. Within this
region the number of points is $\mathcal{O}(n^{2})$ so that the relative
spacing is $\mathcal{O}(n^{-1})$.

Our results indicate that near the origin the two lattices have a discrepancy
(in the re-scaled $s$–plane) of order $\mathcal{O}(n^{-2})$, see Theorem 5.7;
as we move away from the origin, this discrepancy progressively accumulates so
that, in a fixed small disk around $s_{0}\neq 0$, their discrepancy would be
only $\mathcal{O}(n^{-1})$. This would seem at first inconclusive because the
lattices have a separation already of order $\mathcal{O}(n^{-1})$. However, as
shown in Prop. 5.6, the lattices are very regular and have the same local
geometry. This contributes to the impression of them being almost identical.
This drift of the lattices is rather visible in Figure 2. In particular, this
settles the Conjecture of [ST22].

  5. 5.

In Appendix B we report on the numerical verifications of the quantization
conditions we have derived; moreover we observe the curious and mostly
accidental fact that if we scale the ST lattice by $n^{-\frac{2}{3}}$ instead
of the more natural $(n+1)^{-\frac{2}{3}}$ the numerical coincidence become
even more surprising. However, as we document numerically, the order of
discrepancies of the two lattices of the ST problem and VY zeros is the same
with either scaling.

Figure 2: Scaled roots of the Vorob’ev-Yablonsky polynomials
$Y_{n}(\hbar^{-2/3}s)$ with $\hbar^{-1}=n+1/2$ in red, and roots of
$D_{n}(\hbar^{-2/3}s)$ with $\hbar^{-1}=n+1$ in black, for $n=40$. This is the
natural scaling with regards to the “exact” WKB analysis, as it will be
demonstrated below.

##  2 Link to Painlevé

In this section we explain the link between the anharmonic oscillator (1.1)
and the second Painlevé transcendent, following the ideas first pioneered by
Its and Novkshenov [IN86] and later by Masoero [Mas10a, Mas10b]. The PII
equation can be expressed as the compatibility condition of two linear systems
called Lax pairs. There are several such Lax pairs but for our purposes we
will only need the Jimbo-Miwa Lax pair [JM81, FIKN06], namely:

$\begin{dcases}\partialderivative{\Phi}{x}=A(x,t)\Phi,\\\
\partialderivative{\Phi}{t}=B(x,t)\Phi,\end{dcases}$ (2.1)

with the coefficient functions

$\displaystyle
A(x,t;\alpha)=\left(x^{2}+w+\frac{t}{2}\right)\sigma_{3}+(x-u)v\sigma_{+}-\frac{2}{v}\left(xw+uw+-\alpha+\frac{1}{2}\right)\sigma_{-}$
(2.2) $\displaystyle
B(x,t)=\frac{x}{2}\sigma_{3}+\frac{v}{2}\sigma_{+}-\frac{w}{v}\sigma_{-}$
(2.3)

where $\alpha\in\mathbb{C}$ is a parameter, $u=u(t),v=v(t),w=w(t)$ are
meromorphic functions of $t$ and

$\sigma_{3}=\begin{bmatrix}1&0\\\
0&-1\end{bmatrix},\quad\sigma_{+}=\begin{bmatrix}0&1\\\
0&0\end{bmatrix},\quad\sigma_{-}=\begin{bmatrix}0&0\\\ 1&0\end{bmatrix}.$
(2.4)

The compatibility condition
$\partial_{x}\partial_{t}\Phi=\partial_{t}\partial_{x}\Phi$ of the system
(2.1) is equivalent to

$\partialderivative{A}{t}-\partialderivative{B}{x}+[A,B]=0,$ (2.5)

which in turn yields the following system of ODEs for $u,v,w$:

$\derivative{u}{t}=u^{2}+w+\frac{t}{2},\quad\derivative{v}{t}=-uv,\quad\derivative{w}{t}=-2uw+\alpha-\frac{1}{2}.$
(2.6)

Eliminating $w$ from the system $u(t)$ satisfies (PII) with parameter
$\alpha$; similarly eliminating $u(t)$ gives that $w(t)$ satisfies PXXXIV
[Inc44, FIKN06].

Now we find local solutions of the system (2.6) near a pole of the Painlevé
transcent $u(t)$. This will be useful in the upcoming computations.

######  Proposition 2.1.

Let $t=a$ be a pole of residue $+1$ of the Painlevé II function $u(t)$ with
parameter $\alpha$. Then near $t=a$ we have the following Laurent series
expansion of (2.6):

$\displaystyle u(t)$
$\displaystyle=\frac{1}{t-a}-\frac{a}{6}(t-a)-\frac{\alpha+1}{4}(t-a)^{2}+b(t-a)^{3}+\mathcal{O}{((t-a)^{4})},$
$\displaystyle w(t)$
$\displaystyle=-\frac{2}{(t-a)^{2}}-\frac{a}{3}-\frac{1}{2}(t-a)+\left(-\frac{a^{2}}{36}+b\right)(t-a)^{2}+\mathcal{O}{((t-a)^{3})},$
$\displaystyle v(t)$
$\displaystyle=c\left(\frac{1}{t-a}+\frac{a}{12}(t-a)+\frac{\alpha+1}{12}(t-a)^{2}+\mathcal{O}{((t-a)^{3})}\right),$
(2.7)

where $b$ is arbitrary and $c\neq 0$ is a constant of integration.

Similarly, if $t=a$ is a pole of residue $-1$, then near $t=a$ we have the
following Laurent expansions:

$\displaystyle u(t)$
$\displaystyle=\frac{-1}{t-a}+\frac{a}{6}(t-a)-\frac{\alpha-1}{4}(t-a)^{2}+b(t-a)^{3}+\mathcal{O}{((t-a)^{3})},$
(2.8) $\displaystyle w(t)$
$\displaystyle=\left(\frac{1}{2}-\alpha\right)(t-a)+\left(5b-\frac{a^{2}}{36}\right)(t-a)^{2}+\frac{a(2\alpha-1)}{6}(t-a)^{3}+\mathcal{O}((t-a)^{4}),$
(2.9) $\displaystyle v(t)$
$\displaystyle=c\left((t-a)-\frac{a}{12}(t-a)^{3}+\mathcal{O}((t-a)^{4}),\right)$
(2.10)

where, again, $b$ is arbitrary and $c\neq 0$ is a constant of integration.

Proof. It is well known that all poles of PII are simple and with residue $\pm
1$ (see e.g. [IN86]). For every such pole $t=a$, and arbitrary
$b\in\mathbb{C}$ there exists a solution of PII with the prescribed Laurent
expansions, where all the coefficients of order $4$ or higher are polynomials
in $a,b$ (see [GLS02]). Furthermore, the Laurent expansion for $w(t),v(t)$ can
be obtained by substituting $u(t)$ in the system (2.6) and comparing the
coefficients. $\blacksquare$

###  2.1 Gauge transformation

With a particular gauge transformation we can convert the matrix ODE
$\partialderivative{x}\Phi=A(x,t;\alpha)\Phi$ to a scalar ODE; the procedure
results in adding an apparent singularity in the equation at the position
$x=u(t)$. Thus, when the independent variable $t$ tends to one of the poles of
the solution $u(t)$ to the Painlevé equation, the singularity “escapes” to
infinity and in the limit we obtain a polynomial ODE. This gives us an
eigenvalue problem similar to the one of Shapiro and Tater (1.1).

To implement this idea in detail, first we outline the general gauge
transformation as used by [Mas10b]. Take a traceless $2\times 2$ matrix ODE
system $\Phi_{x}=M(x)\Phi$, where

$M(x)=\begin{bmatrix}m_{11}(x)&m_{12}(x)\\\ m_{21}(x)&-m_{11}(x)\\\
\end{bmatrix}.$ (2.11)

We wish to turn this into a scalar ODE of the form $y^{\prime\prime}-V(x)y=0$,
which we can do by using the gauge transformation

$W(x)=G(x)\Phi(x),\qquad G(x):=\begin{bmatrix}{m_{12}(x)}^{-\frac{1}{2}}&&0\\\
&&\\\
-\frac{m_{12}^{\prime}(x)}{2m_{12}(x)^{\frac{3}{2}}}+\frac{m_{11}(x)}{m_{12}(x)^{\frac{1}{2}}}&&m_{12}(x)^{\frac{1}{2}}\end{bmatrix}$
(2.12)

We find that $W(x)$ solves the matrix ODE $\derivative{x}W=\widehat{M}(x)W(x)$
where

$\widehat{M}(x)=G_{x}G^{-1}+GMG^{-1}=\begin{bmatrix}0&1\\\ V(x)&0\\\
\end{bmatrix},$ (2.13)

where the function $V(x)$ is:

$V(x)=m_{12}m_{21}+m_{11}^{2}+m_{11}^{\prime}-m_{11}\frac{m_{12}^{\prime}}{m_{12}}-\frac{m_{12}^{\prime\prime}}{2m_{12}}+\frac{3}{4}\left(\frac{m_{12}^{\prime}}{m_{12}}\right)^{2}.$
(2.14)

Thus the matrix ODE system $W_{x}=\widehat{M}(x)W(x)$ is then equivalent to
the scalar ODE

$y^{\prime\prime}(x)=V(x)y(x)$

and $W$ is the Wronskian matrix of a pair of independent solutions (whence the
choice of symbol).

Let us now apply this transformation to the matrix ODE from the Jimbo-Miwa Lax
pair (2.1), i.e. we set $M(x)=A(x,t)$ and perform the aforementioned gauge
transformation. Since the entries of $A(x,t)$ depend on $t$, our potential $V$
will be a function of both $x$ and $t$, namely:

$\displaystyle\begin{aligned}
V(x,t)=&\frac{1}{\left(x-u(t)\right)^{2}}\Bigg{[}2\,w\left(t\right)u(t)^{4}+\left(-4\,xw\left(t\right)-2\,\alpha+1\right)u(t)^{3}\\\
+&\left(w(t)^{2}+\left(2\,{x}^{2}+\,t\right)w(t)+x^{4}+tx^{2}+(6\alpha-1)x+\frac{{t}^{2}}{4}\right)u(t)^{2}\\\
+&\left(-2\,xw(t)^{2}+(-2\,tx+1)w(t)-2\,x^{5}-2\,tx^{3}-\frac{1}{2}\,t^{2}x-6\,\alpha\,x^{2}+\frac{t}{2}\right)u(t)\\\
+&\,x^{2}w(t)^{2}+\left(\,t{x}^{2}-x\right)w(t)+x^{6}+tx^{4}+\frac{t^{2}x^{2}}{4}+2\,\alpha\,x^{3}-\frac{tx}{2}+\frac{3}{4}\Bigg{]}\,.\end{aligned}$
(2.15)

We are interested in the expression of this potential at a pole $t=a$ of PII
equation.

######  Proposition 2.2.

Let $t=a$ be a pole with residue $\pm 1$ of the PII solution $u(t)$ with
parameter $\alpha$. Then the $\lim_{t\to a}V(x,t)$ exists and we have for a
pole of residue $+1$

$V_{\rm JM}(x;a,b,\alpha):=\lim_{t\to
a}V(x,t)=x^{4}+ax^{2}+(2\alpha-1)x+\left(\frac{7a^{2}}{36}-10b\right)$ (2.16)

and for a pole of residue $-1$

$\lim_{t\to
a}V(x,t)=x^{4}+ax^{2}+(2\alpha+1)x+\left(\frac{7a^{2}}{36}+10b\right)=V_{\rm
JM}(x;a,-b,\alpha+1)$ (2.17)

Proof. The result is obtained from a straightforward computation where we
substitute the Laurent series expansions (2.1) into (2.15) and take the limit
as $t\to a$. The details are left as exercise (it is helpful to use a computer
algebra program for this purpose). $\blacksquare$

The importance of this transformation for us is the following: the gauge
transformation (2.12) introduces only a singularity at the zeros of
$m_{1,2}(x)$, which in our case is only at $x=u(t)$, namely, the value of the
Painlevé transcendent solution $u(t)$. The singularity is clearly only a
square-root type singularity with local monodromy $-1$. Other than this, the
Stokes phenomenon of the ODE $y^{\prime\prime}(x)=V(x;t)y(x)$ is unchanged and
independent of $t$ by construction. This is evident from the fact that the
transformation of the system effected by (2.12) is a left multiplication by a
simple algebraic matrix $G(x)$ with at most square–root singularities around
which the whole matrix has a scalar monodromy of multiplication by $-1$.
However the Stokes’ phenomenon involves right multiplications of the solution
by constant matrices, which are thus the same as for the original system.

This important observation has the consequence that as $t$ approaches one of
the poles of the given solution $u(t)$ of the Painlevé equation, the
additional singularity moves off to infinity. Thus we have the following
simple but essential statement that we formalize in the proposition below.

######  Proposition 2.3.

Let $u(t)$ be a solution of the Painlevé II equation corresponding to a
particular set of Stokes data for the ODE in (2.1). Let $t=a$ be a pole of
$u(t)$ with residue $-1$ and $b$ the coefficient as in (2.8). Then the Stokes
phenomenon of the ODE

$\displaystyle
y^{\prime\prime}(x)=\left(x^{4}+ax^{2}+(2\alpha+1)x+\Lambda\right)y(x),\ \ \
\Lambda:=\frac{7a^{2}}{36}+10b$ (2.18)

is the same as the original Stokes phenomenon of the ODE (2.1).

A similar statement holds for the poles with residue $+1$ but we make the
choice of considering only those with negative residue because as formula
(1.10) shows the poles with positive residue are the zero of $Y_{n-1}(t)$ and
the poles with negative residue are the zeros of $Y_{n}(t)$.  The complete
proof of the above proposition is equivalent to the one presented in [Mas10b]
Appendix B.

##  3 A study of quasi-polynomials

In this section we find a characterization of the (quasi)-polynomials
corresponding to a repeated eigenvalue for the operator (1.4), namely, the
Shapiro-Tater eigenvalue problem:

$\displaystyle\derivative[2]{y}{x}-(x^{4}+tx^{2}+2Jx+\Lambda)y=0$ (3.1)
$\displaystyle y(s{\rm e}^{{k}\pi i/3})\to 0,\ \ \ s\to+\infty,\ {k=1,3,5}.$
(3.2)

We will call the set $(t,J,\lambda)$ for which there is a solution of the
problem (3.1)-(3.2), the Exactly Solvable (ES) spectrum. Our setting is
different from [ST22], where the authors considered a modified eigenvalue
problem (3.1) with only two boundary conditions. In this case only a finite
portion of the spectrum can be computed explicitly, and for this reason it is
a Quasi Exactly Solvable (QES) spectrum, the naming comes from [BB98].

We will see below that with three boundary conditions at infinity as in (3.2),
the whole spectrum can be characterised by the vanishing of a finite
determinant, and it is therefore Exactly Solvable. The first result in
Proposition 3.2 shows that the boundary conditions in (3.2) is compatible only
with $J$ being a positive integer. We then compute the Stokes phenomenon of
the quasi-polynomial solutions explicitly, as shown in (3.3). Finally we
relate the problem to (degenerate) orthogonality for a suitable class of non-
hermitean orthogonal polynomials, see Theorem 3.4 and Theorem 3.6 which are
the main results of this section.

######  Lemma 3.1.

The equation (3.1) admits quasi–polynomial solutions of the form

$y(x)=p(x)e^{\theta(x;t)}\quad\mbox{where
$\;\;\theta(x;t)=\frac{x^{3}}{3}+\frac{tx}{2}$},$ (3.3)

with $p(x)$ a polynomial of degree $n$ if and only if $J=n+1$ and
$\lambda=\Lambda-\frac{t^{2}}{4}$ is an eigenvalue of the operator

$\mathcal{L}_{J}:=\derivative[2]{x}+2\left(x^{2}+\frac{t}{2}\right)\derivative{x}-2(J-1)x$
(3.4)

acting on the space of polynomials of degree up to n.

Proof. Substituting $y(x)=p(x){\rm e}^{\theta(x)}$ in the ODE (3.1) gives an
equivalent differential equation for the function $p(x)$:

$\mathcal{L}_{J}(p(x))=\lambda p(x)\quad\text{ where
}\quad\lambda=\Lambda-\frac{t^{2}}{4}\,,$ (3.5)

and $\mathcal{L}_{J}$ as in (3.4). One can readily see that if $J=n+1$ then
the operator $\mathcal{L}_{n+1}$ in (3.5) preserves the space of polynomials
of degree at most $n$ and then $\Lambda$ is, by definition, an eigenvalue of
(3.1).

Viceversa, if $p(x)$ is a polynomial of degree $n$ and solves (3.5) then one
finds by inspection that the l.h.s is a polynomial of degree $n+1$ whose
leading coefficient is $2(n-J+1)$ while the r.h.s. is a polynomial of degree
$n$. Thus $J=n+1$. Then $\mathcal{L}_{n+1}$ preserves the space of polynomials
of degree $n$ and $\Lambda$ (and the corresponding $\lambda$ as per (3.5)) is
an eigenvalue of the corresponding finite dimensional operator. $\blacksquare$

###  3.1 Stokes phenomenon

Figure 3: Stokes matrices and Stokes sectors for the Shapiro-Tater eigenvalue
problem (3.1) with quasi-polynomial solutions. The Stokes matrices
$\mathbb{S}_{0},\mathbb{S}_{2},\mathbb{S}_{4}$ are all the identity.

######  Proposition 3.2.

The eigenvalue problem (3.1) with the boundary conditions (3.2) require that
$J=n+1,\ n=0,1,\dots$ and that $\lambda$ in (3.5) is an eigenvalue of the
operator $\mathcal{L}_{J}$. In particular the solutions are quasipolynomials
as in Lemma 3.1.

Proof. The equation (3.1) can be written as a first order system for the
Wronkstian matrix $\mathcal{W}(x)$:

$\displaystyle\derivative{\mathcal{W}}{x}=\begin{bmatrix}0&1\\\
x^{4}+tx^{2}+2Jx+\Lambda&0\end{bmatrix}\mathcal{W}(x),\ \
\Lambda=\lambda+\frac{t^{2}}{4}.$ (3.6)

This equation has a singularity at infinity with Poincaré rank $3$. Following
the ordinary asymptotic analysis [Wa87] we see that we have formal-series
solutions of the form

$\displaystyle\mathcal{W}_{\text{form}}(x)\sim$ $\displaystyle
x^{-\sigma_{3}}(\mathbb{I}-\sigma_{+}+\sigma_{-})\left(\mathbb{I}+\mathcal{O}(x^{-1})\right)x^{J\sigma_{3}}e^{\theta(x;a)\sigma_{3}},\quad\theta(x;t)=\frac{x^{3}}{3}+\frac{tx}{2}$
(3.7) $\displaystyle=$
$\displaystyle\begin{bmatrix}x^{J-1}+\mathcal{O}(x^{J-2})&-x^{-J-1}+\mathcal{O}(x^{-J-2})\\\
x^{J+1}+\mathcal{O}(x^{J})&x^{-J+1}+\mathcal{O}(x^{-J})\end{bmatrix}{\rm
e}^{\theta(x;t)\sigma_{3}},\quad x\to\infty.$ (3.8)

We define the Stokes sectors as shown in Figure 3, i.e. $\mathcal{S}_{k}$ is
the sector of opening $\pi/3$ centered around the rays of argument
$\pi/6+(k-1)\pi/3$:

$\mathcal{S}_{k}=\left\\{z\in\mathbb{C}:\absolutevalue{\arg(z)-\frac{\pi}{6}{-(k-1)}\frac{\pi}{3}}<\frac{\pi}{6}\right\\},\quad
k=0,1,2,3,4,5.$ (3.9)

Each open sector
$\Omega_{k}=\left\\{z\in\mathbb{C}:\absolutevalue{\arg(z)-\frac{\pi}{6}{-(k-1)}\frac{\pi}{3}}<\frac{\pi}{6}+\delta\right\\}$
with $\delta>0$, strictly contains a Stokes’ sector $S_{k}$. Then there exists
a unique solution $W^{(k)}(x)$ of (3.6) in $\Omega_{k}$ asymptotic to
$W_{\text{form}}(x)$, namely

$\mathcal{W}^{(k)}(x)\simeq\mathcal{W}_{\text{form}}(x),\quad
x\in\Omega_{k},\;\;k=0,\dots,5.$

Since $\mathcal{W}^{(k)}(x)$ and $\mathcal{W}^{(k+1)}(x)$ have the same
asymptotic expansion in $\Omega_{k}\cap\Omega_{k+1}$, there is a constant
matrix $\mathbb{S}_{k}$ (called Stokes’ matrix), such that

$\mathcal{W}^{(k+1)}(x)=\mathcal{W}^{(k)}(x)\mathbb{S}_{k},\quad
x\in\Omega_{k}\cap\Omega_{k+1},\quad k=0,1,\dots,5,$

where $\mathbb{S}_{2k}=\begin{bmatrix}1&0\\\ s_{2k}&1\end{bmatrix}$ and
$\mathbb{S}_{2k+1}=\begin{bmatrix}1&s_{2k+1}\\\ 0&1\end{bmatrix}$ for
$k=0,1,2$ (the parameters $s_{k}$ are called Stokes’ multipliers). By standard
methods one can see that the Stokes phenomenon consists of the following
relation

$\mathcal{W}^{(6)}(x)=\mathcal{W}^{(0)}(x)\mathbb{S}_{0}\mathbb{S}_{1}\dots\mathbb{S}_{5},\quad\mathcal{W}^{(6)}(x)=\mathcal{W}^{(0)}(x)e^{2\pi
iJ\sigma_{3}}$

which gives the relation

$\begin{bmatrix}1&0\\\ s_{0}&1\end{bmatrix}\begin{bmatrix}1&s_{1}\\\
0&1\end{bmatrix}\begin{bmatrix}1&0\\\
s_{2}&1\end{bmatrix}\begin{bmatrix}1&s_{3}\\\
0&1\end{bmatrix}\begin{bmatrix}1&0\\\
s_{4}&1\end{bmatrix}\begin{bmatrix}1&s_{5}\\\
0&1\end{bmatrix}\begin{bmatrix}{\rm e}^{-2i\pi J}&0\\\ 0&{\rm e}^{2i\pi
J}\end{bmatrix}=\mathbb{I}$ (3.10)

The first six matrices are the Stokes matrices associated with the directions
$\arg(x)=k\frac{\pi}{3}$, $k=0,1,..,5$ and the last matrix is the formal
monodromy matrix. The boundary conditions (3.2) imply that
$s_{0}=s_{2}=s_{4}=0$ because it means that the recessive solution along the
direction $\arg(x)=\frac{\pi}{3}$ is also recessive along the directions
$\arg(x)=k\frac{\pi}{3}$, $k=3,5$. But then the matrix equation (3.10) implies
that

$\displaystyle s_{1}+s_{3}+s_{5}=0$ (3.11)

and ${\rm e}^{2i\pi J}=1$ and hence $J$ must be an integer.

To show that $J=n+1$ is a positive integer and that $\lambda$ is an eigenvalue
of $\mathcal{L_{J}}$ we proceed as follows. Given that now the Stokes matrices
are all upper triangular, the first column of the solution is an entire
function which is asymptotic to the first column of the formal-series solution
(3.8) along all directions. But then the asymptotic (3.8) implies that the
$(1,1)$ entry is of the form $p(x){\rm e}^{\theta(x;t)}$ with $p(x)$ entire
and bounded at infinity by $x^{J-1}$. Then Liouville’s theorem implies that if
$J\geq 1$ then $p(x)$ is a polynomial, and for $J=0,-1,-2,\dots$ $p(x)$ should
vanish at infinity and hence it should be identically zero, leading to a
contradiction.

We have now established that the only solutions of the eigenvalue problem
(3.1) are quasipolynomials and therefore the hypothesis of Lemma 3.1 prevail,
thus showing that $\lambda$ is the claimed eigenvalue. $\blacksquare$

Next we consider the properties of the operator (3.1) (or equivalently of
(3.5)) under the assumption that $J=n+1\in\mathbb{N}$ and $\lambda$
(respectively $\Lambda$) is an eigenvalue. We start from further analysis of
the corresponding Stokes phenomenon.

######  Theorem 3.3.

Let $J=n+1$, $n\in\mathbb{N}$ and $\Lambda$ be an eigenvalue of the boundary
value problem (3.1)-(3.2) with eigenfunction the quasi-polynomial
$F(x)=p_{n}(x){\rm e}^{\theta(x)}$. Let $G_{k}$, be solutions of the ODE (3.1)
linearly independent from $F(x)$ and which can be expressed as

$G_{k}(x)=F(x)\int_{\infty_{k}}^{x}\frac{\differential\zeta}{F(\zeta)^{2}},\quad
k=0,2,4.$ (3.12)

Here $\infty_{k}$ indicates that the contour of integration extends to
infinity along the direction $\arg(z)={k\frac{\pi}{3}}$.

Then the Stokes phenomenon for the solutions $[F,G_{k}]$ is given by the
following equations:

$\begin{split}[F,G_{2}]=&[F,G_{0}]\begin{bmatrix}1&s_{1}\\\
0&1\end{bmatrix},\quad\quad
s_{1}:=\int_{\infty_{2}}^{\infty_{0}}\frac{\differential\zeta}{F(\zeta)^{2}}\\\
[F,G_{4}]=&[F,G_{2}]\begin{bmatrix}1&s_{3}\\\ 0&1\end{bmatrix},\quad\quad
s_{3}:=\int_{\infty_{4}}^{\infty_{2}}\frac{\differential\zeta}{F(\zeta)^{2}}\\\
[F,G_{0}]=&[F,G_{4}]\begin{bmatrix}1&s_{5}\\\ 0&1\end{bmatrix},\quad\quad
s_{5}:=\int_{\infty_{0}}^{\infty_{4}}\frac{\differential\zeta}{F(\zeta)^{2}}.\end{split}$
(3.13)

Furthermore the Stokes parameters $s_{j}$ satisfy (3.11).

Proof. Let $p_{n}(x)$ be a polynomial solution of degree $n$ of (3.5) with
$J=n+1$. We can obtain a second linearly independent solution $q$ of (3.5)
using the Wronskian identity:

$\derivative{p_{n}}{x}q-p_{n}\derivative{q}{x}=e^{-2\theta}.$ (3.14)

The solution is written as:

$q(x):=p_{n}(x)\int_{x_{0}}^{x}\left(p_{n}(\zeta)e^{\theta(\zeta;t)}\right)^{-2}\differential\zeta$
(3.15)

with $x_{0}$ is arbitrary. We have thus found two linearly independent
solutions of (3.5), and in turn we found two particular solutions of the
eigenvalue problem (3.1), namely:

$\displaystyle F(x):=$ $\displaystyle p_{n}(x)e^{\theta(x;t)},$ $\displaystyle
G_{2k}(x):=$ $\displaystyle
F(x)\int_{\infty_{2k}}^{x}\frac{\differential\zeta}{F(\zeta)^{2}},\quad
k=0,1,2.$ (3.16)

Note that since $\exp(\theta(x;t))\to 0$ along rays $\arg(x)=\pi/3,\pi,5\pi/3$
(see Fig. 4), the function $F(x)$ satisfies the boundary conditions (3.2),
while $G_{2k}(x)$ is unbounded as $x\to\infty$ along the same directions.

We can split the integral representation of $G_{2k}$ in (3.16) as follows: for
$k=0,1,2$ we have

$\displaystyle G_{2k+2}(x)=$ $\displaystyle
F(x)\int_{\infty_{2k+2}}^{x}\frac{\differential\zeta}{F(\zeta)^{2}}$ (3.17)
$\displaystyle=$ $\displaystyle
F(x)\left(\int_{\infty_{2k+2}}^{\infty_{2k}}\frac{\differential\zeta}{F(\zeta)^{2}}+\int_{\infty_{2k}}^{x}\frac{\differential\zeta}{F(\zeta)^{2}}\right)$
(3.18) $\displaystyle=$ $\displaystyle s_{2k+1}F(x)+G_{2k}(x)$ (3.19)

were the indices are taken $\mod{6}$ and we have defined $s_{2k+1},k=0,1,2$ to
be the integral of $F(\zeta)^{-2}\differential\zeta$ between $\infty_{2k+2}$
and $\infty_{2k}$, owing to the observation that its contour of integration
crosses the Stokes line of argument $(2k+1)\pi/3$ with $k=0,1,2$.

The specific contour of integration of the Stokes parameters $s_{2k+1}$
defined in (3.13) does not matter as long as it avoids the poles of the
integrand $F(\zeta)^{-2}\differential\zeta$. Indeed the integrand
$F(\zeta)^{-2}\differential\zeta$ has zero residue in the finite complex plane
because if there were a non-zero residue at a pole, the function $G_{2k}(x)$
would have non-trivial monodromy around that pole, but this is not the case
since $G_{2k}$ is a solution to the linear ODE system (3.1) which has analytic
coefficients in the finite complex plane $\mathbb{C}$. One may also verify
directly that the ODE (3.1) implies the vanishing of these residues. To prove
the condition (3.11) for the Stokes parameters $s_{j}$ we observe that sum of
the Stokes parameters corresponds to an integral of the function
$F^{-2}(\zeta)$ on a closed contour in the complex plane and by the residue
theorem it vanishes since the integrand has second order poles with zero
residue. $\blacksquare$

Figure 4: Directions at infinity $\infty_{k}$ of argument $k\frac{i\pi}{3}$
and the oriented contours $\gamma$ and $\tilde{\gamma}$ from $\infty_{1}$ to
$\infty_{3}$ and $\infty_{3}$ to $\infty_{5}$ respectively. The shaded regions
denote the sectors of dominance of ${\rm e}^{\theta(x;t)}$ where
$\theta(x,t)=\frac{x^{3}}{3}+\frac{tx}{2}$. The unshaded regions denote the
sectors of recessiveness.

###  3.2 Connection to degenerate orthogonal polynomials

In this section we characterise the quasi-polynomial solutions to (3.1) as
degenerate orthogonal polynomials in the following sense.

######  Definition 3.1.

Let $\theta(x)$ be a fixed polynomial of degree $d+1$ and positive leading
coefficient, and let $\Gamma=\sum_{j=1}^{d}s_{2j-1}\gamma_{2j-1}$, where
$\gamma_{2j-1}$ are the “wedge” contours extending from $\infty_{2j-1}$ to
$\infty_{2j+1}$ where $\infty_{\ell}$ denotes the point at infinity in the
directions $\arg(x)=\frac{\ell\pi}{d+1}$. Consider the non-hermitian
orthogonal polynomials $p_{n}(x)$ satisfying

$\langle
p_{n},x^{k}\rangle:=\sum_{j=1}^{d}s_{2j-1}\int_{\gamma_{2j-1}}p_{n}(x)x^{k}e^{2\theta(x)}\differential
x=\int_{\Gamma}p_{n}(x)x^{k}e^{2\theta(x)}\differential x=0,\quad
k=0,\dots,n-1.$ (3.20)

We say the polynomial $p_{n}$ is $l$-degenerate orthogonal if we additionally
have

$\langle p_{n},x^{k}\rangle=0\quad k=n,n+1,\dots,n+l-1.$ (3.21)

To prove that the quasi-polynomials are degenerate orthogonal polynomials we
will interpret the jumps in Theorem 3.3 as a Riemann-Hilbert problem for
orthogonal polynomials [Dei99].

######  Theorem 3.4.

Supppose that $F(x)=p_{n}(x)e^{\theta(x;t)}$ is a quasi-polynomial solution of
the ODE (3.1). Then $p_{n}(x)$ is a $1$-degenerate non-hermitian orthogonal
polynomial with respect to the weight $w(x)=e^{2\theta(x,t)}\differential x$
on the contour

$\Gamma=\varkappa\gamma+\widetilde{\varkappa}\widetilde{\gamma}$ (3.22)

where $\varkappa=s_{1}$ and $\widetilde{\varkappa}=s_{5}$ are defined in
(3.13). The contour $\gamma$ is the wedge contour from $\infty_{1}$ to
$\infty_{3}$ and $\widetilde{\gamma}$ is the wedge contour from $\infty_{5}$
to $\infty_{3}$ (see Figure 4).

Proof. Consider the quasi-polynomial $F(x):=p_{n}(x){\rm e}^{\theta(x;t)}$
which solves (3.1) with $J=n+1$. Define $G_{k}(x),k=0,2,4$ as in (3.12), which
also solves the same ODE, and let us set

$\Psi_{k}(x):=[F(x),G_{k}(x)],\quad k=0,2,4,$ (3.23)

to be a fundamental solution of (3.1). We interpret Theorem 3.3 as an $2\times
2$ Riemann-Hilbert problem solved by $\Psi_{k}(x)$ in the sectors of opening
$2\pi/3$ :

$\mathcal{S}_{k}\cup\mathcal{S}_{k+1},\quad k=0,2,4.$ (3.24)

The jump conditions (3.13) imply that $G_{k}$ satisfies:

${\rm e}^{\theta(x;t)}G_{k}(x)={\rm
e}^{\theta(x;t)}G_{k+2}(x)+s_{k+1}p_{n}(x){\rm e}^{2\theta(x;t)},\quad
k=0,2,4.$ (3.25)

The three functions $G_{k}$, $k=0,2,4$ define a piecewise–analytic function
$G(x)$ with discontinuities along the three rays $[0,\infty_{j}]$ with
$j=1,3,5$. The function $G(x)$ satisfies the jump conditions

${\rm e}^{\theta(x;t)}G_{+}(x)={\rm e}^{\theta(x;t)}G_{-}(x)+s_{j}p_{n}(x){\rm
e}^{2\theta(x;t)},\quad x\in[0,\infty_{j}],\quad j=1,3,5,$ (3.26)

where $G_{\pm}(x)$ are the boundary values of the function $G(x)$ on the right
and left boundaries of the oriented rays $[0,\infty_{j}]$, $j=1,3,5$. By the
Sokhotski–Plemelj formula we can express $G$ as the following Cauchy
transform:

$G(x):=\frac{{\rm e}^{-\theta(x;t)}}{2\pi
i}\int_{\Gamma}\frac{p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}}{\zeta-x}\differential\zeta,\quad\Gamma=s_{1}\gamma+s_{5}\tilde{\gamma},$
(3.27)

with $\gamma$ and $\widetilde{\gamma}$ as in Fig.4. Note that the two contours
overlap with the same orientation on the ray $[0,\infty_{3}]$. Note that the
function $G$ defined by (3.27) coincides with $G_{k}$ in each appropriate
sector, since they satisfy the same RHP and $G(x){\rm e}^{\theta(x)}$ and
$G_{k}(x){\rm e}^{\theta(x)}$ behave like $\mathcal{O}(x^{-1})$ at infinity.
To prove this claim it is sufficient to use the formula (3.27) for $G$ and the
asymptotics (3.8) for $G_{k}$.

This considerations imply the degenerate orthogonality of the polynomials
$p_{n}(x)$ as follows. From (3.8) with $J=n+1$ we see the Wronskian matrix of
the ODE (3.6) for the pair of solutions $(F,G_{k})$ has the following
asymptotic expansion (using ${}^{\prime}=\derivative{x}$):

$\displaystyle\begin{bmatrix}F(x)&G_{k}(x)\\\
F^{\prime}(x)&G_{k}^{\prime}(x)\end{bmatrix}\sim$
$\displaystyle\begin{bmatrix}x^{n}{\rm e}^{\theta(x)}&-x^{-n-2}{\rm
e}^{-\theta(x)}\\\ x^{n+2}{\rm e}^{\theta(x)}&x^{-n}{\rm
e}^{-\theta(x)}\end{bmatrix}\left(\mathbb{I}+\mathcal{O}(x^{-1})\right),\quad
x\to\infty.$ (3.28)

In fact, with $F$ and $G$ given as above, it follows that (3.28) is the
asymptotic expansion of $F$ and $G$ in each of the six sectors. Thus we find
that the asymptotic

$\displaystyle G(x)=$ $\displaystyle\frac{{\rm e}^{-\theta(x)}}{2\pi
i}\int_{\Gamma}\frac{p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}}{\zeta-x}\differential\zeta$ (3.29) $\displaystyle=$
$\displaystyle-\frac{1}{x}\frac{{\rm e}^{-\theta(x)}}{2\pi
i}\int_{\Gamma}p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}\left(1+\frac{\zeta}{x}+\dots+\frac{\zeta^{n+1}}{x^{n+1}}+\dots\right)\differential\zeta\sim-x^{-n-2}{\rm
e}^{-\theta(x;a)}$ (3.30)

implies the vanishing of the integrals:

$\int_{\Gamma}p_{n}(\zeta)\zeta^{k}{\rm
e}^{2\theta(\zeta;t)}\differential\zeta,\quad k=0,1,\dots,n.$ (3.31)

Therefore the polynomials $p_{n}(x)$ are degenerate orthogonal as claimed.
$\blacksquare$

One can also prove the converse of this result, but to do so we first need the
following lemma.

######  Lemma 3.5.

Consider the second order ODE in the complex plane

$\derivative[2]{y}{x}-V(x)y=0.$ (3.32)

Suppose that $x=x_{*}$ is a (possible) singularity of the potential $V(x)$
where it has a pole of order at most $2$ . Additionally, assume that $x=x^{*}$
is an apparent singularity ( namely two linearly independent solutions of the
ODE (3.32) are analytic at $x=x^{*}$). Then, $V(x)$ is analytic at $x=x^{*}$.

Proof. We argue by contradiction and consider separately the cases when $V(x)$
has a double pole and when $V(x)$ has a simple pole at $x=x^{*}$. In both
cases we make use of the indicial equation and its basic properties, which can
be found in more detail in [Olv74].

  * •

Case 1: Simple pole. Suppose that near the singular point $x=x_{*}$ the
potential is of the form

$V(x)=\frac{a}{\zeta}+b+\mathcal{O}(\zeta)$ (3.33)

with $a\neq 0$ and $\zeta=x-x_{*}$. The indicial equation of the ODE becomes

$d(d-1)=0\,.$ (3.34)

It has two solutions $d_{1}=1$ and $d_{2}=0$ differing by a non-zero integer
$d_{1}-d_{2}=1$, meaning that that there are two linearly independent
solutions $y_{1}(x),y_{2}(x)$ such that

$\displaystyle y_{1}(x)=\mathcal{O}(\zeta^{d_{1}}),$ (3.35) $\displaystyle
y_{2}(x)=\mathcal{O}(\zeta^{d_{2}})$ (3.36)

as $x\to x_{*}$ [Olv74, p.151]. But a simple computation shows that
$y_{2}(x)=1+\mathcal{O}(\zeta)$ cannot solve the differential equation (3.32)
if $V(x)$ has a simple pole at $x=x_{*}$. This is a contradiction.

  * •

Case 2: Double pole. Suppose that near the singular point $x=x_{*}$ the
potential has the shape

$V(x)=\frac{a}{\zeta^{2}}+\frac{b}{\zeta}+\mathcal{O}(1)$ (3.37)

with $a\neq 0$ and again we denote $\zeta=x-x_{*}$. The indicial equation now
gives

$\displaystyle a=$ $\displaystyle d(d-1)\ .$ (3.38)

Let $d_{1},d_{2}$ be the two solutions of the indicial equation (3.38). By the
assumption that the solutions to the ODE are all analytic, we must have that
$d_{1}$ and $d_{2}$ must be non-negative integer and not equal to each other
(if $d_{1}=d_{2}$ then one of the solutions has a logarithmic singularity, see
[Olv74]). Note that $d_{1}=d_{2}$ if and only if $a=-1/4$, in which case
$d_{1}=d_{2}=1/2$, so we may simply assume that $d_{1},d_{2}$ are non-negative
integers. Rewriting the equation $d_{1}(d_{1}-1)=d_{2}(d_{2}-1)$ as

$(d_{1}-d_{2})(d_{1}+d_{2}-1)=0$ (3.39)

we see that the only non-negative integer solutions are $d_{1}=0$ and
$d_{2}=1$ and vice versa. In either case from the indicial equation we find
that $a=0$, which is a contradiction. Therefore $V(x)$ cannot have a double
pole. $\blacksquare$

We are ready to prove the converse of Theorem 3.4, namely that these
degenerate orthogonal polynomials satisfy the Shapiro-Tater ODE (3.1).

######  Theorem 3.6.

Suppose that $p_{n}(x)$ is a $1$-degenerate orthogonal polynomial with respect
to the weight $w(x)={\rm e}^{2\theta(x;t)}$ on the weighted contour
$\Gamma:=\varkappa\gamma+\tilde{\varkappa}\tilde{\gamma}$ for
$(\varkappa,\tilde{\varkappa})\neq(0,0)$, in the notations of Theorem 3.4.
Then $F(x)=p_{n}(x){\rm e}^{\theta(x;t)}$ is a quasi-polynomial solution of
the boundary problem (3.1)-(3.2) with $J=n+1$.

Proof. Suppose that $p_{n}(x)$ is a degenerate orthogonal polynomial with
respect to the weight $w(x)={\rm e}^{\theta(x;t)}$ on the (weighted) contour
$\Gamma$, i.e.

$\int_{\Gamma}p_{n}(\zeta)z^{k}{\rm
e}^{2\theta(\zeta;t)}\differential\zeta=0\quad k=0,1,\cdots,n.$ (3.40)

Define the functions

$\displaystyle F(x)$ $\displaystyle:=p_{n}(x){\rm e}^{\theta(x;t)}$ (3.41)
$\displaystyle G(x)$ $\displaystyle:=\frac{{\rm e}^{-\theta(x;t)}}{2\pi
i}\int_{\Gamma}\frac{p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}\differential\zeta}{\zeta-x}.$ (3.42)

We claim that the Wronskian $W=FG^{\prime}-F^{\prime}G$ is constant. The
degenerate orthogonality condition implies that

$G(x)\simeq\left(h_{n}x^{-n-2}+\mathcal{O}(x^{-n-3})\right){\rm
e}^{-\theta(x;t)},\quad x\to\infty$ (3.43)

where the leading factor is

$h_{n}:=-\frac{1}{2\pi i}\int_{\Gamma}p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}\zeta^{n+1}\differential\zeta.$ (3.44)

Furthermore, by differentiating (denoting ${}^{\prime}=\derivative{x}$) we
find

$\displaystyle G^{\prime}(x)=$
$\displaystyle-\theta^{\prime}(x)G(x)+\frac{{\rm e}^{-\theta(x;t)}}{2\pi
i}\int_{\Gamma}\frac{p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}}{(\zeta-x)^{2}}\differential\zeta$ (3.45)
$\displaystyle=$
$\displaystyle-\theta^{\prime}(x)G(x)+\frac{1}{x^{2}}\frac{{\rm
e}^{-\theta(x;t)}}{2\pi i}\int_{\Gamma}p_{n}(\zeta){\rm
e}^{2\theta(\zeta;t)}\left(1+2\frac{\zeta}{x}+3\frac{\zeta^{2}}{x^{2}}+\dots\right)\differential\zeta$
(3.46) $\displaystyle=$ $\displaystyle{\rm
e}^{-\theta(x;t)}\left(-h_{n}x^{-n}+\mathcal{O}(x^{-n-1})\right)+\frac{1}{x^{2}}\frac{{\rm
e}^{-\theta(x;t)}}{2\pi i}\mathcal{O}(x^{-n-1})$ (3.47) $\displaystyle=$
$\displaystyle{\rm
e}^{-\theta(x;t)}\left(-h_{n}x^{-n}+\mathcal{O}(x^{-n-1})\right).$ (3.48)

Additionally, since $p_{n}(x)$ is a monic polynomial, we have that $F(x)\sim
x^{n}{\rm e}^{\theta(x;t)}$ and $F^{\prime}(x)\sim x^{n+2}{\rm
e}^{\theta(x;t)}$, which means that the Wronskian is bounded at infinity, i.e.
$W(x)=F^{\prime}(x)G(x)-F(x)G^{\prime}(x)=-2h_{n}+\mathcal{O}(x^{-1})$ as
$x\to\infty$.

The Wronskian $W(x)$ has no jump-discontinuities since the boundary values
$W_{\pm}(x)$ for $x\in\gamma,\widetilde{\gamma}$ satisfied

$W_{+}=F^{\prime}G_{+}-FG_{+}^{\prime}=F^{\prime}(s_{k}F+G_{-})-F(s_{k}F^{\prime}+G_{-}^{\prime})=F^{\prime}G_{-}-FG_{-}^{\prime}=W_{-}.$
(3.49)

Furthermore, since $W$ is built from locally analytic functions, it also
follows that $W(x)$ has no poles. Thus $W$ is an entire function; furthermore
it is bounded at infinity since, in each sector, $F^{\prime}G,FG^{\prime}$ are
bounded. By the Liouville theorem we conclude that $W(x)$ must be a constant,
i.e.

$W(x)=F^{\prime}(x)G(x)-F(x)G^{\prime}(x)\equiv-2h_{n}.$ (3.50)

Differentiating this equation gives that $F^{\prime\prime}/F\equiv
G^{\prime\prime}/G$. Let us denote by $V(x)$ this ratio; then both $F$ and $G$
satisfy a 2nd order linear ODE of the form:

$y^{\prime\prime}-V(x)y=0\quad\text{ with potential
}V(x):=\frac{F^{\prime\prime}(x)}{F(x)}.$ (3.51)

We can rewrite the potential using the defining expression $F=p_{n}{\rm
e}^{\theta}$ in terms of the polynomial $p_{n}$, which gives us:

$V(x)=\theta^{\prime\prime}(x)+(\theta^{\prime}(x))^{2}+2\theta^{\prime}(x)\frac{p_{n}^{\prime}(x)}{p_{n}(x)}+\frac{p^{\prime\prime}_{n}(x)}{p_{n}(x)}.$
(3.52)

Let $c$ be one of the zeros of $p_{n}$ of multiplicity $d$ then we can write
$p_{n}(x)=(x-c)^{d}h(x)$ and then we expand (3.52) near $x=c$:

$V(x)=\frac{d(d-1)}{(x-c)^{2}}+\frac{2d}{x-c}\left(\frac{h^{\prime}(c)}{h(c)}+\theta^{\prime}(c)\right)+\mathcal{O}(1),\
\ \ x\to c.$ (3.53)

This shows that $V(x)$ may have at most a second-order pole: since now $F(x)$
and $G(x)$ are both analytic near $x=c$ and both satisfy the ODE
$y^{\prime\prime}-V(x)y=0$ we deduce that all the singularities of the ODE are
apparent. We can therefore, apply Lemma 3.5, from which it follows that
$V(x)=F^{\prime\prime}/F$ is entire. We conclude that $d=1$, namely all the
zeros $x_{1},\dots,x_{n}$ of the polynomial $p_{n}(x)$ are simple and obtain

$V(x)=\theta^{\prime\prime}(x)+(\theta^{\prime}(x))^{2}+2nx+2\sum_{j=1}^{n}x_{j}+\sum_{j=1}^{n}\frac{1}{x-x_{j}}\left[\theta^{\prime}(x_{j})+\sum_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{n}\frac{1}{x_{j}-x_{k}}\right]$ (3.54)

which implies that

$\theta^{\prime}(x_{j})=\sum_{\begin{subarray}{c}k=1\\\ k\neq
j\end{subarray}}^{n}\frac{1}{x_{k}-x_{j}},\quad j=1,\dots,n.$ (3.55)

This gives the potential

$V(x)=x^{4}+tx^{2}+2x(n+1)+\frac{t^{2}}{4}+2\sum_{j=1}^{n}x_{j}$ (3.56)

which coincides with the potential in (3.1) with
$\Lambda=\frac{t^{2}}{4}+2\sum_{j=1}^{n}x_{j}$ $\blacksquare$

######  Remark 3.7.

The proof of the above theorem shows that $V(x)$ is a quartic polynomial if
and only if all zeros $x_{1},\dots,x_{n}$ of $p_{n}(x)$ are simple and satisfy
a Fekete type equilibrium property (3.55). In fact, this property also hold
true for the zeros of non-hermitean semiclassical degenerate orthogonal
polynomials; this interesting observation is generalized and expanded in
[BCG22].

To summarise, we have shown the following:

######  Corollary 3.8.

The following statements are equivalent:

  1. 1.

the boundary problem in (3.1)-(3.2) admits a quasi-polynomial solution
$F(x)=p_{n}(x){\rm e}^{\theta(x,a)}$ with $p_{n}(x)$ a polynomial of degree
$n$;

  2. 2.

The polynomial $p_{n}(x)$ (of degree $n$) is a $1$–degenerate orthogonal
polynomial with respect to the weight $w(x)={\rm e}^{2\theta(x;t)}$ on the
contour $\Gamma:=\varkappa\gamma+\tilde{\varkappa}\tilde{\gamma}$ for some
values of $(\varkappa,\tilde{\varkappa})\neq(0,0)$ and with $\gamma$ and
$\tilde{\gamma}$ as in Fig.4. Furthermore, the coefficients
$\varkappa=s_{1},\tilde{\varkappa}=s_{5}$ are expressible in terms of
$p_{n}(x)$ by the formulas (3.13).

###  3.3 Repeated eigenvalues

Recall that the characteristic polynomial $C_{n}(t,\lambda)$ in (1.6) gives us
the eigenvalues of (1.4), which in turn characterizes the Exactly Solvable
spectrum of the eigenvalue problem of (3.1) with boundary condtions (3.2).

In this section we prove the following crucial result.

######  Theorem 3.9.

The following statements are equivalent:

  1. 1.

The value of $t\in\mathbb{C}$ is such that the Exactly Solvable spectrum of
(3.1)-(3.2) has a repeated eigenvalue;

  2. 2.

The value $t\in\mathbb{C}$ is such that there exist $n\in\mathbb{N}$ and
$\Lambda\in\mathbb{C}$ so that for $J=n+1$ the problem (3.1), (3.2) has a
nontrivial solution and furthermore $D_{n}(t)=0$, with $D_{n}$ the
discriminant in (1.7);

  3. 3.

there is a quasi-polynomial solution $p_{n}(x){\rm e}^{\theta(x;t)}$ of
(3.1)-(3.2) that satisfies

$\int_{\gamma}\left(p_{n}(x)e^{\theta(x;{t})}\right)^{2}\differential
x=0,\;\;\;\int_{\widetilde{\gamma}}\left(p_{n}(x)e^{\theta(x;{t})}\right)^{2}\differential
x=0$ (3.57)

where $\gamma$ and $\widetilde{\gamma}$ are defined as in Theorem 3.6
(contours from $\infty_{1}$ to $\infty_{3}$ and $\infty_{3}$ to $\infty_{5}$
respectively, see Fig.4).

Proof. [$(1.)\Rightarrow(2.)$]. If $y(x)$ is a solution of the eigenvalue
problem (3.1) then $y$ must be a quasipolynomial according to Proposition 3.2,
and we must have that $J=n+1$ and $\Lambda$
($\lambda=\Lambda-\frac{t^{2}}{4}$) is an eigenvalue.

If this eigenvalue is repeated then the derivative of the characteristic
polynomial $C_{n}(t,\lambda)$ in (1.6) must also vanish.

[$(2.)\Rightarrow(1.)$]. This is immediate consequence of Lemma 3.1 together
with the fact that the derivative of a polynomial vanish at each root of
multiplicity higher than one.

[$(2.)\Rightarrow(3.)$]. The condition $(2.)$ means the eigenvalue has
algebraic multiplicity (at least) $2$ and so we can consider the generalized
eigenvector. Let us remind the reader here that if $v$ is the eigenvector of
the matrix $M_{n}(t)$ in (1.5) with eigenvalue $\lambda$, then the generalized
eigenvector $w$ satisfies the equation $(M-\lambda I)w=v$. It follows that the
generalized eigenvector equation for a polynomial $r(x)$ (of degree $\leq n$)
takes the form of the following differential equation for $r(x)$:

$\derivative[2]{r(x)}{x}+2\left(x^{2}+\frac{t}{2}\right)\derivative{r(x)}{x}-(2(J-1)x+\lambda)r(x)=p_{n}(x).$
(3.58)

The fact that the above non-homogeneous equation admits a polynomial solution
$r(x)$ is essential to bear in mind: such solution is not unique since we can
add to it an arbitrary multiple of $p_{n}(x)$. The associated homogeneous
differential equation has two solutions, one of which is $p_{n}(x)$ (a
polynomial of degree $n=J-1$) and the other one is obtained from (3.12) and
can be written as follows:

$q_{k}(x)=p_{n}(x)\int_{\infty_{k}}^{x}F(\zeta)^{-2}\differential\zeta,\ \ \
F(x):=p_{n}(x){\rm e}^{\theta(x;{t})},$

where we can choose any $k=0,2,4$ for the basepoint of integration.
Consequently, a particular solution $r_{0}(z)$ of (3.58) is found by the
standard “variation of parameters” as follows:

$r_{0}(x):=p_{n}(x)\underbrace{\int_{\infty_{0}}^{x}F(\zeta)^{-2}\left(\int_{\infty_{1}}^{\zeta}F(s)^{2}\differential
s\right)\differential\zeta}_{H(x)}.$ (3.59)

Before proceeding with the proof we claim and prove that $H(x)$ is bounded by
a constant as $x\to\infty_{1}$. To see this it is sufficient to show that the
inner integral in the defintion of $H(x)$ tends to zero as
$\mathcal{O}(\zeta^{2n-2}){\rm e}^{2\theta(\zeta;t)}$ (recall that the real
part of $\theta$ goes to $-\infty$ along the direction $\infty_{1}$ and hence
this is exponentially small). To prove this claim we write
$p_{n}^{2}(s)=2Q_{2n-2}(s)\theta^{\prime}(s)+R(s)$, with $Q_{2n-2}$ the
quotient (of degree $2n-2$ and $R$ the remainder (of degree $1$) polynomials
of division by $2\theta^{\prime}$. Then integrating by parts we obtain

$\displaystyle\int_{\infty_{1}}^{\zeta}F(s)^{2}{\mathrm{d}}s=Q_{2n-2}(\zeta){\rm
e}^{2\theta(\zeta)}+\int_{\infty_{1}}^{\zeta}R(s){\rm
e}^{2\theta(s)}\differential s.$ (3.60)

To estimate the last integral we “force” again an integration by parts:

$\displaystyle\int_{\infty_{1}}^{\zeta}R(s){\rm e}^{2\theta(s)}\differential
s=$ $\displaystyle\frac{R(\zeta)}{2\theta^{\prime}(\zeta)}{\rm
e}^{2\theta(\zeta)}-\int_{\infty_{1}}^{\zeta}\left(\frac{R(s)}{2\theta^{\prime}(s)}\right)^{\prime}{\rm
e}^{2\theta(s)}\differential s.$ (3.61)

In the remaining integral the prefactor to the exponential is
$\mathcal{O}(s^{-2})$ and thus the integral is easily estimated to be
$\mathcal{O}(1){\rm e}^{2\theta}$. With the claim proved, let us proceed: we
know that (3.58) has by assumption a polynomial solution and we are going to
show now that $r_{0}(x)$ itself is such a polynomial. Indeed, the general
solution of (3.58) is obtained by adding an arbitrary linear combination of
$p_{n}(x),q_{0}(x)$ to $r_{0}(x)$ in (3.59). We then must show that here are
constants $A,B$ such that $r_{0}(x)+Ap_{n}(x)+Bq_{0}(x)$ is a polynomial of
degree at most $n$; clearly here only the value of $B$ is relevant (since
$p_{n}$ is already a polynomial). Thus the issue boils down to whether there
is a value of $B$ for which $r_{0}(x)+Bq_{0}(x)$ is a polynomial. We first
observe that the only possible value for $B$ must be zero. Indeed consider the
asymptotic behaviour for $x\to\infty_{1}$. The integral $H(x)$ is bounded by a
constant as $x\to\infty_{1}$. However, $q_{0}(x)$ has a dominant exponential
growth in this direction and hence the expression $r_{0}(x)+Bq_{0}(x)$ has
polynomial growth in this direction only if $B=0$.

We thus conclude that since (3.58) has a (pencil of) polynomial solution(s),
this must be given by $r_{0}(x)+Ap_{n}(x)$ so that $r_{0}$ must itself be a
polynomial.

We now establish (3.57). To this end we consider the behaviour of $r_{0}(x)$
near $\infty_{3}$ and $\infty_{5}$. We can write, for example for
$\infty_{3}$,

$\displaystyle r_{0}(x)$
$\displaystyle=p_{n}(x)\int_{\infty_{0}}^{x}F(\zeta)^{-2}\left(\int_{\infty_{1}}^{\infty_{3}}F(s)^{2}\differential
s+\int_{\infty_{3}}^{\zeta}F(s)^{2}\differential s\right)\differential\zeta$
$\displaystyle=\left(\int_{\infty_{1}}^{\infty_{3}}F(s)^{2}\differential
s\right)q_{0}(x)+p_{n}(x)\int_{\infty_{0}}^{x}F(\zeta)^{-2}\left(\int_{\infty_{3}}^{\zeta}F(s)^{2}\differential
s\right)\differential\zeta$ (3.62)

The second term above is polynomially bounded near $\infty_{3}$, by the same
argument used to show that $H$ is bounded near $\infty_{1}$. But since $q_{0}$
is exponentially dominant also near $\infty_{1}$ we deduce the condition
$\int_{\infty_{1}}^{\infty_{3}}F(s)^{2}\differential s=0$. One similarly
deduces $\int_{\infty_{5}}^{\infty_{3}}F(s)^{2}\differential s=0$, which
establishes (3.57).

[$(3.)\Rightarrow(2.)$]. Consider the expression (3.59). It is easy to see
directly that it satifies the generalized eigenvector equation (3.58); we must
only verify that the conditions (3.57) guarantee that $r_{0}(z)$ is a
polynomial. But this follows again from the Liouville theorem and using (3.3).
$\blacksquare$

######  Remark 3.10.

The Theorem 3.9 seems at first sight nothing short of a miracle; indeed once
we fix $J=n+1\in\mathbb{N}$, then the ODE (3.1) has only two continuous
parameters $t,\lambda$.

However the multiple eigenvalue condition apparently involves now three
equations which are:
(i) the existence of a quasi-polynomial solution which determines $\lambda$ as
a function of $t$;
(ii) the two equations (3.57) for the parameter $t$.

However, the system is actually not overdetermined because of the following
reasoning: if, for given $J=n+1\in\mathbb{N}$ the pair $(t,\Lambda)$ is in the
ES spectrum, then we have shown in Theorem 3.4 that $p_{n}(x)$ is a degenerate
orthogonal polynomial, namely,

$\displaystyle\varkappa\int_{\gamma}p_{n}^{2}(x){\rm
e}^{2\theta(x;t)}{\mathrm{d}}x+\widetilde{\varkappa}\int_{\widetilde{\gamma}}p_{n}^{2}(x){\rm
e}^{2\theta(x;t)}{\mathrm{d}}x=0,$ (3.63)

with $\gamma,\widetilde{\gamma},\varkappa,\widetilde{\varkappa}$ defined in
the same theorem. The coefficients $\varkappa,\widetilde{\varkappa}$ cannot be
both vanishing for otherwise the Stokes phenomenon of the ODE would be trivial
(which is not possible). Then the two equations (3.57) yield only one
additional constraint.

The following corollary is also a nice property, although it will not be used
in the rest.

######  Corollary 3.11.

Suppose that $y(x)=p_{n}(x)^{\theta(x;t)}$ is a solution of the boundary
problem (3.1) and hence $t$, $J=n+1$ and $\Lambda$ are in the ES spectrum.
Suppose also that $\Lambda=\lambda+\frac{t^{2}}{4}$ is a repeated eigenvalue.
Then the antiderivative of $y^{2}(x)$ is also a quasipolynomial,

$\displaystyle\int^{x}p_{n}(\zeta)^{2}{\rm
e}^{2\theta(\zeta;t)}{\mathrm{d}}\zeta=q_{2n+2}(x){\rm e}^{2\theta(x;t)}.$
(3.64)

Proof. Take the integration basepoint from $\infty_{1}$ and consider the
function

$\displaystyle K(x)=\int_{\infty_{1}}^{x}p_{n}(z)^{2}{\rm
e}^{2\theta(z;t)}{\mathrm{d}}z$ (3.65)

Standard asymptotic analysis shows that $K(x)$ is of the form
$\mathcal{O}(x^{2n+2}){\rm e}^{2\theta(x;t)}$ along the (dominant) direction
towards $\infty_{0},\infty_{2},$ and $\infty_{4}$.

Moreover towards $\infty_{1}$ we have

$\displaystyle q(x):=K(x){\rm e}^{-2\theta(x;t)}=\mathcal{O}(x^{2n+2}).$
(3.66)

Thus the function $q(x)=K(x){\rm e}^{-2\theta(x;t)}$ is polynomially bounded
along the sectors containing $\infty_{j},$ $j\in\\{0,1,2,4\\}$. If we can show
that it is also polynomially bounded in the sectors containing $\infty_{3}$
and $\infty_{5}$ then we conclude, by the Liouville’s theorem, that it is a
polynomial.

Consider $x$ tending to infinity along the sector that contains $\infty_{3}$;
we have

$\displaystyle q(x)$ $\displaystyle={\rm
e}^{-2\theta(x;t)}\int_{\infty_{1}}^{x}p_{n}(\zeta)^{2}{\rm
e}^{2\theta(\zeta;t)}{\mathrm{d}}\zeta=$ $\displaystyle={\rm
e}^{-2\theta(x;t)}\int_{\infty_{3}}^{x}p_{n}(\zeta)^{2}{\rm
e}^{2\theta(\zeta;t)}{\mathrm{d}}\zeta+{\rm
e}^{-2\theta(x;t)}\int_{\infty_{1}}^{\infty_{3}}p_{n}(\zeta)^{2}{\rm
e}^{2\theta(\zeta;t)}{\mathrm{d}}\zeta.$ (3.67)

The second term above is just a multiple of ${\rm e}^{-2\theta(x;t)}$ which is
not polynomially bounded near $\infty_{2j+1}$. However the coefficient is
precisely one of the integrals in (3.57) which has been proved to vanish. Thus
$q(x)$ is also polynomially bounded near $\infty_{3}$. Similarly we can prove
that it is polynomially bounded at $\infty_{5}$ and the proof is thus
complete. $\blacksquare$

##  4 Exact WKB analysis

In this section we give a brief overview of the exact WKB method and establish
useful notations for later. The material here is well-known and we follow the
exposition in [IN14] and [KT05]. In order to apply the exact WKB method to
(1.1) and (2.17) we need to scale them appropriately to obtain an ODE of the
form $y^{\prime\prime}-\hbar^{-2}Qy=0$, where $\hbar$ is a small parameter.

We introduce the scaling for the Shapiro-Tater potential

$V_{\rm ST}=x^{4}+tx^{2}+2(n+1)x+\Lambda$ $\displaystyle z$
$\displaystyle=(n+1)^{-1/3}x,\qquad{s}=(n+1)^{-2/3}t,\qquad
E=(n+1)^{-4/3}\Lambda$ (4.1)

which implies $\derivative{y}{x}=(n+1)^{-2/3}\derivative{y}{z}$.

A similar scaling applies in the case of the Jimbo–Miwa potential (for the
poles with residue $-1$),

$V_{\rm JM}(x;a,b)=x^{4}+ax^{2}+(2n+1)x+\left(\frac{7a^{2}}{36}+10b\right)$

$\displaystyle
z=\left(n+\frac{1}{2}\right)^{-1/3}x,\qquad{s}=\left(n+\frac{1}{2}\right)^{-2/3}a,\qquad\widehat{b}=\left(n+\frac{1}{2}\right)^{-4/3}b.$
(4.2)

The result of these scaling is the exact same potential with the
identification of $E=\frac{7{s}^{2}}{36}+10\widehat{b}$, but with different
scaling factors. If we set

$\displaystyle Q(z;{s},E)=z^{4}+{s}z^{2}+2z+E,$ (4.3)

then in either cases we obtain an $n$-independent potential, namely:

$\displaystyle\derivative[2]{y}{z}-\left(n+1\right)^{2}Q(z;{s},E)y=0,\quad\text{for
the Shapiro-Tater case;}$ (4.4)
$\displaystyle\derivative[2]{y}{z}-\left(n+\tfrac{1}{2}\right)^{2}Q(z;{s},E)y=0,\quad\text{for
the Jimbo--Miwa case.}$ (4.5)

###  4.1 Schrödinger, Riccati and WKB

Keeping with the tradition we will denote by $\hbar^{-1}$ the large parameter
in either (4.5), (4.4); namely $\hbar=(n+1)^{-1}$ in the Shapiro-Tater (ST)
case or $\hbar=(n+1/2)^{-1}$ in the Jimbo-Miwa (JM) case.

Consider the following Schrödinger equation with small parameter $\hbar$ and
polynomial potential $Q(x)$ (4.3):

$\derivative[2]{y}{z}-\hbar^{-2}Q(z;{s},E)y=0.$ (4.6)

z We now explain how to construct the WKB solutions of (4.6). The formal
series ansatz

$y(z)=\exp\left(\int^{z}S(\zeta,\hbar)\differential\zeta\right),$ (4.7)

with

$S(z,\hbar):=h^{-1}S_{-1}(z)+\sum_{k\geq 0}\hbar^{k}S_{k}(z)$ (4.8)

implies that $S(z,\hbar^{-1})$ satisfies the following Riccati equation:

$S^{2}(z)+\derivative{S(z)}{z}=\hbar^{-2}Q(z).$ (4.9)

Comparing each power of $\hbar$ we get a recursive relation for the
coefficients $S_{k}(z)$ of the formal series:

$\displaystyle S_{-1}(z)^{2}=$ $\displaystyle\quad Q(z),$ (4.10)
$\displaystyle 2S_{-1}S_{k+1}=$
$\displaystyle-\derivative{S_{k}}{z}-\sum_{\begin{subarray}{c}n+m=k\\\ n,m\geq
0\end{subarray}}S_{n}S_{m}\quad(k\geq-1).$ (4.11)

The first three coefficients $S_{k}(z)$ are:

$\displaystyle S_{-1}(z)=\sqrt{Q(z)},\qquad
S_{0}(z)=-\frac{1}{4}\frac{Q^{\prime}(z)}{Q(z)},\qquad
S_{1}(z)=\frac{4Q(z)Q(z)^{\prime\prime}-5(Q(z)^{\prime})^{2}}{32Q(z)^{\frac{5}{2}}}.$
(4.12)

The choice of sign for $S_{-1}(z)=\pm\sqrt{Q(z)}$ gives two solutions of the
Riccati equation (4.9), denoted $S_{+},S_{-}$ respectively; if we change sign
in $S_{-1}$ then all the odd $S_{2k+1}$ change sign while the even $S_{2k}$
remain unchanged. These correspond to the two linearly independent solutions
of (4.6). We define the odd and even parts of $S(z;\hbar)$:

$\displaystyle S_{\rm{odd}}(z;\hbar):=$
$\displaystyle\frac{1}{2}\big{(}S_{+}(z,\hbar)-S_{-}(z,\hbar)\big{)}$ (4.13)
$\displaystyle S_{\rm{even}}(z;\hbar):=$
$\displaystyle\frac{1}{2}\big{(}S_{+}(z,\hbar)+S_{-}(z,\hbar)\big{)}$ (4.14)

Since our potential $Q(z)$ is independent of $\hbar$, it follows that
$S_{\rm{odd}}$ only contains odd powers of $\hbar$, namely:

$\displaystyle S_{\rm{odd}}(z;\hbar)=\frac{1}{\hbar}\sqrt{Q(z)}+\hbar
S_{1}(z)+\hbar^{3}S_{3}(z)+\dots.$ (4.15)

We will only care about $S_{\rm{odd}}$ since it can be shown [KT05] that
$S_{\text{even}}$ can be written in terms of $S_{\rm{odd}}$ as:

$S_{\rm even}(z,\hbar)=-\frac{1}{2}\derivative{x}\log S_{\rm{odd}}(z,\hbar).$
(4.16)

Another important fact is the following: the differential
$S_{\rm{odd}}(z){\mathrm{d}}z$ has a pole at $z=\infty$ which comes solely
from the term $S_{-1}(z)=\sqrt{Q}(z)$; namely,
$S_{2k+1}(z)=\mathcal{O}(z^{-2})$ as $z\to\infty$. This property can also be
shown inductively.

These facts motivate the following definition. As is customary in exact WKB
theory, we will refer to the zeroes
$\mathrm{T}=\\{\tau_{0},\tau_{1},\tau_{2},\tau_{3}\\}$ of $Q(z)$ as turning
points.

######  Definition 4.1 (WKB solutions).

The  WKB solutions to (4.6) are formal power series in $\hbar$ given in terms
of $S_{\rm{odd}}$ in (4.15). We give two different normalizations that will be
used throughout this paper.

  * •

Near a turning point $\tau$ of the potential $Q(z)$ we define the normalized
WKB solutions to be:

$\psi^{(\tau)}_{\pm}(z,\hbar):=\frac{1}{\sqrt{S_{\rm{odd}}(z,\hbar)}}\exp\left(\pm\int^{z}_{\tau}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta\right)$
(4.17)

  * •

Near infinity we define the normalized WKB solutions to be:

$\psi_{\pm}^{(\infty)}(z,\hbar):=\frac{1}{\sqrt{S_{\rm{odd}}(z,\hbar)}}\exp\Big(\pm
R(z;\hbar)\Big{missing})$ (4.18)

where

$\displaystyle R(z;\hbar):=$
$\displaystyle\frac{1}{\hbar}\lim_{p\to\infty}\left[\int_{p}^{z}\sqrt{Q(\zeta)}\differential\zeta-\left(\frac{p^{3}}{3}+\frac{s}{2}p+\log
p\right)\right]+\sum_{j\geq
0}\hbar^{2j+1}\int_{\infty}^{z}S_{2j+1}(\zeta){\mathrm{d}}\zeta$
$\displaystyle=$
$\displaystyle\frac{1}{\hbar}\left(\frac{z^{3}}{3}+\frac{s}{2}z+\log(z)\right)+\mathcal{O}(z^{-1})\mathcal{O}(\hbar),\quad(z\to\infty),$
(4.19)

where all logarithms are principal. Notice that $R(z;\hbar)$ is the anti-
derivative of $S_{\rm{odd}}(z,\hbar)$ that does not have a constant term in
the expansion as $|z|\to\infty$.

######  Remark 4.1.

The integral in the exponent of $\psi_{\pm}^{(\tau)}$ is to be understood
term-wise in each coefficient of the powers of $\hbar$. Additionally,
$S_{\rm{odd}}(z,\hbar)$ is multivalued on $\mathbb{C}$ with branch points at
the zeros $\mathrm{T}:=\\{\tau_{0},\tau_{1},\tau_{2},\tau_{3}\\}$ of $Q(z)$.
Therefore the integral should be considered on the elliptic compact Riemann
surface $\overline{\Sigma}$ obtained from the affine curve

$\Sigma=\Big{\\{}(w,z)\in\mathbb{C}^{2}:w^{2}=Q(z;s,E)\Big{\\}}$ (4.20)

by adding two points $P_{\infty}^{\pm}$ at infinity. The projection
$\pi:\overline{\Sigma}\mapsto\overline{\mathbb{C}}$ from the Riemann surface
$\overline{\Sigma}$ to the extended complex plane $\overline{\mathbb{C}}$,
maps $(w,z)\mapsto z$. The projection $\pi$ realizes $\overline{\Sigma}$ as a
double cover of $\overline{\mathbb{C}}$ ramified at the zeros of $Q$ (turning
points). The pre-image of any point $z\in\overline{\mathbb{C}}$ are the two
points $\pi^{-1}(z)=(z,\pm w)$ on the two sheets of the Riemann surface where
the numbering of the sheets is such that $(z,w=\sqrt{Q(z)})$ belongs to the
first sheet. Choosing branch cuts and the first sheet of the Riemann surface,
we can talk about the integral from a turning point $\tau$ to $z$ by defining:

$\int_{\tau}^{z}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta:=\frac{1}{2}\int_{\gamma(z)}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta$
(4.21)

where the contour $\gamma(z)$ lives in the Riemann surface $\overline{\Sigma}$
and goes from $\widehat{z}$ to $z$, where $\widehat{z}$ lies on the second
sheet and the second sheet is reached by encircling the branch point $\tau$.
This integral is a well defined convergent integral. Finally, it can be shown
by induction that all the correction terms $S_{2j+1}(z){\mathrm{d}}z$, $j\geq
0$ are differentials on the Riemann surface $\overline{\Sigma}$ with poles of
increasing order at each of the ramification points, but always without
residue. For more details see [IN14].

It is useful to relate the periods of $S_{1}$ to the periods of $S_{-1}$. To
this end we have the following

######  Proposition 4.2.

If $I({s},E)$ denotes the period of $S_{-1}=\sqrt{Q}$ along a closed contour
$\gamma$, then the corresponding period of $S_{1}$ is

$\displaystyle\oint_{\gamma}S_{1}(z)\differential
z=\left(-\frac{\partial^{2}}{\partial{s}\partial
E}-\frac{s}{6}\frac{\partial^{2}}{\partial E^{2}}\right)I({s},E)$ (4.22)

Proof. Using (4.3) we see that

$\displaystyle\partial_{E}\sqrt{Q(z;{s},E)}=\frac{1}{2\sqrt{Q(z;{s},E)}}.$
(4.23)

Now, we can write (we drop the indication of the dependence on $x,{s},E$)

$\displaystyle
S_{1}=\frac{1}{48}\frac{Q^{\prime\prime}}{Q^{\frac{3}{2}}}-\frac{5}{24}\left(\frac{1}{\sqrt{Q}}\right)^{\prime\prime}.$
(4.24)

The periods of the second term in (4.24) vanish because this gives an exact
differential; the first term reads

$\displaystyle\frac{1}{48}\frac{Q^{\prime\prime}}{Q^{\frac{3}{2}}}=\frac{6z^{2}+{s}}{24Q^{\frac{3}{2}}}=-(6\partial_{s}+{s}\partial_{E})\left(\frac{1}{12Q^{\frac{1}{2}}}\right)=-(6\partial_{s}+{s}\partial_{E})\partial_{E}\left(\frac{1}{6}Q^{\frac{1}{2}}\right).$
(4.25)

Integrating (4.24) along $\gamma$ and using the identity (4.25) completes the
proof. $\blacksquare$

###  4.2 Stokes graphs and connection formulæ.

The WKB series are asymptotic to actual solutions of (4.6) in certain regions
of the complex plane that we presently define. We start making by fixing a
choice of $\sqrt{Q(z)}$ and introducing the notion of Stokes’ graph.

######  Definition 4.2 (Square root of $Q$).

Choosing the branch cuts $\mathcal{B}$ of $\sqrt{Q(z)}$ in the finite part of
the complex plane, the function $\sqrt{Q(z)}$ becomes single valued in the
complement of the branch cuts and we fix it in such a way that
$\sqrt{Q(z)}\sim z^{2}$ as $|z|\to\infty$. This choice identifies the first
sheet of the Riemann surface $\overline{\Sigma}$. The second sheet correspond
to the choice $\sqrt{Q(z)}\sim-z^{2}$ as $|z|\to\infty$.

To minimize confusion when performing integrations along the branch cuts we
will make the following explicit definition.

######  Definition 4.3 (Branch cut integration).

Let $\tau$ and $\widetilde{\tau}$ be two zeroes of $Q(z;s,E)$ joined by a
branch cut. We denote by

$\int_{\tau}^{\widetilde{\tau}}\sqrt{Q(z_{+};s,E)}\differential z$ (4.26)

to be the integral of $\sqrt{Q(z;s,E)}$ along the $+$ side of the branch cut
oriented from $\tau$ to $\widetilde{\tau}$. We denote with a minus sign $-$
the corresponding integral along the $-$ side of the branch cut. As usual, the
$+$ and $-$ sides correspond to left side and right side of the oriented
contour, respectively.

Now we introduce the notion of Stokes’ curve, which we will use to build
Stokes’ graphs.

######  Definition 4.4 (Stokes’ curve).

A Stokes curve of the potential $Q(z)$ is a horizontal trajectory of the
quadratic differential $Q(z)\differential z^{2}$ where one of the end-points
is a turning point. In other words, in a local coordinate $z$ of
$\overline{\Sigma}$ it is a curve emanating from a turning point $\tau$ and
satisfying

$\imaginary\int_{\tau}^{z}\sqrt{Q(\zeta)}\differential\zeta=0.$ (4.27)

The support of the curve is independent of choice of determination of
$\sqrt{Q}$. Furthermore:

  1. 1.

The orientation of a Stokes curve is defined by the direction in which
$\real\int_{\tau}^{z}\sqrt{Q(u)}\differential u$ is increasing.

  2. 2.

The directions near the point at $\infty$ are indicated by $\oplus$ or
$\ominus$ if the function $\real\int_{\tau}^{z}\sqrt{Q(u)}\differential u$ is
increasing or decreasing, respectively, along the corresponding Stokes curve.
We say the Stokes curve is oriented towards $\oplus$ or oriented away from
$\ominus$, respectively.

######  Definition 4.5 (Stokes graph).

The Stokes’ graph $\mathcal{G}$ associated to the potential $Q(z)$ is the
graph embedded in $\overline{\mathbb{C}}$ where the edges are the Stokes
curves and the vertices are the turning points $\mathrm{T}=\\{\tau_{j}\\}$ of
the potential and the directions near the point at infinity
$\infty_{j},j=0,\dots,6$. The graph is oriented according to definition 4.4.
The Stokes curves are indicated with black lines in Fig. 5 and with blue and
red lines in Fig.6.

Furthermore, we decorate the Stokes graph with the following additional
contours.

  1. 1.

Branch cuts. We draw additional paths called branch-cuts between all pairs of
turning points in such a way that they do not intersect any Stokes curve.
Their orientation is fixed in an arbitrary way that we shall specify in each
case. These branch cuts are indicated by green lines in Fig. 5 for each
possible Stokes graph configuration.

  2. 2.

Ideal paths. We draw arbitrary (smooth) paths connecting the different
$\infty$ in all possible ways that do not intersect any of the Stokes curves.
These paths are indicated by dashed lines in Fig. 5 and determine an ideal
triangulation of an ideal hexagon. We call these paths ideal paths and the
resulting partition of the plane the ideal triangulation. We orient the outer
ideal paths that form the hexagon in the clockwise way (see Fig. 6), while the
remaining inner ideal paths are oriented towards the $\oplus$ directions. Note
that the ideal paths separating Stokes regions intersect one and only one
branch-cut, and the orientation of an ideal path on the two sides of a branch
cut is opposite. Refer to Fig. 6.

  3. 3.

Stokes’ regions and external regions. Consider the connected components of the
complement of the Stokes graph, the ideal paths, and the branch-cuts. Amongst
them the components that have at least one Stokes curve on the boundary will
be called Stokes’ regions. The remaining ones (unbonded regions bounded by a
ideal path) will be called external regions 444These are the regions on the
outside of the ideal “hexagon” in Fig. 5 and Fig. 6..

The construction in Def. 4.5, under the Assumptions 4.5, is crafted so that
each Stokes region $\mathscr{D}$ has precisely one and only one turning point
on its boundary.

In each Stokes region $\mathscr{D}$ one can select a fundamental basis of
solutions of the ODE (4.6) as we now explain. Consider a Stokes curve $\gamma$
originating at the turning point $\tau$ and oriented towards $\oplus$ (with
similar considerations for the $\ominus$ curves ): since the real part of
$\int^{z}_{\tau}\sqrt{Q(u)}du$ is increasing, there is a region around
$\gamma$ near $\infty$ where the real part is positive; then the formal
solution $\psi_{-}^{(\tau)}$ is recessive (i.e. exponentially small as
$\hbar\to 0_{+}$). Then there is a unique solution $\Psi(z;\hbar)$ which is
asymptotic to this $\psi_{-}^{(\tau)}$ in both Stokes regions on the two sides
of $\gamma$. If $\mathscr{D}$ is one of these regions, we will denote this
first selected solution by $\Psi^{(\mathscr{D})}_{-}$. To uniquely determine
the other solution $\Psi^{(\mathscr{D})}_{+}$, it is not sufficient to examine
its asymptotic behaviour near $\gamma$ because its asymptotics there is
dominant (i.e. exponentially large as $\hbar\to 0_{+}$). However, the same
Stokes region must be bounded also by either a $\ominus$ trajectory or a
branch-cut. In the former case, in a neighbourhood of the $\ominus$ trajectory
the formal solution $\psi_{+}^{(\tau)}$ is now recessive and this allows to
uniquely determine $\Psi^{(\mathscr{D})}_{+}$.

If, instead, the other boundary is a branch-cut, we need to consider the
Stokes region, $\mathscr{D}^{\prime}$ on the other side of the cut: in this
region, due to having crossed the branch-cut, the formal solution
$\psi_{+}^{(\tau)}$ is now somewhere recessive and this allows to fix
$\Psi^{(\mathscr{D})}_{+}$ uniquely. We refer to [BK19], Section 5 for more
details. The following theorem of Vöros [Vör83], [KT05] relates the WKB
solutions of different Stokes regions near the same turning point.

######  Theorem 4.3 ( [Vör83], [KT05]).

Let $\mathscr{D}$ be a Stokes region. Then there exist unique solutions
$\Psi^{(\mathscr{D})}_{\pm}$, of (4.6), that we refer to as normalized
solutions, that are asymptotic to the WKB solutions in Def. 4.1 uniformly in
$z$ in the Stokes region, that is:

$\Psi^{(\mathscr{D})}_{\pm}(z)\sim\psi_{\pm}^{(\tau)}(z;\hbar)\quad\hbar\to
0,\quad x\in\mathscr{D}.$ (4.28)

Furthermore, let $\mathscr{D}_{\ell},\mathscr{D}_{r}$ be two adjacent Stokes
regions separated by the Stokes curve $\gamma$ oriented as in Def. 4.5 with
$\mathscr{D}_{\ell}$ on the left and $\mathscr{D}_{r}$ on the right of
$\gamma$. Then the corresponding solutions
$\Psi_{\pm}^{\mathscr{D}_{\ell}},\Psi_{\pm}^{\mathscr{D}_{r}}$ are related by:

$\begin{split}&\left[\Psi^{(\mathscr{D}_{\ell})}_{+}(z),\Psi^{(\mathscr{D}_{\ell})}_{-}(z)\right]=\left[\Psi^{(\mathscr{D}_{r})}_{+}(z),\Psi^{(\mathscr{D}_{r})}_{-}(z)\right]\begin{cases}B\quad\text{
if $\gamma$ is oriented towards $\oplus$}\\\\[10.0pt] R\quad\text{ if $\gamma$
is oriented away from $\ominus$},\end{cases}\\\ &B:=\begin{bmatrix}1&0\\\
-i&1\end{bmatrix},\quad R:=\begin{bmatrix}1&i\\\ 0&1\end{bmatrix}\end{split}$
(4.29)

The relationship between solutions in regions separated by an ideal path is
simply a different scaling as given by the following proposition.

######  Proposition 4.4.

Let $\mathscr{D}_{\ell},\mathscr{D}_{r}$ be two Stokes regions separated by an
ideal path $\sigma$, with $\mathscr{D}_{\ell}$ on the left and
$\mathscr{D}_{r}$ on the right of $\sigma$ with the orientation given in Def.
4.5. Let $\tau_{\ell},\tau_{r}$ be the (unique) turning points on the
boundaries of $\mathscr{D}_{\ell},\mathscr{D}_{r}$, respectively. Then we have
the following connection formula

$[\psi_{+}^{(\tau_{\ell})},\psi_{-}^{(\tau_{\ell})}]=[\psi_{+}^{(\tau_{r})},\psi_{-}^{(\tau_{r})}]\exp(\sigma_{3}v_{\ell
r})$ (4.30)

where

$v_{\ell r}=v_{\ell
r}(\hbar):=\int_{\tau_{\ell}}^{\tau_{r}}S_{\rm{odd}}(z_{+},\hbar)\differential
z$ (4.31)

with the integration along the branch cut according to Def. 4.3.

The corresponding actual solutions
$\Psi^{(\mathscr{D}_{\ell})}_{\pm},\Psi^{(\mathscr{D}_{r})}_{\pm}$ given in
Thm. 4.3 are similarly related:

$\displaystyle[\Psi^{(\mathscr{D}_{\ell})}_{+},\Psi^{(\mathscr{D}_{\ell})}_{-}]=[\Psi^{(\mathscr{D}_{r})}_{+},\Psi^{(\mathscr{D}_{r})}_{-}]\exp(\sigma_{3}\hat{v}_{\ell
r})$ (4.32)

where now $\hat{v}_{\ell r}(\hbar)$ is a function of $\hbar$ that is
asymptotic, in the Poincaré sense, to $v_{\ell r}(\hbar)$ in (4.31).

A similar proposition can be stated for the relation between a Stokes region
$\mathscr{D}$ and an external region $\mathscr{B}$. Indeed, if $\tau$ is the
turning point on the boundary of $\mathscr{D}$, we have

$[\psi_{+}^{(\infty)},\psi_{-}^{(\infty)}]=[\psi_{+}^{(\tau)},\psi_{-}^{(\tau)}]\exp\left(\sigma_{3}w_{\tau}\right)$
(4.33)

where $\tau$ is a turning point in the boundary of a Stokes region extending
to infinity, and

$w_{\tau}(\hbar):=R(z;\hbar)-\int_{\tau}^{z}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta$
(4.34)

where $R(z;\hbar)$ is defined in (• ‣ 4.1). Note that $w_{\tau}(\hbar)$ is a
constant with respect to $z$ since $\dfrac{\differential}{\differential
z}w_{\tau}(\hbar)=0$

In the case of the ODE (4.6) we will work under the assumption that the
turning points are simple and there are no “saddle trajectories” as specified
below.

######  Assumption 4.5.

The following assumptions shall prevail.

  * •

Simplicity. The roots of the potential are all simple. In the case of the
potential $Q(z;s,E)$, there are no repeated roots if and only if $(s,E)$
satisfy:

$Es^{4}-8E^{2}s^{2}+16E^{3}-s^{3}+36Es-27\neq 0.$ (4.35)

  * •

Genericity. There are no saddle trajectories i.e. there are no Stokes’ curves
connecting two turning points. Saddle trajectories can only occur if there is
$\gamma$ in the homology group of the Riemann surface $\overline{\Sigma}$ for
which

$\imaginary\oint_{\gamma}\sqrt{Q(z)}\differential z=0.$ (4.36)

The simplicity assumption means there are exactly three Stokes’ curves
emanating from each turning point. The genericity assumption implies that all
the Stokes’ curves must extend to $\infty$. With these assumptions we can
classify all the possible Stokes graphs.

In what follows we draw the trajectories on one sheet of the Riemann surface
$\Sigma$ relative to the choice of branch as in Def. 4.2.

######  Proposition 4.6 (Classification of generic Stokes graph).

Under the Assumptions 4.5 and with the choice of branch cuts in Definition
4.2, the Stokes graphs are in one-to-one correspondence with the
triangulations of the hexagon, so there are 14 such configurations. Three of
them are topologically distinct (as graphs), they are depicted in Fig. 5 and
named $E,D$ and $Z$. The remaining configurations can be obtained from the
graphs of types $E,D$ and $Z$ shown in Fig. 5 by a $\mathbb{Z}_{6}$ rotation
and by a reflection along the line from $\infty_{3}$ to $\infty_{0}$.

(a) Graph of type $D$

(b) Graph of type $E$

(c) Graph of type $Z$

Figure 5: Generic Stokes graph configurations for a quartic polynomial $Q(z)$.
Solid lines depict the Stokes curves emanating from the turning points
$\tau_{j}$, dashed lines denote the triangulation of the hexagon and green
curvy lines correspond to our choice of branch cuts. We remark that for
various $(s,E)\in\mathbb{C}^{2}$ the potential $Q(z;s,E)$ may have a different
different Stokes lines in $\overline{\mathbb{C}}$ but the underlying Stokes
graph will be topologically identical to one of the types depicted.

Proof. Assuming simplicity, from each turning point there are exactly three
Stokes curves emanating from it. These Stokes curves end either at infinity or
at another turning point. Assuming genericity the latter cannot happen,
therefore the Stokes curves determine an ideal triangulation of the Riemann
sphere with a small disk around infinity removed and with 6 marked points in
the boundary (corresponding to the asymptotic directions at infinity). These
triangulations correspond to triangulations of the hexagon, and there are 14
such triangulations.

Indeed, we can view these 6 marked points as determining a (topological)
hexagon. Furthermore, each of the turning points $\tau_{j}$ determines a
triangle inside the hexagon by connecting the asymptotic directions at
infinity with a Stokes curve originating from $\tau_{j}$. Up to rotations and
reflections of the hexagon there are 3 distinct such triangulations, named
$E$, $D$ and $Z$ as depicted in Fig. 5. The 14 possible configuration can be
obtained as follows. From configuration $D$ we obtain one other configuration
rotating by $2\pi/6$. From configuration $Z$ we obtain 6 distinct
configurations, 3 of them corresponding to a $\mathbb{Z}_{3}$ (i.e. $2\pi/3$)
rotation, and another 3 corresponding to a reflection followed by a
$\mathbb{Z}_{3}$ rotation. Finally from configuration $E_{+}$ we obtain 6
distinct configurations corresponding to a $\mathbb{Z}_{6}$ (i.e. $2\pi/6$)
rotation. Finally, $2+3+3+6=14$ as claimed. $\blacksquare$

Finally we can coordinate the Theorem 4.3 and Proposition 4.4 into a Riemann
Hilbert problem as follows.

######  Riemann-Hilbert Problem 4.7 (Quartic WKB jumps).

Fix a quartic potential $Q$ and suppose that it has one of the generic Stokes
graph configuration from Fig. 5. The Riemann-Hilbert problem for the vector
valued function $\Psi$ such that
$\Psi|_{\mathscr{D}}=(\Psi^{(\mathscr{D})}_{+},\Psi^{(\mathscr{D})}_{-})$ with
$\Psi|_{\mathscr{D}}$ defined in Theorem 4.3, consists of the following
oriented contours in $\overline{\mathbb{C}}$ with their associated jump
matrices.

  1. 1.

Square root branch cuts: To each square-root branch cut (coloured green in
Fig. 6 ), with the orientation indicated in Fig. 6, we associate the jump
matrix

$G:=\begin{bmatrix}0&i\\\ i&0\end{bmatrix}.$ (4.37)

  2. 2.

Fourth-root branch cuts: in configurations $E$ and $Z$ there is an extra jump
contour corresponding to the fact that

$\sqrt{S_{\rm{odd}}(z,\hbar)}\sim
h^{-1/2}Q(z)^{1/4}+\mathcal{O}(\hbar^{1/2}).$ (4.38)

To each fourth-root branch cut we associate the jump matrix $Y=-\mathbf{1}$.
The fourth-root branch cuts are coloured in yellow in Fig. 6. Note there is no
need to specify the orientation of these contours.

  3. 3.

Stokes curves: along each Stokes curve oriented towards $\oplus$ and away from
$\ominus$ we assign the jump matrices $B$ and $R$ (respectively)

$\displaystyle B:=$ $\displaystyle\begin{bmatrix}1&0\\\
-i&1\end{bmatrix},\qquad R:=\begin{bmatrix}1&i\\\ 0&1\end{bmatrix}.$ (4.39)

The Stokes curves are coloured blue or red, respectively, in Fig. 6.

  4. 4.

Inner ideal paths: along the inner ideal paths with the orientation in Def.
4.5 (and indicated in Fig. 6) we associate the following jump matrix
corresponding to the connection formula in Prop. 4.4:

$V_{jk}:=\exp\left(\sigma_{3}v_{jk}\right)=\begin{bmatrix}e^{v_{jk}}&0\\\
0&e^{-v_{jk}}\end{bmatrix},\quad
v_{jk}(\hbar)=\int_{\tau_{j}}^{\tau_{k}}S_{\rm{odd}}(z_{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+},\hbar)\differential
z$ (4.40)

where $j,k\in\\{0,1,2,3\\}$ and $\tau_{j},\tau_{k}$ are turning points. These
contours are denoted by a brown dashed line in Fig. 6.

  5. 5.

Outer ideal paths: along the outer ideal paths separating the external regions
from the Stokes regions we associate the jump matrix corresponding to the
connection formula (4.33) between turning points and infinity:

$W_{j}:=\exp(\sigma_{3}w_{j})=\begin{bmatrix}e^{w_{j}}&0\\\
0&e^{-w_{j}}\end{bmatrix},\quad
w_{j}(\hbar):=R(z;\hbar)-\int_{\tau_{j}}^{z}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta.$
(4.41)

where $j\in\\{0,1,2,3\\}$, $\tau_{j}$ is a turning point and $R(z;\hbar)$ is
the particular antiderivative in (• ‣ 4.1). These contours are denoted by a
black dashed line in Fig. 6.

######  Remark 4.8.

We observe that $w_{j}(\hbar)$ is a constant in $x$. It can be thought as the
regularized integral
$\int^{\tau_{j}}_{\infty}S_{\rm{odd}}(z,\hbar)\differential z.$

This construction gives us, up to rotation and reflections, three distinct WKB
Riemann-Hilbert problems as shown in Fig. 6.

(a) RHP for graphs of type $D$.

(b) RHP for graphs of type $E$.

(c) RHP for graphs of type $Z$.

Figure 6: Generic WKB Riemann-Hilbert problem corresponding to each Stokes’
graph configuration of the quartic polynomial potential $Q(z)$.

In order to simplify the upcoming computations we make some notational
definitions. Given turning points $\tau_{j}$ and $\tau_{k}$ connected by a
branch cut we denote

$\xi_{jk}:=\exp(2v_{jk})=\exp\left(2\int_{\tau_{j}}^{\tau_{k}}S_{\rm{odd}}(z_{+},\hbar)\differential
z\right),$ (4.42)

with the determination of $S_{\rm{odd}}$ given Def. 4.2, and the boundary
value $z_{+}$ is in accordance to Def. 4.3. If the two branchpoints are not
connected by a branch-cut we take the integration on the main sheet (i.e. with
the determination as above). We call the above parameters the Fock–Goncharov
parameters.

Additionally we will use the same notation to indicate the “exact”
Fock–Goncharov parameters, namely, the result of the Borel resummation. To
phrase it differently we will not distinguish in the notation the Borel
resummation by its asymptotic expansion in $\hbar$. Now we obtain the
following results about the Stokes matrices for each configuration shown in
Fig. 6.

######  Theorem 4.9.

The Stokes matrices

$\displaystyle\mathbb{S}_{j}=$ $\displaystyle\begin{bmatrix}1&0\\\
s_{j}&1\end{bmatrix},\;\;j=0,2,4,\hskip
25.60747pt\mathbb{S}_{j}=\begin{bmatrix}1&s_{j}\\\
0&1\end{bmatrix},\;\;\;j=1,5,\hskip
25.60747pt\mathbb{S}_{3}=\begin{bmatrix}1&s_{3}\\\ 0&1\end{bmatrix}{\rm
e}^{\frac{2i\pi}{\hbar}\sigma_{3}},$ (4.43)

in each of the WKB Riemann-Hilbert problems in Fig. 6 are expressed in terms
of the contour integrals $v_{jk}$, $w_{j}$ and $\xi_{jk}$ as defined in
(4.40), (4.41) and (4.42) as follows:

$\displaystyle\begin{array}[]{c|c|c}\hbox{Configuration
$D$}&\hbox{Configuration $E$}&\hbox{Configuration $Z$}\\\\[4.0pt]
\begin{array}[]{l}s_{0}=-ie^{2w_{1}}(\xi_{10}\xi_{30}+\xi_{10}+1),\\\\[4.0pt]
s_{1}=-ie^{-2w_{1}}\\\\[4.0pt]
s_{2}=-ie^{2w_{2}}(\xi_{20}\xi_{10}+\xi_{20}+1),\\\\[4.0pt]
s_{3}=-ie^{-2w_{2}}\\\\[4.0pt]
s_{4}=-ie^{2w_{3}}(\xi_{30}\xi_{20}+\xi_{30}+1),\\\\[4.0pt]
s_{5}=-ie^{-2w_{3}}.\end{array}&\begin{array}[]{l}s_{0}=-ie^{2w_{0}}(\xi_{01}\xi_{12}\xi_{32}+\xi_{01}\xi_{12}+\xi_{01}+1),\\\\[4.0pt]
s_{1}=-ie^{-2w_{0}}\\\\[4.0pt] s_{2}=-ie^{2w_{1}}(\xi_{01}+1),\\\\[4.0pt]
s_{3}=-ie^{-2w_{2}}(\xi_{12}+1)\\\\[4.0pt]
s_{4}=-ie^{2w_{3}}(\xi_{32}+1),\\\\[4.0pt]
s_{5}=-ie^{-2w_{3}}.\end{array}&\begin{array}[]{l}s_{0}=-ie^{2w_{0}}(\xi_{01}\xi_{12}+\xi_{01}+1),\\\\[4.0pt]
s_{1}=-ie^{-2w_{0}}\\\\[4.0pt] s_{2}=-ie^{2w_{1}}(\xi_{01}+1),\\\\[4.0pt]
s_{3}=-ie^{-2w_{3}}(\xi_{32}\xi_{12}+\xi_{32}+1)\\\\[4.0pt]
s_{4}=-ie^{2w_{3}}\\\\[4.0pt]
s_{5}=-ie^{-2w_{2}}(\xi_{32}+1).\end{array}\end{array}$ (4.49)

In the configuration $D$ the following identities hold

$\displaystyle{\rm e}^{2(w_{3}-w_{1})}=\xi_{10}\xi_{30},\ \ \ {\rm
e}^{2(w_{1}-w_{2})}=\xi_{10}\xi_{20},\ \ \ {\xi_{10}\xi_{20}\xi_{30}}={\rm
e}^{\frac{2i\pi}{\hbar}}.$ (4.50)

In configuration $E$ and $Z$ the following identity holds

$\xi_{01}\xi_{23}={\rm e}^{\frac{2i\pi}{\hbar}}.$ (4.51)

Proof. We compute the Stokes matrices associated to Configuration $D$ of the
Riemann-Hilbert problem in Fig. 6. Consider the vertex $\infty_{0}$ in
Configuration D, there are seven edges incident on it. The Stokes matrix
$\mathbb{S}_{0}$ is given by the clockwise (about the vertex $\infty_{0}$)
product of all the jump matrices corresponding to the edges incident to
$\infty_{0}$ as follows

$\mathbb{S}_{0}=W_{3}^{-1}BV_{30}BV_{10}BW_{1}$ (4.52)

Notice that the jump matrix $V_{30}$ follows from the fact that we have

$\exp\left(\int_{\tau_{3}}^{\tau_{0}}S_{\rm{odd}}(z_{+},\hbar)\differential
z\right)=\exp\left(\int_{\tau_{0}}^{\tau_{3}}S_{\rm{odd}}(z_{-},\hbar)\differential
z\right),$

where $z_{\pm}$ are the boundary values of $S_{\rm{odd}}(z,\hbar)$ on the left
and right side of the segment $[\tau_{0},\tau_{3}]$ oriented from $\tau_{3}$
to $\tau_{0}$. Therefore we have

$\mathbb{S}_{0}=\exp\left[\sigma_{3}(w_{1}-w_{3}+v_{10}+v_{30})\right]\
\begin{bmatrix}1&0\\\
-ie^{2w_{1}}(\xi_{10}\xi_{30}+\xi_{10}+1)&1\end{bmatrix}.$

We observe that

$\begin{split}w_{1}-w_{3}&=-\int_{\tau_{1}}^{z}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta+\int_{\tau_{3}}^{z}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta\\\
&={-\int_{\tau_{1}}^{\tau_{0}}S_{\rm{odd}}(\zeta_{+},\hbar)\differential\zeta-\int_{\tau_{0}}^{z}S_{\rm{odd}}(\zeta,\hbar)\differential\zeta+\int_{\tau_{3}}^{\tau_{0}}S_{\rm{odd}}(\zeta_{-},\hbar)\differential\zeta+\int_{\tau_{0}}^{z}S_{\rm{odd}}(\zeta,\hbar)}\differential\zeta\\\
&={-\int_{\tau_{3}}^{\tau_{0}}S_{\rm{odd}}(\zeta_{+},\hbar)\differential\zeta-\int_{\tau_{1}}^{\tau_{0}}S_{\rm{odd}}(\zeta_{+},\hbar)\differential\zeta}\end{split}$
(4.53)

and

$v_{10}+v_{30}=\int_{\tau_{1}}^{\tau_{0}}S_{\rm{odd}}(z_{+},\hbar)\differential
z+\int_{\tau_{3}}^{\tau_{0}}S_{\rm{odd}}(z_{+},\hbar)\differential z$ (4.54)

where the integral from $\tau_{1}$ to $\tau_{3}$ is with the determination in
Def.4.2. Therefore we obtain

$w_{1}-w_{3}+v_{10}+v_{30}=0,$ (4.55)

so that

$\exp\left(\sigma_{3}[w_{1}-w_{3}+v_{10}+v_{30}]\right)=\mathbb{I},$ (4.56)

which is equivalent to the first identity in (4.50). Similar calculations
compute the remaining Stokes matrices for each possible configuration. The
extra ${\rm e}^{2i\pi\hbar^{-1}\sigma_{3}}$ in the form of $\mathbb{S}_{3}$ is
due to our choice of branch-cut for the logarithm in (• ‣ 4.1). The fact that
the product of all Stokes matrices is trivial follows by construction.

The second equation in (4.50) for configuration $D$ follows from the
definitions of $\xi_{ij}$ in (4.42) and $w_{j}$ in (4.41). The last equation
follows from the fact that $\gamma_{10}+\gamma_{20}+\gamma_{30}$ is homologous
to $\Gamma_{\infty}$, the contour at infinity with winding number equal to
one. The residue at infinity of the (formal series in $\hbar$)
$S_{\rm{odd}}(z,\hbar)\differential z$ is $2i\pi/\hbar$ coming from the
leading term of the series only (it is shown in [KT05], which can be also
verified by induction, that all the higher terms in the $\hbar$–series
expansion of $S_{\rm{odd}}$ vanish to order $\mathcal{O}(z^{-2})$ and hence do
not contribute to the residue). Regarding the relation (4.51) since the Stokes
graph of $E$ and $Z$ have the same topology of branch points, we consider both
cases simultaneously. One verifies that the contour $\gamma_{01}+\gamma_{23}$
is homologous to $\Gamma_{\infty}$. Thus in a similar way as the previous case

$\left(\oint_{\gamma_{01}}+\oint_{\gamma_{23}}\right)S_{\rm{odd}}(z,\hbar){\mathrm{d}}z=\text{Res}_{z=\infty}\hbar^{-1}\sqrt{Q(z,s,e)}dz=\dfrac{2\pi
i}{\hbar}.$

$\blacksquare$

##### Gauge arbitrariness.

In terms of the Stokes phenomenon, we must point out that the Stokes
parameters $s_{0},...,s_{5}$ are not intrinsically defined because we can
conjugate the fundamental matrix by an arbitrary diagonal matrix. This freedom
translates to the following scaling equivalence for the Stokes parameters

$\displaystyle s_{2j+1}\mapsto\lambda s_{2j+1},\ \ \
s_{2j}\mapsto\lambda^{-1}s_{2j},\ \ \ \ \lambda\in\mathbb{C}^{\times}.$ (4.57)

Using this freedom we can rewrite (in all cases) the Stokes parameters in such
a way that only the Vorös’ symbols $v_{jk}$ enter. For example, in the $D$
configuration (which will be the only one that is relevant later on), we have

$\displaystyle s_{0}=$ $\displaystyle
i(\xi_{10}\xi_{30}+\xi_{10}+1),\qquad\qquad s_{1}=i$ $\displaystyle s_{2}=$
$\displaystyle
ie^{2w_{2}-2w_{1}}(\xi_{20}\xi_{10}+\xi_{20}+1)=i\frac{1}{\xi_{20}\xi_{10}}(\xi_{20}\xi_{10}+\xi_{20}+1),$
$\displaystyle s_{3}=$ $\displaystyle ie^{-2w_{2}+2w_{1}}=i\xi_{10}\xi_{20}$
$\displaystyle s_{4}=$ $\displaystyle
ie^{2w_{3}-2w_{1}}(\xi_{30}\xi_{20}+\xi_{30}+1)=i\xi_{30}\xi_{10}(\xi_{30}\xi_{20}+\xi_{30}+1)$
$\displaystyle s_{5}=$ $\displaystyle
ie^{-2w_{3}+2w_{1}}=i\frac{1}{\xi_{10}\xi_{30}}.$ (4.58)

##  5 Quantization conditions

In this section we derive “exact” quantisation conditions for both the zeroes
of Vorob’ev–Yablonskii polynomials (the Jimbo-Miwa case) and for the points of
the ES spectrum corresponding to repeated eigenvalues (the Shapiro-Tater
case). To leading order these quantisation conditions yield a system that
describes both sets of points in terms of contour integrals of
$\sqrt{Q(z;s,E)}$.

In the Jimbo-Miwa case this is achieved by matching the Stokes’ phenomenon of
the Jimbo-Miwa Lax pair representation of PII with the Stokes’ phenomenon of
the quartic WKB Riemann-Hilbert problem 4.9. This result is contained in
Theorem 5.2.

In the Shapiro-Tater case it is not enough matching Stokes phenomenon in
Theorem 3.3 to obtain quantisation conditions, this can seen from Theorem 5.3.
For this reason we additionally impose the condition (3.57) which yields the
correct quantisation equations (5.23) from an application of Theorem 5.4.

###  5.1 The Jimbo-Miwa case

The parameters $s,E$ in the potential $Q(z;s,E)$ i.e. the parameters
$t,\Lambda$ corresponding to the pole with residue $-1$ of the rational
solution $u_{n}$ and “eigenvalue” of the (2.17), are determined by the
implicit requirement that the Stokes phenomenon for the ODE matches the one
indicated below.

Indeed it was shown in [BM14] that rational solutions of the Painlevé II
equation correspond to a particular Stokes phenomenon as shown in Fig. 7.

We recall that the map to the Stokes’ data for general solution of the
Painlevé II transcendent was obtained originally in [IN86], see also [FIKN06].

Figure 7: Stokes data for the Lax pair (2.1) corresponding to rational
solutions of Painlevé II

######  Theorem 5.1 ([BM14]).

The matrices $\mathbb{S}_{0},\dots,\mathbb{S}_{5}$ in Fig. 7 form the
monodromy data of the Jimbo-Miwa Lax pair (2.1) corresponding to the rational
solutions of PII.

Consider a rational solution $u_{n}(t)$ (1.8) of the PII equation with
$\alpha=n\in\mathbb{Z}$. Let $a$ be a pole of residue $-1$, namely a zero of
the Vorob’ev-Yablonskii polynomial $Y_{n}(t)$. Let $b$ be the coefficient in
the expansion as in (2.8), and recall $\Lambda=\frac{7a^{2}}{36}+10b$ as in
(2.17). Then, according to Prop. 2.3, the Stokes phenomenon of the ODE (2.18)
must be as in Fig. 7. Viceversa [Mas10b], if the anharmonic potential (2.18)
exhibit the Stokes phenomenon as in Fig. 7, then the pair of values $a,b$
characterize (uniquely) the (rational, in this case) solution $u(t)$ of the
Painlevé II equation with a pole at $t=a$.

Thus, to find the positions of a pole, we can find for which values of $a,b$
in (2.18) (with $\alpha=n$) the Stokes phenomenon matches Fig. 7. Of course
the map that associates to the parameters $a,b$ in (2.18) the Stokes data is
highly transcendental.

It is the nature of our problem, however, that we are interested in the
behaviour when $n$ is large and also in the re-scaled plane. Thus we can apply
our (exact) WKB analysis to obtain asymptotic information on the Stokes
parameters and set up an implicit equation for $a,b$, or rather their rescaled
counterparts $s,E$ as in (4.5) with the identifications (4.2).

In order to apply the exact WKB method, we set our large parameter to be
$\hbar^{-1}=(n+1/2)$. By the general theory in Section 4, we construct the WKB
solutions associated to (4.6) as formal power series in $\hbar$:

$\psi_{\pm}^{(\tau)}(z,\hbar)=\frac{1}{\sqrt{S_{\rm{odd}}(z;\hbar)}}\exp\left(\pm\int_{\tau}^{z}S_{\rm{odd}}(u;\hbar)\differential
u\right)$ (5.1)

normalized near a turning point $\tau$.

######  Theorem 5.2 (Vorob’ev-Yablonskii quantisation).

Suppose that $(a,b)\in\mathbb{C}^{2}$ determine a rational solution of PII
with $\alpha=n$ and Laurent expansion (2.1). Let $s=\hbar^{2/3}a$ and
$E=7s^{2}/36+\hbar^{4/3}b$ with $\hbar^{-1}=n+1/2$ as in the scaling (4.2).
Then the Stokes’ graph for the potential $Q(z;s,E)$ must be of type $D$ and
the corresponding Fock–Goncharov parameters (4.42) must satisfy

$\xi_{10}=\xi_{20}=\xi_{30}=-1.$ (5.2)

In terms of the Vorös symbols we obtain the following leading order estimates:

$\displaystyle\pi
i(2k_{1}+1)=2\int_{\tau_{1}}^{\tau_{0}}S_{\rm{odd}}(z_{+},\hbar)\differential
z\simeq(2n+1)\int_{\tau_{1}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}\differential
z+\mathcal{O}(\hbar)$ (5.3) $\displaystyle\pi
i(2k_{2}+1)=2\int_{\tau_{2}}^{\tau_{0}}S_{\rm{odd}}(z_{+},\hbar)\differential
z\simeq(2n+1)\int_{\tau_{2}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}\differential
z+\mathcal{O}(\hbar)$ (5.4) $\displaystyle\pi
i(2k_{3}+1)=2\int_{\tau_{3}}^{\tau_{0}}S_{\rm{odd}}(z_{+},\hbar)\differential
z\simeq(2n+1)\int_{\tau_{3}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}\differential
z+\mathcal{O}(\hbar)$ (5.5)

where $k_{1},k_{2},k_{3}$ are integers and $\hbar=(n+1/2)^{-1}$

Proof. We take the Stokes data for the Jimbo-Miwa Lax pair corresponding to
the rational solutions of PII in Fig. 7 and we equate it to the Stokes
matrices from the WKB Riemann-Hilbert problem in each configuration of Theorem
4.9 i.e.:

$\mathbb{S}_{2j}=\begin{bmatrix}1&0\\\
i&1\end{bmatrix},\quad\mathbb{S}_{2j+1}=(-1)^{\delta_{j1}}\begin{bmatrix}1&i\\\
0&1\end{bmatrix}\quad(j=0,1,2)$ (5.6)

From each configuration we obtain a system of 6 equations involving the
exponentials of the periods $v_{j}$, for example in Configuration $D$ we have
from (4.58) the conditions:

$\displaystyle s_{0}=$ $\displaystyle i=i(\xi_{10}\xi_{30}+\xi_{10}+1),\qquad
s_{1}=i$ $\displaystyle s_{2}=$ $\displaystyle
i=i\frac{1}{\xi_{20}\xi_{10}}(\xi_{20}\xi_{10}+\xi_{20}+1),\quad
s_{3}=i=i\xi_{10}\xi_{20}$ $\displaystyle s_{4}=$ $\displaystyle
i=i\xi_{30}\xi_{10}(\xi_{30}\xi_{20}+\xi_{30}+1),\quad
s_{5}=i=i\frac{1}{\xi_{10}\xi_{30}}\,.$

The only solution of the above system is $\xi_{20}=\xi_{10}=\xi_{30}=-1$ which
is consistent with the fact that the product is $-1$. Direct inspection of the
formulas in Theorem 4.9 shows that it is impossible to satisfy the constraints
(5.6) in all of the other configurations. $\blacksquare$

###  5.2 The Shapiro-Tater case

######  Theorem 5.3.

Suppose that the quasi-polynomial $y=p_{n}(x)e^{\theta(x;t)}$ is a solution to
the boundary problem (3.1)-(3.2) with $J=n+1$. Then the “exact” Fock–Goncharov
parameters $\xi_{jk}$ in (4.42) satisfy one of the following systems,
depending on the indicated Stokes graph configuration.

  * •

Configuration $D$ (Fig. 5) or its $\mathbb{Z}_{3}$ rotations:

$\begin{cases}\xi_{10}\xi_{30}+\xi_{10}+1=0\\\
\xi_{20}\xi_{10}+\xi_{20}+1=0\\\ \xi_{30}\xi_{20}+\xi_{30}+1=0\end{cases}$
(5.7)

This systems cuts an affine rational curve in $\mathbb{C}^{3}$ given by:

$\xi_{10}=\rho,\quad\xi_{20}=-\frac{1}{\rho+1},\quad\xi_{30}=-\frac{\rho+1}{\rho}\,.$
(5.8)

  * •

Configuration $E$ (Fig. 5) or its $\mathbb{Z}_{3}$ rotations:

$\begin{cases}\xi_{01}\xi_{12}(\xi_{23}+1)+\xi_{01}+1=0\\\ \xi_{01}+1=0\\\
\xi_{23}+1=0\,.\end{cases}$ (5.9)

Furthermore all other configurations cannot occur.

Proof. In Proposition 3.2 it was shown that if $J=n+1$ and $(t,\Lambda)$
belong to the ES spectrum (i.e. there is a solution of the boundary problem
(3.1) and (3.2)) then the Stokes parameters $s_{0},s_{2},s_{4}$ all vanish
simultaneously.

In Theorem 4.9 and using the gauge transformation (4.58) we have expressed the
Stokes parameters in terms of the exact Fock-Goncharov parameters $\xi_{jk}$:
thus we have to see which configurations are compatible with the three
equations $0=s_{0}=s_{2}=s_{4}$. Note that equations (5.7) and (5.9)
corresponds exactly to the condition $0=s_{0}=s_{2}=s_{4}$ in configuration
$D$ and $E$ as in Fig.(6) or its $Z_{3}$ reflections. Note that
$\hbar=1/(n+1)$ and hence ${\rm e}^{\frac{2i\pi}{\hbar}}=1$ in the Theorem.
Direct inspection and simple algebra then allows us to rule out configuration
$Z$ as well as all the other configurations obtained from it by a
$\mathbb{Z}_{3}$ rotation or reflection, and those obtained from $D$ or $E$ by
a reflection around the imaginary axis followed by any $\mathbb{Z}_{3}$
rotation. $\blacksquare$

####  5.2.1 Repeated eigenvalue condition

Theorem 5.3 establishes the conditions for the Vöros symbols to yield a point
in the ES spectrum; together with Theorem 5.3 the conditions are equivalent to
the statement that the Stokes’ graph is either of $D$ or $E$ type and the
Fock-Goncharov parameters $\xi_{jk}$ satisfy the corresponding conditions
specified in Theorem 5.3.

In addition we must now impose the condition that the eigenvalue is a repeated
one: as proved in Theorem 3.9 this requires that all the integrals of
$p_{n}^{2}{\rm e}^{2\theta}$ between $\infty_{2k+1}$ vanish. It was also
explained in the theorem that it suffices to impose one of the two vanishing
conditions and the other one will follow. Thus the strategy now is to compute

$\displaystyle I_{13}=\int_{\infty_{1}}^{\infty_{3}}p_{n}(x)^{2}{\rm
e}^{2\theta(x;t)}{\mathrm{d}}x\,,$ (5.10)

using the asymptotic expansion in terms of formal WKB solutions obtained so
far. In order to simplify the upcoming computations we will label the regions
in the WKB RH problem for configuration $D$ and $E$ according to Fig. 8. This
will help us distinguish between the entire solutions of the differential
equation (4.6) that are asymptotics to WKB solutions $\psi_{\pm}^{\tau_{j}}$
in Def. 4.1 in different regions in accordance with Thm. 4.3. Namely

$\psi_{\pm}^{\mathcal{B}}(z)\simeq\psi_{\pm}^{\tau_{j}}(z;\hbar)=\frac{1}{\sqrt{S_{\rm{odd}}(z,\hbar)}}\exp\left(\pm\int^{z}_{\tau}S_{\rm{odd}}(u,\hbar)\differential
u\right),\quad u\in\mathcal{B},\;\hbar\to 0$

where $\mathcal{B}$ is one of the labelled regions in Fig. 8 ad $\tau_{j}$ is
one of the turning points in its boundary. We observe that to compute $I_{13}$
we can equivalently compute the integral of $\Psi_{+}^{(A_{1})}$ because this
function is proportional to $p_{n}(x)^{2}{\rm e}^{2\theta(x;t)}$ since they
are both recessive in the direction $\infty_{1}$. Thus, the main aim of this
section is to prove the following Theorem.

######  Theorem 5.4.

Let $(t,\Lambda)$ belong to the ES spectrum. Let $s=t\hbar^{\frac{2}{3}}$ and
$E=\Lambda\hbar^{\frac{4}{3}}$, with $\hbar=(n+1)^{-1}$. Then

  1. 1.

In configuration $D$ as in Fig. (8) we have:

$\int_{\infty_{1}}^{\infty_{3}}\Psi_{+}^{(A_{1})}(z;\hbar)^{2}{\mathrm{d}}z\simeq
2i\hbar\left({\rm
e}^{2v_{12}}\int_{\tau_{2}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}-\int_{\tau_{1}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}\right)+\mathcal{O}(\hbar^{2})\,$
(5.11)

where $v_{12}$ is defined in (4.31). The rescaled parameters $(s,E)$ of the ES
spectrum are asymptotic to a repeated eigenvalue provided that

$\displaystyle\exp\left[{\frac{2}{\hbar}\int_{\tau_{1}}^{\tau_{2}}\sqrt{Q(z;s,E)}{\mathrm{d}}z}\right]=\boldsymbol{\tau}(s,E)+\mathcal{O}(\hbar)$
(5.12)

where

$\displaystyle\boldsymbol{\tau}(s,E)=\frac{\displaystyle\int_{\tau_{1}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}}{\displaystyle\int_{\tau_{2}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}},\quad\imaginary(\boldsymbol{\tau}(s,E))>0.$
(5.13)

  2. 2.

In configuration $E$ as in Fig. (6) we have

$\displaystyle\int_{\infty_{1}}^{\infty_{3}}\Psi_{+}^{(A_{1})}(z;\hbar)^{2}{\mathrm{d}}z\simeq
2i\int_{\tau_{2}}^{\tau_{0}}\frac{\hbar{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}+\mathcal{O}(\hbar^{2}).$
(5.14)

Then the rescaled parameters $(s,E)$ of the ES spectrum cannot be a double
eigenvalue for large $n$.

The proof is relatively straightforward but significantly delicate and
technical and it is postponed to the Appendix A.1. We can take the naive
approach in which we replace $\Psi_{+}^{(A_{1})}$ with the appropriate
combinations of the WKB formal solutions along the pieces of the contour of
integration that traverse each Stokes region. In fact this approach yields the
correct result, but over-estimates the error term. It is however appealing for
its simplicity to give here the heuristic derivation and defer technicalities
to later. Consider the case of configuration $D$: we split the integration at
the turning points $\tau_{1},\tau_{0},\tau_{2}$. The unbounded integrals are
then along steepest descent paths for the integrand and can be neglected. In
the integrations along $[\tau_{1},\tau_{0}]$ and $[\tau_{0},\tau_{2}]$ we
observe that, as a consequence of the Riemann-Hilbert problem 4.7 we can
express $\Psi^{(A_{1})}_{+}$ as suitable linear combinations of
$\psi_{\pm}^{(\tau_{j})}$; indeed

$\displaystyle\Psi_{+}^{(A_{1})}(z;\hbar)$
$\displaystyle\simeq\left\\{\begin{array}[]{ll}\psi_{+}^{(\tau_{1})}(z;\hbar)+i\psi_{-}^{(\tau_{1})}(z;\hbar)&z\in(B_{1}^{r})\cup(B_{1}^{\ell})\\\\[4.0pt]
{\rm
e}^{v_{12}}\left(\psi_{+}^{(\tau_{2})}(z;\hbar)-i\psi_{-}^{(\tau_{2})}(z;\hbar)\right)&z\in(C_{1}^{\ell})\cup(C_{1}^{r})\end{array}\right.$
(5.17)

where we have used already the fact that we are on the ES spectrum so that
Theorem 5.3 applies and the Stokes matrices $\mathbb{S}_{j},\,j\in\\{0,2,4\\}$
are trivial. Computing the square of $\Psi_{+}^{(A_{1})}$, we have that the
cross–products yield non-oscillatory functions that contribute to the leading
order while the squares of the “pure” WKB solutions give oscillatory integrals
which can be neglected to leading order.

Thus one is lead to the rough estimate

$\displaystyle\int_{\infty_{1}}^{\infty_{3}}\left(\Psi_{+}^{(A_{1})}\right)^{2}{\mathrm{d}}z\simeq
2i\int_{\tau_{1}}^{\tau_{0}}\psi_{+}^{(\tau_{1})}\psi_{-}^{(\tau_{1})}{\mathrm{d}}z-2i{\rm
e}^{2v_{12}}\int_{\tau_{0}}^{\tau_{2}}\psi_{+}^{(\tau_{2})}\psi_{-}^{(\tau_{2})}{\mathrm{d}}z\simeq$
$\displaystyle\simeq
2i\int_{\tau_{1}}^{\tau_{0}}\frac{\hbar{\mathrm{d}}z}{\sqrt{Q(z_{-};s,E)}}-2i{\rm
e}^{2v_{12}}\int_{\tau_{0}}^{\tau_{2}}\frac{\hbar{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}$
(5.18)

where the boundary values of $\sqrt{Q}$ are due to our choice of orientations
for the branch-cuts in Fig. 8 (i.e. not according to Def. 4.3). Rearranging
the endpoints and boundary values yields (5.11).

The reason why the above reasoning is defective is that it replaces the formal
WKB expansions also in the neighbourhoods of the turning points, where the
formal WKB solutions have a singularity. One may still make sense of the
resulting integrals because they involve a singularity of type
$(z-\tau)^{-\frac{1}{2}}$ which is integrable. However, approaching the
integrals in this way and using a (formal) application of the Laplace method
would suggest that the subleading order is $\mathcal{O}(\hbar^{\frac{4}{3}})$.
The careful analysis, instead, of the contribution near the turning points
reveals that the subleading correction is of order is
$\mathcal{O}(\hbar^{2})$.

Unfortunately we could not find a reference in the vast literature on exact
WKB analysis that helps us on this issue. For this reason we have postponed
this part of the proof to the appendix.

(a) Configuration $D$

(b) Configuration $E$

Figure 8: Labelled regions in the WKB Riemann-Hilbert problem.

###  5.3 Comparison of quantization conditions

##### The ST case.

In view of Theorem 5.4 and in particular (5.12), we can now express, to within
the leading order, the quantization conditions that characterize those points
$(s,E)$ in the ES spectrum with a double eigenvalue. Indeed from (5.8) it
follows that

$\displaystyle-\frac{\xi_{10}}{1+\xi_{10}}=\frac{1}{\boldsymbol{\tau}}\qquad\Rightarrow\qquad\left\\{\begin{array}[]{l}\displaystyle\xi_{10}=\frac{-1}{1+\boldsymbol{\tau}}\\\\[8.0pt]
\displaystyle\xi_{20}=-1-\frac{1}{\boldsymbol{\tau}}\\\\[8.0pt]
\displaystyle\xi_{30}={\boldsymbol{\tau}}\end{array}\right.$ (5.22)

In terms of the Voros’ symbol we then have, to leading order

$\displaystyle
2(n+1)\int_{\tau_{1}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=\ln\left(\frac{-1}{1+\boldsymbol{\tau}(s,E)}\right)-2i\pi(m_{1}+1)$
$\displaystyle
2(n+1)\int_{\tau_{2}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=\ln\left(-1-\frac{1}{\boldsymbol{\tau}(s,E)}\right)-2i\pi(m_{2}+1)$
$\displaystyle
2(n+1)\int_{\tau_{3}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=\ln\left({\boldsymbol{\tau}(s,E)}\right)-2i\pi(m_{3}+1)$
(5.23)

where the three integers satisfy $m_{1}+m_{2}+m_{3}=n-1$ due to the fact that
the sum of the three integrals on the left is $-2i\pi(n+1)$ while the sum of
the three logarithms is $2i\pi$ (principal determination) due to the
definition of $\boldsymbol{\tau}(s,E)$ as in (5.13).

##### The JM case.

On the other hand, the quantization conditions for the Vorob’ev-Yablonskii
zeroes, to the same order of approximation, read

$\displaystyle\left(2n+1\right)\int_{\tau_{j}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=-{i\pi}-2i\pi
k_{j}$ $\displaystyle k_{1}+k_{2}+k_{3}=n-1.$ (5.24)

Both conditions (5.23), (5.24) involve a triple of positive integers adding to
$n-1$ but they differ notably in the multiplicative factor $2(n+1)$ vs.
$(2n+1)$ on the left side, and on the values on the right side. We now analyze
the two lattices to explain their similarity which is apparent from the
numerical experiments.

###  5.4 The elliptic region

Figure 9: The phase space in the $s$–plane; depicted are the curves defined by
(5.33). The inserted vignettes depict the critical graph $\mathfrak{C}$
corresponding to a sample value of $s$ in their approximate location (which is
at the centre of each vignette). The types of topologies in each region is
marked by the corresponding letter in the figure. The central triangular
region corresponds to the $D$ topology; the three “corridors”, adjacent to the
sides, in the directions $\arg(s)=0,\tfrac{\pm 2\pi}{3}$ correspond to the
three $\mathbb{Z}_{3}$ rotation of a $E$, while in the six remaining regions
we find the various $\mathbb{Z}_{2}\times\mathbb{Z}_{3}$ versions of the $Z$
topologies. Also indicated some of the topologies with coalesced turning
points, which are found when $s$ is on any of the boundaries between regions.

By the term “elliptic region” we refer, with a nod to the terminology in
[BM14, BM15, BB15], to the triangular region containing the (scaled) zeros.
The boundary of this region is determined in loc. cit. for the zeroes of the
Vorobev–Yablonski polynomials. We want to show that (at least asymptoticallyl)
also the lattice (described in detail in the next section) of rescaled zeroes
of $D_{n}(t)$ are contained in the same region.555We thank the anonymous
referee for pointing out this aspect of the conjecture. In order to argue that
the (rescaled) zeros of the discriminant $D_{n}(t)$ and of the Vorob’ev-
Yablonski polynomials $Y_{n}(t)$ fill the same region of the $s$–plane, we
should show that the $D$–configuration occurs only within the elliptic region.
For the Vorob’ev-Yablonski polynomials this issue has been addressed in the
above literature. For the Shapiro-Tater problem we give only a semi-rigorous
argument in this section (with numerical support).

We have established that only Stokes’ graphs of topology $D$ can give rise to
a solution of the repeated eigenvalue condition. In fact the equations (5.23)
imply that, to leading order as $n\to\infty$, the three periods
$\int_{\tau_{j}}^{\tau_{0}}\sqrt{Q(z)}\differential z$ are all purely
imaginary:

$\displaystyle\int_{\tau_{j}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=-i\pi\frac{m_{j}+1}{n+1}+\mathcal{O}(n^{-1}),\quad
j=1,2,3.$ (5.25)

We take these equations to leading order;

$\displaystyle\int_{\tau_{j}}^{\tau_{0}}\sqrt{Q(z_{+};s,E)}{\mathrm{d}}z=-i\pi\mu_{j},\
\ \ \mu_{1}+\mu_{2}+\mu_{3}=1,\ \ \mu_{j}\in\mathbb{R}.$ (5.26)

This means that the differential $\sqrt{Q}\differential z$ is (to leading
order) an imaginary normalized differential and the algebraic elliptic curve
$w^{2}=Q(z;s,E)$ is termed, in the literature [Ber11], a Boutroux curve, with
the condition that all the periods (5.26) are imaginary being the Boutroux
condition. Under this condition we know that the critical graph,
$\mathfrak{C}$, is connected. This is, by definition, the graph determined by
the “vertical critical trajectories666The ”vertical” here refers to the fact
that in the plane of the locally defined conformal map
$\xi(z):=\int^{z}\sqrt{Q}\differential w$ these are vertical straight lines.”
[Str84], namely, the maximal integral curves of the direction field

$\displaystyle\real\sqrt{Q(z;s,E)}\differential z=0$ (5.27)

that issue from each of the turning points. Observe that the integral curves
of (5.27) form a foliation that is orthogonal to the Stokes’ graph (4.27), in
the sense of the Riemannian metric
$\differential\sigma^{2}:=|Q(z;s,E)||\differential z|^{2}$ (with conical
singularities at the turning points). Briefly this can be seen by considering
the real–valued function
$\Phi(z):=\real\left(\int_{\tau_{0}}^{z}w\differential z\right)$, which is
easily shown to be

  1. –

harmonic on the Riemann surface $\Sigma$ (4.20) minus the two points above
$z=\infty$;

  2. –

odd under the holomorphic involution $(w,z)\to(-w,z)$ (i.e. under the flip of
sign of the square-root).

The connectedness of $\Phi^{-1}(\\{0\\})$ on $\Sigma$ (and hence also of its
$z$–projection) follows from the fact that $\Phi$ is a Morse function whose
critical points (the turning points) belong all to the zero level-set (due to
the Boutroux conditions). These two properties imply that the zero level set
$\Phi^{-1}(\\{0\\})$ is a well–defined subset of the complex $z$–plane (i.e.
of the projection from $\Sigma$ to the $z$–plane). It is then not hard,
following similar arguments to those used in analyzing the topology of the
Stokes’ graph, that $\mathfrak{C}$ can only have three different topologies
corresponding to the three configurations $D,E,Z$ (and “dual” to the Stokes’
graph). The graph $\mathfrak{C}$ is however also a metric space with the
distance induced by the Riemannian metric $\differential\sigma^{2}$; in fact
the geodesic length between two turning points (along the graph) is precisely
the (absolute value) of the purely imaginary integral of
$\int_{\tau_{i}}^{\tau_{j}}\sqrt{Q}\differential z$ along any path homotopic
to the edge of $\mathfrak{C}$ that connects the two turning points. In other
words the numbers $\pi|\mu_{j}|$ are distances, and the vanishing of any of
them indicates that the two turning points have coalesced; therefore this
happens when the topology of the critical graph undergoes a change. In
particular this can only happen if at least two turning points coincide,
namely on the critical locus where the discriminant of $Q$ (4.35) vanishes.
Summarizing the above discussion, the boundary separating between different
configurations is given by system of equations (of mixed
algebraic/transcendental nature)

$\displaystyle\left\\{\begin{array}[]{cc}Es^{4}-8E^{2}s^{2}+16E^{3}-s^{3}+36Es-27=0\\\
\real\int_{\tau_{0}}^{\mu}\sqrt{Q}\differential z=0\end{array}\right.$ (5.30)

where $\mu$ denotes the multiple root of $Q(z;s,E)$ and $\tau_{0}$ a remaining
simple root.

It is possible to simplify the system (5.30): the first equation can be
rationally parametrized

$\displaystyle E=\frac{8a^{3}+1}{16a^{4}},\ \ s=\frac{4a^{3}-1}{2a^{2}};\ \ \
Q(z;s({a}),E({a}))=\big{(}z+\tfrac{1}{2a}\big{)}^{2}\left(z^{2}-\frac{z}{a}+\frac{1+8a^{3}}{4a^{2}}\right)$
(5.31)

A further degeneration is when three turning points coalesce, and from the
above we see that this happens if and only if $2a^{3}+1=0$. The three
corresponding values of $s$ are

$\displaystyle s_{0}=-\frac{3}{2^{\frac{1}{3}}}\simeq-2.381101,\ \ \
s_{k}=s_{0}{\rm e}^{\frac{2i\pi k}{3}},\ \ k=1,2.$ (5.32)

and the boundary of the regions are determined by the image in the $s$–plane
(under the map $s=\frac{4a^{3}-1}{2a^{2}}$) of the curves defined implicitly
in the $a$–plane by

$\displaystyle\real\int_{-\frac{1}{2a}}^{\tau_{0}}\sqrt{Q(z;s(a),E(a))}\differential
z=\real\left[\frac{4a^{3}-1}{6a^{3}}\sqrt{2a^{3}+1}+\ln\left(\frac{1-\sqrt{2a^{3}+1}}{a^{\frac{3}{2}}\sqrt{2}}\right)\right]=0.$
(5.33)

We observe that the above condition does not depend on the choice of
$\tau_{0}$, namely one of the two remaining roots of the quartic polynomial,
$Q(z;s(a)),E(a))$. We remark that these are the same equations that appeared
in [BB15] (Definition 1.3 ibidem) defining the boundary of the elliptic
region. The curves defined by (5.33) partition the $s$–plane into 10 regions
(see Fig. 9); in each region the topology of the critical graph of the
Boutroux curve is the same, and the boundary between regions corresponds to
two turning points coalescing. The “elliptic region” is the triangular shaped
region containing the origin and it is the only region where the topology of
the critical graph is of type $D$. This was shown in [BM14] and we will not
repeat the arguments here because it is a rather lengthy affair. We rather
offer a numerical investigation of the “phase space” summarized in Fig. 9. A
complete proof can be obtained following the ideas in [BT16, Ber11] of
performing a “continuation in parameter space” and carefully tracking the
change in topology of the critical graph $\mathfrak{C}$.

The important issue for us is that since only the elliptic region $D$ is
compatible with the solution of the Shapiro-Tater problem, we deduce that in
any compact sets of the complement of the elliptic region there are (at least
for $n$ large enough) none of the rescaled zeros of $D_{n}(t)$.

This also addresses another aspect of the Shapiro-Tater conjecture.

###  5.5 Analysis of the two lattices

Figure 10: Canonical basis of cycles for the homology of the elliptic Riemann
surface $\overline{\Sigma}$ with Stokes graph configuration of type $D$. The
points $\tau_{j}$ are the branch points of $\sqrt{Q(z;s,E)}$ and the lines
connecting them are the branch cuts.

Both lattices involve implicit equations for the parameters $(s,E)$ via the
periods of the differential $\sqrt{Q(z;s,E)}{\mathrm{d}}z$. We introduce a
canonical basis of cycles $\mathcal{A},\mathcal{B}$ of the elliptic Riemann
surface $\overline{\Sigma}$ as in the Fig.10 so that

$\displaystyle\boldsymbol{\tau}(s,E)=\frac{\displaystyle\int_{\tau_{1}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}}{\displaystyle\int_{\tau_{2}}^{\tau_{0}}\frac{{\mathrm{d}}z}{\sqrt{Q(z_{+};s,E)}}}=\frac{\displaystyle\int_{\mathcal{B}}\frac{{\mathrm{d}}z}{\sqrt{Q(z;s,E)}}}{\displaystyle\int_{\mathcal{A}}\frac{{\mathrm{d}}z}{\sqrt{Q(z;s,E)}}}.$
(5.34)

Notice that the quantities $v_{j0}$ defined in (4.40) can be written in the
form

$\displaystyle
v_{10}=\frac{1}{2}\int_{\mathcal{B}}\sqrt{Q(z;s,E)}{\mathrm{d}}z,\;\;v_{20}=\frac{1}{2}\int_{\mathcal{A}}\sqrt{Q(z;s,E)}{\mathrm{d}}z,\;\;v_{30}=-\frac{1}{2}\int_{\mathcal{A}+\mathcal{B}}\sqrt{Q(z;s,E)}{\mathrm{d}}z$
(5.35)

The Jacobian determinant is a constant as we prove in the next lemma.

######  Lemma 5.5.

Let $\mathcal{A}$ and $\mathcal{B}$ be the canonical homology basis as in Fig.
10 and consider the periods

$I_{\mathcal{A}}=\oint_{\mathcal{A}}\sqrt{Q(z;s,E)}{\mathrm{d}}z,\quad
I_{\mathcal{B}}=\oint_{\mathcal{B}}\sqrt{Q(z;s,E)}{\mathrm{d}}z.$ (5.36)

Then