Exactly Solvable Anharmonic Oscillator, Degenerate Orthogonal Polynomials and Painlevé II

Using WKB analysis, 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 set the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t\in \mathbb {C}$$\end{document}t∈C for which the spectrum of the quartic anharmonic oscillator in the complex plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d}^{2} y}{\textrm{d} x^{2}} - \left( x^4 + tx^2 + 2Jx \right) y = \Lambda y, \end{aligned}$$\end{document}d2ydx2-x4+tx2+2Jxy=Λ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 (monic) polynomials.


Introduction and results
The second Painlevé equation is an ODE in the complex domain given by d 2 u dt 2 = 2u 3 + tu + α, (PII) with α ∈ C and t ∈ C. It admits rational solutions when α = n ∈ Z and this was recognized by Vorob'ev and Yablonskii in two separate papers [Vor65,Yab59]. For α = 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 (VY) 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 d 2 y(x) dx 2 − x 4 + tx 2 + 2Jx y(x) = Λy(x) (1.1) y(x) → 0 as x → ∞ and arg(x) = π, ±π/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)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 ∈ 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 ∈ C, the corresponding eigenvalues of these quasi-polynomial eigenfunctions are the zeros of a secular equation of degree n + 1 in the shifted variable λ = Λ − t 2 4 . There are particular discrete values of t ∈ 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 ∈ N, these particular values of the parameter t ∈ 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). 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 quasipolynomials solutions of (1.1), that is, solutions of the form where θ(x) = x 3 3 + tx 2 (1.3) 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].
and p(x) is a polynomial. We may sometimes use the notation θ(z; t) to emphasize the dependence of θ on the parameter t. A substitution into (1.1) leads to d 2 p dx 2 + (2x 2 + t) dp dx − 2(J − 1)xp = λ p, with λ = Λ − t 2 4 .
Upon setting J = n + 1, one notices that the ODE d 2 p dx 2 + (2x 2 + t) dp dx − 2nxp = λ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) × (n + 1) matrix M n (t) with respect to the usual basis (1, x, x 2 , . . . , x n ) t , namely (1.5) We observe that the matrix M n (t) is a 4-diagonal matrix. Eigenvalues of this matrix correspond to the eigenvalues λ 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, λ) := det(M n (t) − λI).
Now consider the discriminant of C n (t, λ) with respect to λ, i.e.
D n (t) := Disc λ (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 α = n ∈ 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: The VY polynomials are constructed recursively [Vor65,Yab59] as follows: 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) = − 4 3 n(n+1)/6 n k=1 (2k − 1)!! S (n,n−1,...,1) − 3 4 1 3 t, 0, 1, 0, 0, . . . . (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) = Disc λ (C n )(n 2 3 s) in (1.7) and the roots of the rescaled Vorob'ev-Yablonsky polynomials Y n (n 2 3 s) form two coinciding lattices as n → ∞.
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.
Remark 1.2. When α = 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 (1.10) Therefore for fixed n ∈ 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. The boundary problem (1.1)-(1.2) has solutions if and only if J = n + 1 ∈ N; in this case Λ takes at most n + 1 values for each t ∈ C. See Proposition 3.2. For fixed J ∈ N these values (t, Λ) are called the Exactly-Solvable (ES) spectrum.
2. For given J = n + 1 ∈ N and (t, Λ) in the ES spectrum, the corresponding solution of (1.1) is of the form y(x) = p n (x)e 2θ(x;t) with p n (x) a polynomial of degree n and θ(x; t) = x 3 3 + tx 2 . This polynomial is degenerate orthogonal in the sense that κ ∞3 ∞1 + κ ∞5 ∞3 p n (x)x j e 2θ(x;t) dx = 0, j = 0, . . . , n, (1.11) where κ andκ are some constants. Here ∞ k denotes a contour that tends to infinity with asymptotic direction of argument iπ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. In fact more is true: given t ∈ C suppose that κ, κ are such that (1.11) admits a nontrivial polynomial p n for solution. Then y(x) = p n (x)e x 3 3 +tx is a solution of the boundary problem (1.1)-(1.2). This is proved in Theorem 3.6. 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. 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: ∞2j+3 ∞2j+1 p n (x) 2 e 2θ(x;t) dx = 0, 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 e 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 The reason for this scaling is that it yelds the same WKB-type equation with a small parameter ℏ and a n-independent quartic potential ℏ 2 d 2 y dz 2 − Q(z; s, E)y = 0, Q(z; s, E) = z 4 + sz 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. 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 ℏ. This section does not contain new results and is mostly a preparation for the two subsequent ones.
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 VY polynomials in terms of certain contour integrals.
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.18).
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.18) and (5.19), that we report hereafter (to leading order): , Im(τ (s, E)) > 0.
where τ 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π(n + 1) while the sum of the three logarithms is 2iπ (principal determination) due to the definition of τ (s, E). On the other hand, the quantization conditions for the Vorob'ev-Yablonskii zeroes, to the same order of approximation, read 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ℏ 2 3 plane (where ℏ −1 = n + 1 2 for the VY polynomials and ℏ −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 = e iπ 3 j 3 2 1/3 . Within this region the number of points is O(n 2 ) so that the relative spacing is O(n −1 ).
Our results indicate that near the origin the two lattices have a discrepancy (in the re-scaled s-plane) of order 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 ̸ = 0, their discrepancy would be only O(n −1 ). This would seem at first inconclusive because the lattices have a separation already of order 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. 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 − 2 3 instead of the more natural (n + 1) − 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.

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: where α ∈ C is a parameter, u = u(t), v = v(t), w = w(t) are meromorphic functions of t and The compatibility condition ∂ x ∂ t Φ = ∂ t ∂ x Φ of the system (2.1) is equivalent to which in turn yields the following system of ODEs for u, v, w: Eliminating w from the system u(t) satisfies (PII) with parameter α; 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 α. Then near t = a we have the following Laurent series expansion of (2.6): where b is arbitrary and c ̸ = 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: where, again, b is arbitrary and c ̸ = 0 is a constant of integration.
Proof. It is well known that all poles of PII are simple and with residue ±1 (see e.g. [IN86]). For every such pole t = a, and arbitrary b ∈ 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. ■

Gauge transformation
With a particular gauge transformation we can convert the matrix ODE ∂ ∂x Φ = A(x, t; α)Φ 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 × 2 matrix ODE system Φ (2.11) We wish to turn this into a scalar ODE of the form y ′′ − V (x)y = 0, which we can do by using the gauge transformation where the function V (x) is: (2.14) Thus the matrix ODE system W x = M (x)W (x) is then equivalent to the scalar ODE 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: (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 ±1 of the PII solution u(t) with parameter α. Then the lim t→a V (x, t) exists and we have for a pole of residue +1 and for a pole of residue −1 Proof. The result is obtained from a straightforward computation where we substitute the Laurent series expansions (2.7) into (2.15) and take the limit as t → a. The details are left as exercise (it is helpful to use a computer algebra program for this purpose). ■ 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 ′′ (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).
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.

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: We will call the set (t, J, λ) 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.
with p(x) a polynomial of degree n if and only if J = n + 1 and λ = Λ − t 2 4 is an eigenvalue of the operator acting on the space of polynomials of degree up to n.
Proof. Substituting y(x) = p(x)e θ(x) in the ODE (3.1) gives an equivalent differential equation for the function p(x): and L J as in (3.4). One can readily see that if J = n + 1 then the operator L n+1 in (3.5) preserves the space of polynomials of degree at most n and then Λ 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 L n+1 preserves the space of polynomials of degree n and Λ (and the corresponding λ as per (3.5)) is an eigenvalue of the corresponding finite dimensional operator.

Stokes phenomenon
Proposition 3.2. The eigenvalue problem (3.1) with the boundary conditions (3.2) require that J = n + 1, n = 0, 1, . . . and that λ in (3.5) is an eigenvalue of the operator 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 W(x): 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 We define the Stokes sectors as shown in Figure 3, i.e. S k is the sector of opening π/3 centered around the rays of argument π/6 + (k − 1)π/3: Each open sector Ω k = z ∈ C : arg(z) − π 6 −(k − 1) π 3 < π 6 + δ with δ > 0, strictly contains a Stokes' sector S k . Then there exists a unique solution W (k) (x) of (3.6) in Ω k asymptotic to W form (x), namely Since W (k) (x) and W (k+1) (x) have the same asymptotic expansion in Ω k ∩ Ω k+1 , there is a constant matrix S k (called Stokes' matrix), such that where S 2k = 1 0 s 2k 1 and S 2k+1 = 1 s 2k+1 0 1 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 The first six matrices are the Stokes matrices associated with the directions arg(x) = k π 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) = π 3 is also recessive along the directions arg(x) = k π 3 , k = 3, 5. But then the matrix equation (3.10) implies that s 1 + s 3 + s 5 = 0 (3.11) and e 2iπJ = 1 and hence J must be an integer.
To show that J = n+1 is a positive integer and that λ is an eigenvalue of 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)e θ(x;t) with p(x) entire and bounded at infinity by x J−1 . Then Liouville's theorem implies that if J ≥ 1 then p(x) is a polynomial, and for J = 0, −1, −2, . . . 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 λ is the claimed eigenvalue. ■ Next we consider the properties of the operator (3.1) (or equivalently of (3.5)) under the assumption that J = n + 1 ∈ N and λ (respectively Λ) is an eigenvalue. We start from further analysis of the corresponding Stokes phenomenon.
Theorem 3.3. Let J = n + 1, n ∈ N and Λ be an eigenvalue of the boundary value problem (3.1)-(3.2) with eigenfunction the quasi-polynomial F (x) = p n (x)e θ(x) . Let G k , be solutions of the ODE (3.1) linearly independent from F (x) and which can be expressed as (3.12) Here ∞ k indicates that the contour of integration extends to infinity along the direction arg(z) = k π 3 . Then the Stokes phenomenon for the solutions [F, G k ] is given by the following equations: (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: (3.14) The solution is written as: 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: Figure 4: Directions at infinity ∞ k of argument k iπ 3 and the oriented contours γ andγ from ∞ 1 to ∞ 3 and ∞ 3 to ∞ 5 respectively. The shaded regions denote the sectors of dominance of e θ(x;t) where θ(x, t) = x 3 3 + tx 2 . The unshaded regions denote the sectors of recessiveness.
We can split the integral representation of G 2k in (3.16) as follows: for k = 0, 1, 2 we have were the indices are taken mod 6 and we have defined s 2k+1 , k = 0, 1, 2 to be the integral of F (ζ) −2 dζ between ∞ 2k+2 and ∞ 2k , owing to the observation that its contour of integration crosses the Stokes line of argument (2k + 1)π/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 (ζ) −2 dζ. Indeed the integrand F (ζ) −2 dζ 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 nontrivial 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 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 (ζ) 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. ■

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 θ(x) be a fixed polynomial of degree d + 1 and positive leading coefficient, and let Γ = d j=1 s 2j−1 γ 2j−1 , where γ 2j−1 are the "wedge" contours extending from ∞ 2j−1 to ∞ 2j+1 where ∞ ℓ denotes the point at infinity in the directions arg(x) = ℓπ d+1 . Consider the non-hermitian orthogonal polynomials p n (x) satisfying (3.20) We say the polynomial p n is l-degenerate orthogonal if we additionally have ⟨p n , x k ⟩ = 0 k = n, n + 1, . . . , 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 θ(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θ(x,t) dx on the contour where κ = s 1 and κ = s 5 are defined in (3.13). The contour γ is the wedge contour from ∞ 1 to ∞ 3 and γ is the wedge contour from ∞ 5 to ∞ 3 (see Figure 4).
Proof. Consider the quasi-polynomial F (x) := p n (x)e θ(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 to be a fundamental solution of (3.1). We interpret Theorem 3.3 as an 2 × 2 Riemann-Hilbert problem solved by Ψ k (x) in the sectors of opening 2π/3 : The jump conditions (3.13) imply that G k satisfies: The three functions G k , k = 0, 2, 4 define a piecewise-analytic function G(x) with discontinuities along the three rays [0, ∞ j ] with j = 1, 3, 5. The function G(x) satisfies the jump conditions where G ± (x) are the boundary values of the function G(x) on the right and left boundaries of the oriented rays [0, ∞ j ], j = 1, 3, 5. By the Sokhotski-Plemelj formula we can express G as the following Cauchy transform: with γ and γ as in Fig.4. Note that the two contours overlap with the same orientation on the ray [0, ∞ 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)e θ(x) and G k (x)e θ(x) behave like 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 ′ = d dx ): 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 implies the vanishing of the integrals: Γ p n (ζ)ζ k e 2θ(ζ;t) dζ, k = 0, 1, . . . , n.
(3.31) Therefore the polynomials p n (x) are degenerate orthogonal as claimed. ■ 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 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 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 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 . But a simple computation shows that y 2 (x) = 1 + O(ζ) 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 with a ̸ = 0 and again we denote ζ = x − x * . The indicial equation now gives 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 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. ■ We are ready to prove the converse of Theorem 3.4, namely that these degenerate orthogonal polynomials satisfy the Shapiro-Tater ODE (3.1).
Proof. Suppose that p n (x) is a degenerate orthogonal polynomial with respect to the weight w(x) = e θ(x;t) on the (weighted) contour Γ, i.e. Γ p n (ζ)z k e 2θ(ζ;t) dζ = 0 k = 0, 1, · · · , n. (3.40) Define the functions where the leading factor is 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 ′ G, F G ′ are bounded. By the Liouville theorem we conclude that W (x) must be a constant, i.e.
Differentiating this equation gives that F ′′ /F ≡ G ′′ /G. Let us denote by V (x) this ratio; then both F and G satisfy a 2nd order linear ODE of the form: We can rewrite the potential using the defining expression F = p n e θ in terms of the polynomial p n , which gives us: 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: 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 ′′ − 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 ′′ /F is entire. We conclude that d = 1, namely all the zeros x 1 , . . . , x n of the polynomial p n (x) are simple and obtain This gives the potential which coincides with the potential in (3.1) with Λ = t 2 4 + 2 n j=1 x j ■ 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 , . . . , 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: a polynomial of degree n; 2. The polynomial p n (x) (of degree n) is a 1-degenerate orthogonal polynomial with respect to the weight w(x) = e 2θ(x;t) on the contour Γ := κγ +κγ for some values of (κ,κ) ̸ = (0, 0) and with γ andγ as in Fig.4. Furthermore, the coefficients κ = s 1 ,κ = s 5 are expressible in terms of p n (x) by the formulas (3.13).

Repeated eigenvalues
Recall that the characteristic polynomial C n (t, λ) 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. The value of t ∈ C is such that the Exactly Solvable spectrum of (3.1)-(3.2) has a repeated eigenvalue; 2. The value t ∈ C is such that there exist n ∈ N and Λ ∈ 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. there is a quasi-polynomial solution p n (x)e θ(x;t) of (3.1)-(3.2) that satisfies where γ and γ are defined as in Theorem 3.6 (contours from ∞ 1 to ∞ 3 and ∞ 3 to ∞ 5 respectively, see Fig.4).
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 Λ (λ = Λ − t 2 4 ) is an eigenvalue. If this eigenvalue is repeated then the derivative of the characteristic polynomial C n (t, λ) in (1.6) must also vanish.
[(2.) ⇒ (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.) ⇒ (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 λ, then the generalized eigenvector w satisfies the equation (M − λI)w = v.
It follows that the generalized eigenvector equation for a polynomial r(x) (of degree ≤ n) takes the form of the following differential equation for r(x): 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: 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: . (3.59) Before proceeding with the proof we claim and prove that H(x) is bounded by a constant as x → ∞ 1 . To see this it is sufficient to show that the inner integral in the defintion of H(x) tends to zero as O(ζ 2n−2 )e 2θ(ζ;t) (recall that the real part of θ goes to −∞ along the direction ∞ 1 and hence this is exponentially small). To prove this claim we write p 2 n (s) = 2Q 2n−2 (s)θ ′ (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θ ′ . Then integrating by parts we obtain (3.60) To estimate the last integral we "force" again an integration by parts: In the remaining integral the prefactor to the exponential is O(s −2 ) and thus the integral is easily estimated to be O(1)e 2θ . 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 → ∞ 1 . The integral H(x) is bounded by a constant as x → ∞ 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 ∞ 3 and ∞ 5 . We can write, for example for ∞ 3 , The second term above is polynomially bounded near ∞ 3 , by the same argument used to show that H is bounded near ∞ 1 . But since q 0 is exponentially dominant also near ∞ 1 we deduce the condition ∞3 ∞1 F (s) 2 ds = 0. One similarly deduces ∞3 ∞5 F (s) 2 ds = 0, which establishes (3.57).
[(3.) ⇒ (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.62). ■ Remark 3.10. The Theorem 3.9 seems at first sight nothing short of a miracle; indeed once we fix J = n + 1 ∈ N, then the ODE (3.1) has only two continuous parameters t, λ. However the multiple eigenvalue condition apparently involves now three equations which are: (i) the existence of a quasi-polynomial solution which determines λ 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 ∈ N the pair (t, Λ) is in the ES spectrum, then we have shown in Theorem 3.4 that p n (x) is a degenerate orthogonal polynomial, namely, with γ, γ, κ, κ defined in the same theorem. The coefficients κ, κ 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) θ(x;t) is a solution of the boundary problem (3.1) and hence t, J = n + 1 and Λ are in the ES spectrum. Suppose also that Λ = λ + t 2 4 is a repeated eigenvalue. Then the antiderivative of y 2 (x) is also a quasipolynomial, (3.64) Proof. Take the integration basepoint from ∞ 1 and consider the function Standard asymptotic analysis shows that K(x) is of the form O(x 2n+2 )e 2θ(x;t) along the (dominant) direction towards ∞ 0 , ∞ 2 , and ∞ 4 . Moreover towards ∞ 1 we have Thus the function q(x) = K(x)e −2θ(x;t) is polynomially bounded along the sectors containing ∞ j , j ∈ {0, 1, 2, 4}. If we can show that it is also polynomially bounded in the sectors containing ∞ 3 and ∞ 5 then we conclude, by the Liouville's theorem, that it is a polynomial. Consider x tending to infinity along the sector that contains ∞ 3 ; we have x ∞3 p n (ζ) 2 e 2θ(ζ;t) dζ + e −2θ(x;t) ∞3 ∞1 p n (ζ) 2 e 2θ(ζ;t) dζ. (3.67) The second term above is just a multiple of e −2θ(x;t) which is not polynomially bounded near ∞ 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 ∞ 3 . Similarly we can prove that it is polynomially bounded at ∞ 5 and the proof is thus complete. ■

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 ′′ − ℏ −2 Qy = 0, where ℏ is a small parameter. We introduce the scaling for the Shapiro-Tater potential V ST = x 4 + tx 2 + 2(n + 1)x + Λ z = (n + 1) −1/3 x, s = (n + 1) −2/3 t, E = (n + 1) −4/3 Λ (4.1) which implies dy dx = (n + 1) −2/3 dy dz . A similar scaling applies in the case of the Jimbo-Miwa potential (for the poles with residue −1), The result of these scaling is the exact same potential with the identification of E = 7s 2 36 + 10 b, but with different scaling factors. If we set Q(z; s, E) = z 4 + sz 2 + 2z + E, z We now explain how to construct the WKB solutions of (4.6). The formal series ansatz implies that S(z, ℏ −1 ) satisfies the following Riccati equation: (4.9) Comparing each power of ℏ we get a recursive relation for the coefficients S k (z) of the formal series: S n S m (k ≥ −1). (4.11) The first three coefficients S k (z) are: (4.12) The choice of sign for S −1 (z) = ± 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; ℏ): Since our potential Q(z) is independent of ℏ, it follows that S odd only contains odd powers of ℏ, namely: We will only care about S odd since it can be shown [KT05] that S even can be written in terms of S odd as: Another important fact is the following: the differential S odd (z)dz has a pole at z = ∞ which comes solely from the term S −1 (z) = √ Q(z); namely, S 2k+1 (z) = O(z −2 ) as z → ∞. 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 T = {τ 0 , τ 1 , τ 2 , τ 3 } of Q(z) as turning points.
Definition 4.1 (WKB solutions). The WKB solutions to (4.6) are formal power series in ℏ given in terms of S odd in (4.15). We give two different normalizations that will be used throughout this paper.
• Near a turning point τ of the potential Q(z) we define the normalized WKB solutions to be: • Near infinity we define the normalized WKB solutions to be: where all logarithms are principal. Notice that R(z; ℏ) is the anti-derivative of S odd (z, ℏ) that does not have a constant term in the expansion as |z| → ∞.
Remark 4.1. The integral in the exponent of ψ (τ ) ± is to be understood term-wise in each coefficient of the powers of ℏ. Additionally, S odd (z, ℏ) is multivalued on C with branch points at the zeros T := {τ 0 , τ 1 , τ 2 , τ 3 } of Q(z). Therefore the integral should be considered on the elliptic compact Riemann surface Σ obtained from the affine curve Σ = (w, z) ∈ C 2 : w 2 = Q(z; s, E) (4.20) by adding two points P ± ∞ at infinity. The projection π : Σ → C from the Riemann surface Σ to the extended complex plane C, maps (w, z) → z. The projection π realizes Σ as a double cover of C ramified at the zeros of Q (turning points). The pre-image of any point z ∈ C are the two points π −1 (z) = (z, ±w) on the two sheets of the Riemann surface where the numbering of the sheets is such that (z, w = 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 τ to z by defining: where the contour γ(z) lives in the Riemann surface Σ and goes from z to z, where z lies on the second sheet and the second sheet is reached by encircling the branch point τ . This integral is a well defined convergent integral. Finally, it can be shown by induction that all the correction terms S 2j+1 (z)dz, j ≥ 0 are differentials on the Riemann surface Σ 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 Proof. Using (4.3) we see that (4.23) Now, we can write (we drop the indication of the dependence on x, s, E) (4.24) The periods of the second term in (4.24) vanish because this gives an exact differential; the first term reads 1 48 (4.25) Integrating (4.24) along γ and using the identity (4.25) completes the proof. ■

Stokes graphs and connection formulae.
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 Q(z) and introducing the notion of Stokes' graph. To minimize confusion when performing integrations along the branch cuts we will make the following explicit definition. to be the integral of Q(z; s, E) along the + side of the branch cut oriented from τ to τ . 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)dz 2 where one of the end-points is a turning point. In other words, in a local coordinate z of Σ it is a curve emanating from a turning point τ and satisfying The support of the curve is independent of choice of determination of √ Q. Furthermore: 1. The orientation of a Stokes curve is defined by the direction in which Re z τ Q(u)du is increasing.
2. The directions near the point at ∞ are indicated by ⊕ or ⊖ if the function Re z τ Q(u)du is increasing or decreasing, respectively, along the corresponding Stokes curve. We say the Stokes curve is oriented towards ⊕ or oriented away from ⊖, respectively. Furthermore, we decorate the Stokes graph with the following additional contours.
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. Ideal paths. We draw arbitrary (smooth) paths connecting the different ∞ 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 ⊕ 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. 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 4 .
The construction in Def. 4.5, under the Assumptions 4.5, is crafted so that each Stokes region D has precisely one and only one turning point on its boundary.
In each Stokes region D one can select a fundamental basis of solutions of the ODE (4.6) as we now explain. Consider a Stokes curve γ originating at the turning point τ and oriented towards ⊕ (with similar considerations for the ⊖ curves ): since the real part of z τ Q(u)du is increasing, there is a region around γ near ∞ where the real part is positive; then the formal solution ψ (τ ) − is recessive (i.e. exponentially small as ℏ → 0 + ). Then there is a unique solution Ψ(z; ℏ) which is asymptotic to this ψ + , it is not sufficient to examine its asymptotic behaviour near γ because its asymptotics there is dominant (i.e. exponentially large as ℏ → 0 + ). However, the same Stokes region must be bounded also by either a ⊖ trajectory or a branch-cut. In the former case, in a neighbourhood of the ⊖ trajectory the formal solution ψ If, instead, the other boundary is a branch-cut, we need to consider the Stokes region, D ′ on the other side of the cut: in this region, due to having crossed the branch-cut, the formal solution ψ (τ ) + is now somewhere recessive and this allows to fix Ψ (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. ± , 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: (4.28) Furthermore, let D ℓ , D r be two adjacent Stokes regions separated by the Stokes curve γ oriented as in Def. 4.5 with D ℓ on the left and D r on the right of γ. Then the corresponding solutions Ψ D ℓ ± , Ψ Dr ± are related by: The relationship between solutions in regions separated by an ideal path is simply a different scaling as given by the following proposition. where R(z; ℏ) is defined in (4.19). Note that w τ (ℏ) is a constant with respect to z since d dz w τ (ℏ) = 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. (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 τ 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) ∈ C 2 the potential Q(z; s, E) may have a different different Stokes lines in C but the underlying Stokes graph will be topologically identical to one of the types depicted.
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 ∞. 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 Σ relative to the choice of branch as in Def. 4.2. 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 τ j determines a triangle inside the hexagon by connecting the asymptotic directions at infinity with a Stokes curve originating from τ 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π/6. From configuration Z we obtain 6 distinct configurations, 3 of them corresponding to a Z 3 (i.e. 2π/3) rotation, and another 3 corresponding to a reflection followed by a Z 3 rotation. Finally from configuration E + we obtain 6 distinct configurations corresponding to a Z 6 (i.e. 2π/6) rotation. Finally, 2 + 3 + 3 + 6 = 14 as claimed. ■

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 (c) RHP for graphs of type Z. To each fourth-root branch cut we associate the jump matrix Y = −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. Stokes curves: along each Stokes curve oriented towards ⊕ and away from ⊖ we assign the jump matrices B and R (respectively) The Stokes curves are coloured blue or red, respectively, in Fig. 6. 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: where j, k ∈ {0, 1, 2, 3} and τ j , τ k are turning points. These contours are denoted by a brown dashed line in Fig. 6.

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: (4.41) where j ∈ {0, 1, 2, 3}, τ j is a turning point and R(z; ℏ) is the particular antiderivative in (4.19). These contours are denoted by a black dashed line in Fig. 6.
Remark 4.8. We observe that w j (ℏ) is a constant in x. It can be thought as the regularized integral τj ∞ S odd (z, ℏ)dz. This construction gives us, up to rotation and reflections, three distinct WKB Riemann-Hilbert problems as shown in Fig. 6.
In order to simplify the upcoming computations we make some notational definitions. Given turning points τ j and τ k connected by a branch cut we denote with the determination of S 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 ℏ. Now we obtain the following results about the Stokes matrices for each configuration shown in Fig. 6.
Theorem 4.9. The Stokes matrices S j = 1 0 s j 1 , j = 0, 2, 4, S j = 1 s j 0 1 , j = 1, 5, S 3 = 1 s 3 0 1 e 2iπ ℏ σ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 ξ jk as defined in (4.40), (4.41) and (4.42) as follows: Configuration D Configuration E Configuration Z s 0 = −ie 2w1 (ξ 10 ξ 30 + ξ 10 + 1), s 2 = −ie 2w2 (ξ 20 ξ 10 + ξ 20 + 1), s 0 = −ie 2w0 (ξ 01 ξ 12 ξ 32 + ξ 01 ξ 12 + ξ 01 + 1), s 0 = −ie 2w0 (ξ 01 ξ 12 + ξ 01 + 1), In the configuration D the following identities hold e 2(w3−w1) = ξ 10 ξ 30 , e 2(w1−w2) = ξ 10 ξ 20 , ξ 10 ξ 20 ξ 30 = e 2iπ ℏ . (4.44) In configuration E and Z the following identity holds (4.45) Proof. We compute the Stokes matrices associated to Configuration D of the Riemann-Hilbert problem in Fig. 6. Consider the vertex ∞ 0 in Configuration D, there are seven edges incident on it. The Stokes matrix S 0 is given by the clockwise (about the vertex ∞ 0 ) product of all the jump matrices corresponding to the edges incident to ∞ 0 as follows S 0 = W −1 3 BV 30 BV 10 BW 1 (4.46) Notice that the jump matrix V 30 follows from the fact that we have where z ± are the boundary values of S odd (z, ℏ) on the left and right side of the segment [τ 0 , τ 3 ] oriented from τ 3 to τ 0 . Therefore we have We observe that where the integral from τ 1 to τ 3 is with the determination in Def.4.2. Therefore we obtain which is equivalent to the first identity in (4.44). Similar calculations compute the remaining Stokes matrices for each possible configuration. The extra e 2iπℏ −1 σ3 in the form of S 3 is due to our choice of branch-cut for the logarithm in (4.19). The fact that the product of all Stokes matrices is trivial follows by construction. The second equation in (4.44) for configuration D follows from the definitions of ξ ij in (4.42) and w j in (4.41). The last equation follows from the fact that γ 10 + γ 20 + γ 30 is homologous to Γ ∞ , the contour at infinity with winding number equal to one. The residue at infinity of the (formal series in ℏ) S odd (z, ℏ)dz is 2iπ/ℏ 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 ℏ-series expansion of S odd vanish to order O(z −2 ) and hence do not contribute to the residue). Regarding the relation (4.45) 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 γ 01 + γ 23 is homologous to Γ ∞ . Thus in a similar way as the previous case ■ 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 (4.51) 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 s 0 =i(ξ 10 ξ 30 + ξ 10 + 1), s 1 = i s 2 =ie 2w2−2w1 (ξ 20 ξ 10 + ξ 20 + 1) = i 1 ξ 20 ξ 10 (ξ 20 ξ 10 + ξ 20 + 1), s 3 =ie −2w2+2w1 = iξ 10 ξ 20 s 4 =ie 2w3−2w1 (ξ 30 ξ 20 + ξ 30 + 1) = iξ 30 ξ 10 (ξ 30 ξ 20 + ξ 30 + 1) (4.52)

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 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.18) from an application of Theorem 5.4.

The Jimbo-Miwa case
The parameters s, E in the potential Q(z; s, E) i.e. the parameters t, Λ 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]. Thus, to find the positions of a pole, we can find for which values of a, b in (2.18) (with α = 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 ℏ −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 ℏ: normalized near a turning point τ .

(5.2)
In terms of the Vorös symbols we obtain the following leading order estimates: where k 1 , k 2 , k 3 are integers and ℏ = (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.: 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.52) the conditions: s 0 =i = i(ξ 10 ξ 30 + ξ 10 + 1), s 1 = i s 2 =i = i 1 ξ 20 ξ 10 (ξ 20 ξ 10 + ξ 20 + 1), s 3 = i = iξ 10 ξ 20 s 4 =i = iξ 30 ξ 10 (ξ 30 ξ 20 + ξ 30 + 1), The only solution of the above system is ξ 20 = ξ 10 = ξ 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. ■

The Shapiro-Tater case
Theorem 5.3. Suppose that the quasi-polynomial y = p n (x)e θ(x;t) is a solution to the boundary problem (3.1)-(3.2) with J = n + 1. Then the "exact" Fock-Goncharov parameters ξ jk in (4.42) satisfy one of the following systems, depending on the indicated Stokes graph configuration.

(5.9)
Furthermore all other configurations cannot occur.
Proof. In Proposition 3.2 it was shown that if J = n + 1 and (t, Λ) 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.52) we have expressed the Stokes parameters in terms of the exact Fock-Goncharov parameters ξ 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 ℏ = 1/(n + 1) and hence e 2iπ ℏ = 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 Z 3 rotation or reflection, and those obtained from D or E by a reflection around the imaginary axis followed by any Z 3 rotation. ■

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 ξ 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 2 n e 2θ between ∞ 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 (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 ψ τj ± in Def. 4.1 in different regions in accordance with Thm. 4.3. Namely where B is one of the labelled regions in Fig. 8 ad τ j is one of the turning points in its boundary. We observe that to compute I 13 we can equivalently compute the integral of Ψ (A1) + because this function is proportional to p n (x) 2 e 2θ(x;t) since they are both recessive in the direction ∞ 1 . Thus, the main aim of this section is to prove the following Theorem. 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 , Im(τ (s, E)) > 0. (5.13) 2. In configuration E as in Fig. (6) we have 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 Ψ (A1) + 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 τ 1 , τ 0 , τ 2 . The unbounded integrals are then along steepest descent paths for the integrand and can be neglected. In the integrations along [τ 1 , τ 0 ] and [τ 0 , τ 2 ] we observe that, as a consequence of the Riemann-Hilbert problem 4.7 we can express Ψ where we have used already the fact that we are on the ES spectrum so that Theorem 5.3 applies and the Stokes matrices S j , j ∈ {0, 2, 4} are trivial. Computing the square of Ψ (A1) + , 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 where the boundary values of √ 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 − τ ) − 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 O(ℏ

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 In terms of the Voros' symbol we then have, to leading order 2(n + 1) 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π(n + 1) while the sum of the three logarithms is 2iπ (principal determination) due to the definition of τ (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 Both conditions (5.18), (5.19) 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.

The elliptic region
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. 5 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.18) imply that, to leading order as n → ∞, the three periods We take these equations to leading order; τ0 τj Q(z + ; s, E)dz = −iπµ j , µ 1 + µ 2 + µ 3 = 1, µ j ∈ R. The central triangular region corresponds to the D topology; the three "corridors", adjacent to the sides, in the directions arg(s) = 0, ±2π 3 correspond to the three Z 3 rotation of a E, while in the six remaining regions we find the various Z 2 × 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. This means that the differential √ Qdz 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.21) are imaginary being the Boutroux condition. Under this condition we know that the critical graph, C, is connected. This is, by definition, the graph determined by the "vertical critical trajectories 6 " [Str84], namely, the maximal integral curves of the direction field Re Q(z; s, E)dz = 0 (5.22) that issue from each of the turning points. Observe that the integral curves of (5.22) form a foliation that is orthogonal to the Stokes' graph (4.27), in the sense of the Riemannian metric dσ 2 := |Q(z; s, E)||dz| 2 (with conical singularities at the turning points). Briefly this can be seen by considering the real-valued function Φ(z) := Re z τ0 wdz , which is easily shown to be -harmonic on the Riemann surface Σ (4.20) minus the two points above z = ∞; -odd under the holomorphic involution (w, z) → (−w, z) (i.e. under the flip of sign of the square-root).
The connectedness of Φ −1 ({0}) on Σ (and hence also of its z-projection) follows from the fact that Φ 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 Φ −1 ({0}) is a well-defined subset of the complex z-plane (i.e. of the projection from Σ to the z-plane). It is then not hard, following similar arguments to those used in analyzing the topology of the Stokes' graph, that C can only have three different topologies corresponding to the three configurations D, E, Z (and "dual" to the Stokes' graph). The graph C is however also a metric space with the distance induced by the Riemannian metric dσ 2 ; in fact the geodesic length between two turning points (along the graph) is precisely the (absolute value) of the purely imaginary integral of τj τi √ Qdz along any path homotopic to the edge of C that connects the two turning points. In other words the numbers π|µ 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 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 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.

Analysis of the two lattices
Both lattices involve implicit equations for the parameters (s, E) via the periods of the differential Q(z; s, E)dz.
We introduce a canonical basis of cycles A, B of the elliptic Riemann surface Σ as in the Fig.10  The Jacobian determinant is a constant as we prove in the next lemma.
Lemma 5.5. Let A and B be the canonical homology basis as in Fig. 10  where we have used that . The lemma is useful in that it allows us to explore the geometry of the quantization conditions; Proposition 5.6. Let (s 0 , E 0 ) correspond to the first-order quantization conditions (5.18) or (5.19) in the bulk, namely, m j /n ≃ c j ̸ = 0. Then the neighbour points in the s-plane form a slowly modulated hexagonal lattice in the sense that the six closest neighbours of s 0 are where ω and ω ′ are the half periods of the holomorphic differentials in (5.33) and ∆m j ∈ {−1, 0, 1}, |∆m 1 + ∆m 2 | ≤ 1, |∆m 1 | + |∆m 2 | ≥ 1.
We are now going to show that two quantization conditions yield the same lattices to order ℏ 2 .
Then, using the ODE, we get Thus we seek polynomials such that The latter equation for B has a unique polynomial solution obtained by inverting the (finite-dimensional) linear operator on the space of polynomials of degree ≤ n. The relationship of the leading coefficients follows by plugging into (A.1). ■ Proposition A.2. Let f (x) be an analytic function in a neighbourhood of the origin.
where √ x is analytic in C\(−∞, 0] and positive for x > 0 and √ x + denotes the boundary value from the upper half plane.
[2] Let a, b > 0. Then (A.12) Here ω = e 2πi 3 and √ x and √ x + are defined as above and Proof. We will consider f (x) = x n , with the proof being completed simply by summing the series.
[1] Then we need to compute b −a Ai 2 (ℏ − 2 3 x)x n dx. We can use Lemma A.1 and setting α = aℏ − 2 It is then clear that the evaluation at β = ℏ − 2 3 b yields exponentially small terms thanks to (A.15). For the evaluation at −α = −ℏ − 2 3 a we use instead (A.16); thanks to (A.5) we see that Finally, we observe that [2] All these integrals involve oscillatory direction for each of the arguments. The computation is of entirely similar nature, where it is particularly important to pay great care at the relative phases. Using Lemma A.1 we need to estimate integrals of the form (Ai j (ξ) := Ai(ω j ξ)) At the point α = ℏ − 2 3 a a direct computation using the asymptotic properties (A.15), (A.16) yields: Similarly one finds with q jk , r jk as in (A.13). Now the integral (A.12) with f (x) = x n becomes Putting together the two terms we obtain the proof. ■ As a consequence of the Lemma A.1 we can estimate the contribution of integrals involving Ψ 2 in a neighbourhood of a turning point.
Theorem A.3. Let τ be a simple turning point for the potential Q(x), γ d a steepest descent ⊖ path and γ o the oscillatory path on the opposite side of τ : we choose two points τ + ∈ γ d , τ − ∈ γ o , in the neighbourhood of τ at finite positive distance. Let Ψ be a solution of the ODE Suppose that Ψ is asymptotic to the recessive formal solution ψ with the branch-cut of √ Q running along the oscillatory path and the determination the one that has regative real part on the descent path γ d .
The relationship between the two equations is as follows The formal series ξ(z; ℏ 2 ) is Borel-resummable to an analytic conformal mapping whose asymptotic expansion as ℏ → 0 + coincides with the formal series [KT05]. The arc [τ, τ + ] on the contour γ d is mapped to a segment [0, ξ + ] ⊂ R + in the ξ-plane, and similarly the arc The function Ψ is the recessive solution on γ + and hence it must map to a multiple of the Ai function within a whole neighbourhood of τ . Thus the integral to be estimated translates to where C(ℏ) is an appropriate proportionality constant that depends on the chosen normalization of Ψ. Our choice is to normalize Ψ so that it is asymptotic to the recessive formal solution ψ (τ ) + . To fix C(ℏ) we can consider the asymptotics of the Airy function at τ + (A.15) and compare it to the WKB solution. Indeed we find that so that C(ℏ) = 4ℏ Now from (A.25) we get and hence, to within O(ℏ 2 ), we have Using the estimate of C(ℏ) and substituting into (A.27) we obtain the statement of the theorem. ■ A.1 Proof of Theorem 5.4 We will consider in detail only the case D because it contains all the intricacies.
Since we want to impose that the integral I 13 vanishes, it suffices to determine the integral of a function that is proportional to y n = p n e 2θ . Keeping this in mind we then observe that in the region (A 1 ) along the Stokes curve (τ 1 , ∞ 1 ) the function is recessive (as n → ∞ and as x → ∞ 1 ) and hence it must be proportional to Ψ (A1) + (z; ℏ) and thus asymptotic to ψ (τ1) + . We will then forget about y n and integrate instead the function Ψ (A1) + (z; ℏ) 2 along γ: recall that this function is entire since it is the solution of the ODE (4.6). Let us refer to Fig. 8. We choose the contour γ to follow the steepest descent from ∞ 1 to τ 1 , then the branch-cut to τ 0 , then to τ 2 and the descent path to ∞ 3 . The branch-cut will be chosen to follow the anti-Stokes curve in a neighbourhood of each turning point.
We then split our contour γ = 7 j=1 γ j in several pieces Observe next that the solution Ψ := Ψ (A1) + is recessive both on the Stokes line ∞ 1 → τ 1 as well as on the Stokes line τ 2 → ∞ 3 because we are already on a point of the ES spectrum, where the Stokes matrices satisfy S 0 = S 2 = S 4 = 1. More precisely it follows that Ψ Contributions of γ 1 , γ 7 . On both these contours the integral is exponentially small as ℏ → 0 + since Ψ = Ψ The boundary value is the − because the orientation of the branch-cut that we have chosen is the one egressing from the turning point, while in Theorem A.3 it was the one towards the turning point (and the boundary value was the + one).
Contributions of γ 3 , γ 5 . Along either paths we are away from the turning points and can substitute the asymptotic expressions in terms of the WKB solutions. We can use the Riemann-Hilbert problem 4.7 to see that the function Ψ = Ψ z ∈ (B ℓ 1 ) e v12 e v20 ψ (τ0) z ∈ (C r 1 ) e v12 ψ (τ2) The different expressions in (B ℓ,r 1 ) are due simply to the different choice of normalization point for the WKB formal solutions, but we can use either of the two expressions also as asymptotic expansion in the other. Ditto for the pair (C ℓ,r 1 ). Consequently we can estimate the contribution of the whole γ 3 using the asymptotic expression for (B r 1 ) as follows: where the − boundary value in the last integral is due to the fact that the regions (B ℓ,r 1 ) lie on the − side of the branch-cut, with the orientation that we have chosen. In particular this contribution together with the contribution (A.33) yield a single integration from τ 1 to τ + 0 . Similar considerations apply to the integral along γ 5 which yield To see that the first integral in (A.36) is sub-dominant, for example consider the integration of (ψ (τ1) + ) 2 . We can deform the path of integration (at fixed endpoints) into the region (C 3 ) where the formal solution is recessive. Then a simple estimate shows that the only contribution come from the neighbourhood of the endpoints of integration (at which points the real part of the exponential is zero). Then we can use Laplace's method to easily estimate their contribution of order 7 O(ℏ 2 ). A similar reasoning applies to the integration of (ψ (τ1) − ) 2 where instead we deform the contour within the region (B 1 ).
Contribution of γ 4 . We are going to show that this integral is estimated as follows: .
(A.40) 7 Note that the function being integrated has already an ℏ in front, and the integral contributes another order O(ℏ).
Refer now to Fig. 12: the term marked (a) is recessive in the blue-shaded region, the term (b) in the pink-shaded region while (c) is oscillatory throughout. Consequently we will split the integration of Ψ 2 into three paths from τ + 0 to τ − 0 and integrate (a) along the blue path, (b) along the red path and (c) along the green path.
The whole contribution of (b) is then sub-leading and of order O(ℏ 2 ) where the main contributions come from the neighourhoods of τ ± 0 . The contribution of (a) near τ + 0 and ξ 0 similarly is of order O(ℏ) 2 (with the main contribution coming solely from the neighbourhood of τ − 0 ). The integration from ξ 0 to τ − 0 is achieved by Theorem A.3 and we have We then have the contribution of (c) along the green path; for this we need to use a similar reasoning as in the proof of Theorem A.3. The function Ψ (B r 1 ) + is recessive in the blue-shaded region of Fig. 12; thus in the coordinate ξ(z; ℏ) along the direction arg ξ = 0 and is therefore proportional to Ai(ℏ − 2 3 ξ)/ √ ξ ′ after the change of coordinates. Similarly the function Ψ (B r 1 ) − is recessive in the pink-shaded region and hence proportional to Ai(ℏ − 2 3 ω 2 ξ)/ √ ξ ′ , with ω = e 2iπ/3 where above each term we have recalled with contribution it comes from. We now see that the contribution indicated with (γ 4 ) ′ and (γ 5 ) combine to give a single integral: this is due to the fact that e −2v10 = ξ −1 10 and from (5.8) we see 1 + ξ −1 10 = −ξ 30 . Then, using the identity τ0 τ3 S odd (z + ; ℏ)dz = 2iπ(n + 1) − τ0 τ1 S odd (z + ; ℏ)dz −

B Numerical verifications
We performed numerical verifications of the formulas in view of the adage, valid in many walks of life, of "trust but verify". For the JM case we proceed as indicated below.
1. Compute numerically in arbitrary precision the list of zeros, Z n , of the VY polynomial Y n (t), i.e., the poles of the rational solution u n (1.8).
For each of a ∈ Z n , compute the corresponding value of the parameter b appearing in (2.8). Thus we can define, for each a ∈ Z n a corresponding value of Λ = 7a 2 36 + 10b as in (2.17). Let us call L n the list of the corresponding values of Λ.
2. Let S n := Z n /(n + 1 2 ) 2 3 and E n := L n /(n + 1 2 ) 4 3 as per scaling (4.2). For each pair (s ℓ , E ℓ ) in S n , E n we construct the corresponding potential Q(z; s ℓ , E ℓ ) = z 4 + s ℓ z 2 + 2z + E ℓ and evaluate numerically the Voros symbols including the subleading term, using Prop. 4.2, along the relevant cycles of the Riemann surface of √ Q. The numerical verification consists in checking (2n + 1) Q(z; s ℓ , E ℓ )dz + 1 2n + 1 S 1 (z)dz ≃ iπ + 2iπm γ , (B.1) i.e. an odd multiple of iπ. We tested up to n = 26 and the numerics indeed supports the formula. It is interesting to observe that the leading order computation yields (unsurprisingly) less accurate approximations of odd multiples of iπ, in particular in that it has a small but still non-negligible real part.
For the ST case we proceed similarly as follows.  Figure 13: The ES spectrum (black) and the VY zeros (red) for n = 40. On the left the zeros of the ST discriminant are scaled as (n + 1) − 2 3 , which is the "natural" scaling. On the right they are scaled by n − 2 3 .
To effect this test, we choose a point, s 0 in the elliptic region and find the closest points of either lattices to s 0 . Let ∆ n (s 0 ) denote the distance between these two points belonging to the different lattices; we call it the "local discrepancy". Then we plot ln ∆ n (s 0 ) against ln n. If we choose s 0 at finite positive distance from the origin, the slope of this graph is −1 while for s 0 ≃ 0 the slope is −2. This plot of local discrepancy ∆ n (s 0 ) is reported in Fig. 14 in the two regimes of s 0 ≃ 0 and s 0 ̸ = 0, verifying that the decrease of the discrepancy is O(n −2 ) in the first case and O(n −1 ) in the second.  Figure 14: The local discrepancies ∆ n (s 0 ) against n in a log-log plot for two values of s 0 . The red dots correspond to the discrepancy with the "wrong" scaling n − 2 3 . Although the "wrong" scaling is better, the rate of convergence is the same in the local lattices. In the left picture s 0 = 0 and the slope is −2 (the dotted line is plotted for reference). In the right picture s 0 = 1 + i and the slope is −1 as expected.