High order perturbation theory for difference equations and Borel summability of quantum mirror curves

We adapt the Bender-Wu algorithm to solve perturbatively but very efficiently the eigenvalue problem of"relativistic"quantum mechanical problems whose Hamiltonians are difference operators of the exponential-polynomial type. We implement the algorithm in the function BWDifference in the updated Mathematica package BenderWu. With the help of BWDifference, we survey quantum mirror curves of toric fano Calabi-Yau threefolds, and find strong evidence that not only are the perturbative eigenenergies of the associated 1d quantum mechanical problems Borel summable, but also that the Borel sums are exact.


Introduction
It is usually very rare to have an exact solution to a quantum mechanical problem. Most quantum mechanical systems are either solved numerically or using some approximation scheme, typically relying on some small parameter. The most famous and general approximation scheme is the perturbative expansion around the Planck constant . Perhaps surprisingly however, the generic expansion coefficients grow factorially with the order, rendering the series badly divergent, which calls into question the meaning of the perturbative expansion itself. Enter the resurgence theory ofÉcalle, an idea that a proper definition of the complete solution requires the inclusion of terms non-perturbative in the coupling which, upon proper definition, are believed to cure all ambiguities and pathologies associated with the pathological series expansion. See for example [1,2], and more comprehensive references in [3].
Early connection of this interplay were noticed independently by Zinn-Justin and Bogomolny, when considering the contributions of instanton-anti-instanton pair to the partition function [4,5]. They proposed that such a pair is ill-defined itself, and upon a certain -somewhat ad hoc-prescription (the Bogomolny-Zinn-Justin or BZJ prescription in the literature), contains an ambiguity of the same kind that exists in the Borel summation of the perturbation theory. They showed that indeed this ambiguity between perturbative and non-perturbative contributions cancel to leading order. Recently however, the ad-hoc BZJ prescription found an explanation in terms of Lefshetz thimble decomposition [6][7][8][9][10][11]. Furthermore these ideas led to methods for solving the Schrödinger equation, such as uniform WKB [12][13][14][15][16], exact WKB [17][18][19][20]. We also mention a fresh perspective on the problem of Borel summation [21,22] in which it was shown that in quantum mechanics perturbation theory can be recast in a form which completely captures nonperturbative physics.
On the other hand, resurgence in quantum field theory was discouraged due to the discovery of another source of factorial growth of the perturbation series: the 't Hooft renormalons [23], which occurs because of the running of the coupling, and has no analogue in the quantum mechanical systems and ordinary differential Schödinger equation. Furthermore, the ambiguities coming from the renormalons did not seem to be a result of semiclassical configurations such as instantons. This stymied works in this direction for a long time, and it became widely believed that resurgence is not operative in QFTs on general grounds. This changed recently due to two parallel but distinct ideologies. On the one hand, Unsal and Argyres [24,25] conjectured that renormalon singularities have a semiclassical explanation if the problem is approached from the regime of weakly coupled theory via the idea of adiabatic continuity [26][27][28]. Indeed in such regimes it was shown that renormalon singularities disappear [29], and resurgence is likely operative. However this is difficult to test as no access to high orders of perturbation theory is typically available in QFTs. Nevertheless certain 1+1D models, when dimensionally reduced to quantum mechanics via the special kind of compactification, has weakstrong coupling adiabaticity and resurgent structure [30][31][32]. Resurgence is likewise useful in quantum field theories without renormalon singularities, for instance the Chern-Simons theory [2,33] and certain supersymmetric field theories. Relatedly resurgence also finds its use in topological string theories, where Borel resummation and resurgence techniques have been used to explore non-perturbative contributions and to turn the asymptotic series of topological string free energy into a finite function [34][35][36][37][38][39][40][41][42][43], culminating in [44][45][46][47][48].
Since resurgence is tightly connected with high orders of perturbation theory, it is of immense practical use to have an efficient way to computer high orders of perturbation theory. Recently in [49] a Mathematica package called BenderWu was developed using the method originally used by C. M. Bender and T. T. Wu [50] for an anharmonic oscillator, which efficiently computes symbolic perturbative solutions to a generic one dimensional quantum mechanical problem with the Hamiltonian of the form (1.1) a second order differential operator, where V (x) is an arbitrary non-singular potential, around one of its harmonic minima.
Many quantum mechanical problems also exist whose Hamiltonians are difference operators. They can be regarded as the relativistic version of ordinary quantum mechanical systems, for instance, the relativistic Toda lattices [51], the elliptic Ruijnaars-Schneider systems [52,53], the cluster integral systems [54], and etc. A particular type of relativistic quantum mechanical systems that has recently attracted a lot of attention is quantum mirror curves, and their studies have been extremely fruitful. Consider topological string theory whose target space is a toric Calabi-Yau threefold. The mirror curve to the threefold is the moduli space of the branes compatible with the toric structure [55]. The quantisation of the mirror curve gives rise to Hamiltonian operators of the type H(x, p) = (r,s)∈I a r,s e rx+sp , a r,s ∈ R , (1.2) where I is a finite set of integer pairs, and x, p satisfy the canonical commutation relation [x, p] = i . The wave-functions to these Hamiltonians are related to the open topological string partition function associated to the branes [56] 1 . It is later understood that the quantum mirror curve is more closely related to the refined topological string in the Nekrasov-Shatashvili limit [58]. The quantum mirror curve defines a spectral problem, whose quantum-corrected WKB periods coincide with the quantum deformation of the periods of the Calabi-Yau, while the latter determine the NS topological string free energy F NS via the so-called quantum special geometry relation [59][60][61].
The exact solution to the spectral problem, however, remained elusive until [62]. Naively one would conjecture that the spectral problem is solved by the Sommerfeldtype quantisation condition where k i are the levels of the eigenenergies, and a = (a i ) are the quantum periods. The equation (1.3), nevertheless, cannot be the full story, as the l.h.s., which can be understood as the quantum phase space, have poles whenever is 2π multiplied by a rational number. Important non-perturbative corrections were first found in [63] to cancel the poles, which, after the numerical work [64] that reveals more subtle corrections are needed, led to the exact spectral theory for quantum mirror curves [62,65], followed by a detailed study of wave-functions [66,67], especially in the special case when = 2π (see related works [68][69][70]). One amazing feature of the spectral theory is that it also defines conjecturally a non-perturbative completion of topological string free energy in the conifold frame, which coincides with the results of resurgence analysis [48]. This conjecture was proved in a special example in certain limit in [71]. See review [72] and related works [73][74][75][76][77][78][79]. Furthermore, it has recently become clear that the quantum mirror curve is the quantum Baxter equation of the cluster integrable system [54] associated to the toric Calabi-Yau threefold. Inspired by an elegant reformulation [80] of the quantisation condition in [62], a conjectural exact quantisation condition for the cluster integrable system is also written down [81,82]. The interplay between the quantisation conditions for quantum mirror curve and those for cluster integrable system led to an interesting set of relations for BPS invariants of the Calabi-Yau [83], and they were proved in a special category of examples in [84].
To study these systems, we will generalise the algorithm presented in [49] to difference equations of type (1.2) and study their spectrum. We added to the mathematica package BenderWu 2 of [49] a function called BWDifference, which computes efficiently perturbative solutions to one dimensional quantum mechanical problems whose Hamiltonian is a difference operator of the exponential-polynomial type given in (1.2). This allows us to study the spectral problem of quantum mirror curve perturbativelly to a very high order (≥ 100) in .
When the toric Calabi-Yau threefold is fano, the Hamiltonian operator arising from the quantisation of mirror curve is unique. Y. Hatsuda [85] argued that in the case of one particular toric fano Calabi-Yau threefold, the local F 0 , the perturbative eigenenergies of the Hamiltonian operator are Borel summable and that the Borel sums of the perturbative eigenenergies agree well with the numerical spectrum. The study in [85] was up to 36 orders in . With the BWDifference function we are able to extend the study of the local F 0 to 100 orders in , and confirm that the Bore-Pade partial sums continue to converge to the exact (numerical) result.
Furthermore we study the perturbative solutions to the Hamiltonian operator associated to all toric fano Calabi-Yau threefolds using the function BWDifference in the BenderWu package, and find strong evidence that the spectrum of all of them is Borel summable and that the Borel sum gives the correct answer.
The paper is organized as follows. In the next section we describe the adapted Bender-Wu algorithm that solves perturbatively the Hamiltonian difference operators, and how to use the Mathematica function that implements the algorithm. In Sec. 3, we explain the Hamiltonian operators arising from the quantisation of mirror curve in topological string theory on a toric Calabi-Yau threefold, especially when the Calabi-Yau is fano, before proceeding to provide evidence that the perturbative eigenenergies of Hamiltonians associated to all toric fano Calabi-Yau threefolds are Borel summable. Finally in Sec. 4 we conclude and discuss possible future directions. We relegate to the Appendix the derivation of the adapted Bender-Wu algorithm, as well as the explanation of the technical observation that all the Hamiltonians we have considered have a unique classical minimum.
The Hamiltonian operator is an self-adjoint operator over the domain D which consists of wave-functions Ψ(x) that are not only themselves L 2 (R) integrable but that e x ψ and e p ψ are also L 2 (R) integrable. This constraint can be translated to the condition in the coordinate representation (see for instance [86]) that the wave-function Ψ(x) admits an analytic continuation into the strip where it is L 2 (R) along the x-axis for any fixed value of y, and that the limit exists.
To make the analysisà la Bender-Wu, it is convenient to rescale x = √ x, p = √ p. This scaling would not change the eigenvalue E nor the eigenfunction Ψ, provided that x,p satisfy the commutation relation [x,p] = i. In the coordinate representation,p is the differential operator −i∂ x . The Hamiltonian operator now reads H √ x, Let us further assume that the Hamiltonian as a function has a local minimum at the origin; in other words, H in small expansion has no linear term inx orp. If this is not the case we can always use a canonical transformation which takes (x, p) → (x + x 0 , p + p 0 ) to achieve this, which amounts to the redefinition of a r,s 3 . Now let us expand the operator in powers ofx andp. Up to an overall constant, we get We wish to solve this equation perturbatively in the expansion of small √ . We show in the Appendix that the energy E and the wave-function Ψ(x) have the following expansion and where ψ k (x) is the level k normalized wave-function of a harmonic oscillator with unit mass and frequency. The prefactor e iαx 2 /2 of wave-function expansion comes from another canonical transformation that makes the second term in the small expansion of H( √ x, √ p) (2.6) into the Hamiltonian of a harmonic oscillator. In the Appendix we give the detailed derivation of an algorithm that solves recursively the expansion coefficients E l−2 ,Ã k l . To summarise, we find that in the lowest orders, is the classical energy, E −1 = 0, and where the non-negative integer ν specifies the level of the eigenenergy. Fixing the level ν, one finds in the lowest orders for the wave-functioñ where settingÃ ν 0 to unity is a normalization choice. Furthermore, we can normalize the wave-function so thatÃ ν l = 0 , l ≥ 1 . (2.14) To obtain higher order solutions, we first define and β(r, s) = (r − CB −1 s)ξ + isξ −1 .
(2. 16) Then assuming all the coefficientsÃ k l and l are known for l < l, the coefficientsÃ k l and l can be computed from the following recursive relations respectively, +β n F (−k, −q; 1 + n; 2)Ã k+n l+2−n−2q From the recursion relation (2.17) and the initial condition (2.13) one also finds that A k l = 0 whenever k > 3l + ν.
We have in fact programmed a function called BWDifference for Mathematica which computes the expansion coefficientsÃ k l , l automatically and added it to the updated BenderWu package [49]. Before we proceed to explain how the function can be used, we would like to make three claims here about the structure of the perturbative eigenenergies and wave-functions: (i) There is a unique perturbative solution (up to the normalization constant) of the form (2.9) for any given level number.
(ii) Energy expansion contains only powers of , not powers of √ .
(iii) The perturbative wave-function can always be constructed to obey 19) to every order in perturbation theory.
To prove claim (i), consider the difference equation of the form (2.8). Let us show that this equation cannot have two solutions with the same eigenvalue, both of which reduce to harmonic oscillator solutions as → 0. Indeed if this were the case, the two solutions must be orthogonal to each other. But this would mean that in the → 0 limit, the two solutions reduce to orthogonal harmonic oscillator solutions with different eigenenergies. This violates the assumption that they have the same eigenvalue. Hence we conclude that only one such solution exists. We can also see that this is the case from the recursion equations (2.17), (2.18), as choosing the coefficients 4 A ν l uniquely fixes the solution. Now let us go back to (2.8) and prove the claim (ii). One easy way to see this is to notice that the eigenenergies of H(x, p) have an expansion in not in By complex conjugating the second equation, and subtracting from the first we get that either Ψ |Ψ = 0 or E( √ ) = E(− √ ). However we also know that Ψ(−x, − √ ) and Ψ(x, √ ) cannot be orthogonal, because they reduce to the same harmonic-oscillator solution in the → 0 limit. Hence we must have E( i.e. energy must be an even function in √ , which means that the eigenvalue series expansion is in even powers of √ only. Claim (iii) immediately follows from the above. Since Ψ ν (x, ) and Ψ ν (−x, − ) are wave-functions of the same eigenenergy of level ν, we can construct a new wavefunction of again the same eigenenergy, and it satisfies the condition (2.19). This parity condition implies that which is compatible with the initial condition (2.13) that we choose. From the point of view of the recursion calculation, if the above condition onÃ k l is satisfied for all l <l, then by virtue of (2.17) we have thatÃ k l for (−1)l +k+ν = −1 is given entirely by coefficients which vanish, and hence they vanish themselves.
Incidentally, from (2.18) we can see that if l is odd, the r.h.s. contains coefficients which all vanish by (2.23), confirming the claim that only even powers of √ appear in the expansion of E.

How to use the BWDifference function
Here we present the BWDifference function which is incorporated into the updated BenderWu [49] package of Mathematica. This function solves perturbatively the difference equation of the form were ν is the level number and H(X, P ) is the "Hamiltonian" which depends on the momentum and coordinate displacement operators for integer r and s (note that these can be negative as well). The . . . indicates an ordering of X and P . A conventional ordering which renders the operator H(X, P ) Hermitian is given by This ordering is assumed by the BWDifference function. Furthermore the BWDifference function assumes that at X = P = 1 (i.e. x = p = 0) the classical function H(X, P ) attains (at least a local) minimum.
The BWDifference function produces a perturbative expansion of "energy" E at level ν and an unnormalized wave-function Ψ(x), of the form given in (2.9). As we have shown in the previous section, the energy is always in powers of , not √ . This means that all E l in equation (2.9) vanish whenever the l is odd. For this reason the code returns only even coefficients of E, i.e. returns E 2n . From now on when we talk about the "order" of the perturbative expansion we will mean the number n, rather than the order of √ , for which we reserve the letter l. Now note that n = −1 is the leading order (i.e. classical energy) which is identical to E −2 = H(X = 1, P = 1) and is of order 1/ in our convention.
In order to access the BWDifference function, one must first install the BenderWu package bundled with this work. Alternatively the most up-to-date version can be downloaded at http://library.wolfram.com/infocenter/MathSource/9479/ After following the installation instructions, the package must be loaded via the command

<<"BenderWu'"
This allows the user to access all the functions in the BenderWu package, in particular the BWDifference function relevant for this work. Now let us see how the BWDifference function works. It takes in four essential arguments: the form of the Hamiltonian H(X, P ), the name of the two variables X, P as a list of two elements, i.e. {X,P}, the level ν, and the order l max to which the energy E ν shall be computed. The typical syntax is given by which computes the perturbative expansion of the second level, to the 5 th order in .
Once the computation is done, the function returns a list with three elements. The first element is the list of coefficients {E −2 , E 0 , E 2 , E 4 , E 6 , E 8 , E 10 }, while the second is a matrix of coefficients A k l where the l-index denotes the rows and the k-index the columns. The third element is not important for the user, and only serves for proper functioning of the function BWProcess, which was introduced in [49]. Hence if we execute the command we will get a list of perturbation series coefficients l , which in this case is Alternatively one can use an option Output->"Energy" instead, i.e. with the same outcome as before. However a better way to use the code is to assign the output to a variable, and use the function BWProcess introduced already in the original BenderWu package [49] to control the output without having to recompute the expansion. In other words the benefit of using the BWProcess function is that one can make a computation to a high order once, and use the BWProcessto analyze the result without having to recompute the expansion. For example, if we call the line it assigns the output of BWDifference to a variable BW, and hence contains all the perturbative information to order 20 in . In order to output the energy coefficients, we can simply call

BWProcess[BW]
which produces the output where the dots stand for the terms not written. Often the computation will involve many terms, and the output can be quite bulky. Therefore the BWProcess function has an option which allows the user to display only the limited order, for example

BWProcess[BW,Order->5]
Furthermore the BWProcess function can be used to specify the lower and upper bounds of the perturbative order, as in which gives an output To obtain the wave-function coefficients, all we need to do is to use the option Output->"WaveFunction". For instance, calling the line The first row is A k 0 , which is zero except for k = ν = 2. The second row is A k 1 , the third row A k 1 , etc. To get the element A 5 3 , all we need is take the (4, 6) element of this output, i.e. calling BWProcess[BW, Output -> "WaveFunction"][ [4,6]] which returns 20 9 .
We can also use an option OutputStyle->"Series" to output the series (2.9) for the wave-function. For example writing BWProcess[BW, Output -> "WaveFunction",OutputStyle->"Series",Order->1] produces the following output where g is 5 √ , so g 3 is of order 3 . Note that the prefactor of (2.9) is not included. To include it use the option Prefactor->True Let us define the wave-function and energy to the 10th order of with the commands The difference equation for the difference operator H = X + P + 1/(XP ) , with X = e igx , P = e igp , explicitly reads  21 ], so that the equation is satisfied at least to the 20th order in g = √ . Finally we discuss briefly the option Imaginary. The solution of the difference equation ψ(x) need not be real (up to a constant phase), and the coefficientsÃ k l can have imaginary parts. The example we studied so far returns purely real coefficients A k l (see Appendix A). When the coefficients are not real, the algorithm may slow down significantly, especially if large orders need to be computed. In order to improve this, a refined algorithm is built into the BWDifference function which speeds up the computation when the coefficients are complex by separating the real and the imaginary parts of the coefficients. To switch to the refined algorithm, one needs only to add Imaginary->True in the option list of the BWDifference function. For concrete examples, see the example notebook included in the BenderWu package.

Application: quantum mirror curves
We describe here the Hamiltonian operators arising from the quantisation of mirror curves in topological string theory on toric fano Calabi-Yau threefolds, and then apply our Bender-Wu algorithm to solve perturbatively the eigenvalue problem of the Hamiltonians.

Quantum mirror curves
Consider topological string theory on a toric Calabi-Yau threefold [87][88][89][90]. A toric Calabi-Yau threefold X Σ can be succinctly described by its toric fan Σ. The toric fan consists of n Σ + 3 1-cones and the triangulation of the convex hull of the 1-cones. The 1-cones are subject to n Σ linear relations The Calabi-Yau condition demands that one can always rotate the toric fan so that the endpoints of the 1-cones have coordinates It is therefore enough to present the toric fan by the image of the projection onto the plane (1, •, •), a triangulated convex integral polygon whose vertices are v α = (r α , s α ) , α = 1, . . . , n Σ + 3 .
which preserves the linear relation vectors Mirror symmetry dictates that the free energies of topological string theory on the Calabi-Yau threefold X Σ can be computed from the mirror curve C Σ , a noncompact Riemann surface, whose Newton polygon coincides with the fan support of X Σ . Therefore given the fan support N Σ with vertices v α , the equation of C Σ reads n Σ +3 α=1 a α e rαx+sαy = 0 , x, y ∈ C . (3.6) The coefficients a α in the equation (3.6) parametrise the complex structure moduli space of the mirror curve. They are not all independent, as three of them can be scaled to one through the C * scalings on e x , e y and an overall scaling. It is customary to set to 1 three coefficients associated to vertices on the boundary; the number of internal vertices gives the genus g Σ of the mirror curve. Due to physics consideration, the g Σ coefficients associated to internal vertices are called the true moduli, while the remaining coefficients associated to boundary vertices after fixing the (C * ) 3 scaling are called mass parameters 6 .
In this paper for simplicity we restrict ourselves to fano Calabi-Yau threefolds whose fan supports are reflexive, in other words convex Newton polygons with only one internal vertex. Reflexive 2d polygons have been classified up to the SL(2, Z) isometry, and they are listed in Fig. 3.1 (see for instance the construction in [91,92]). Since they have a single internal vertex, and it allows for a canonical way of writing down the curve equation by putting the only internal vertex at the origin. For instance, the canonical equation for the first polygon in Fig. 3.1 is while the second polygon in Fig. 3.1 gives In these equations u is the true modulus of the model. Note the canonical form still enjoys the SL(2, Z) isometry acting on the exponents To quantise the mirror curve, we simply promote the coordinates x, y to quantum operators x, p satisfying the canonical commutation relation [x, p] = i through the Weyl quantisation prescription e r i x+s i y → e r i x+s i p . (3.10) Here is assumed to be real. For a genus g Σ mirror curve, one can in principle construct g Σ mutually non-commutative Hamiltonian operators, each associated to a different true modulus [65]. The mirror curve of a fano Calabi-Yau threefold is always of genus one, and thus the associated Hamiltonian is unique. It is obtained by taking the l.h.s. of the canonical equation of curve, removing the true modulus u, and then performing the quantisation procedure. In the example of (3.7), we get The SL(2, Z) isometry of the Newton polygon then corresponds to canonical transformations on x, p.
In this paper, we are interested in the eigenvalue problem of the Hamiltonian operator associated to a toric fano Calabi-Yau threefold, in the following form 7 where k is the level number. In [62] a conjectural quantisation condition was given using the (refined) topological string free energies to solve exactly the spectrum of H(x, p). In this paper, we are interested in the perturbative solution to the Hamiltonian eigenvalue problem, and we will not need the input of topological string. Clearly the Hamiltonian operator is of the form (2.1), and so its eigenvalue problem can be treated by our BWDifference function. We also call the polynomial of e x , e y before quantisation the Hamiltonian function H(x, y), and it is the analogue of the classical potential in a nonrelativistic quantum mechanical problem. Consider the perturbative expansion of E (ν) in terms of which is an asymptotic series with zero radius of convergence. Hatsuda in [85] gave evidence that for the second geometry in the list of Fig. 3.1 with the mass parameter set to 1, the Borel sum of the perturbative eigenenergies for finite values of agrees with the numerical results, implying the Borel summability of the eigenenergy series. We want to expand the exploration in [85] to other reflexive geometries with higher precision. The precision of Borel resummation depends crucially on the order of asymptotic series that is included. [85] fixed the coefficients of the perturbative eigenenergies by comparing the asymptotic series with numerical eigenenergies computed by numerous small values of , and in this way, [85] could only obtain up to order 36 of the perturbative eigenenergies for the said geometry. Our BWDifference function provides a far more efficient way to compute perturbative eigenenergies. For instance, for the same geometry the BWDifference function can easily compute the eigenenergy series at level 0 up to order 100 within 240 seconds on an ordinary desktop computer. This results in an agreement between the Borel sums with the numerical results for = π up to more than 25 digits, compared to only 12 matching digits in [85]. We analysed all sixteen reflexive Newton polygons listed in Fig. 3.1, corresponding to all possible toric fano Calabi-Yau threefolds, for appropriately chosen values of mass parameters. We find that for each model the poles of the Borel transforms of the perturbative eigenenergies are never located on the positive real axis of the Borel plane, indicating Borel summability. Besides, the Borel sums of the eigenenergies have very good agreement with the numerical results, and the degree of agreement increases consistently when more orders of perturbative series are used in resummation. We therefore confirm and expand to all toric fano Calabi-Yau threefolds the observation in [85] that the Borel-Padé resummation captures the exact eigen-energies. The details of the results are discussed in the next section.

Results
We first write down in Tab. 3.1 the Hamiltonian operators for each of the 16 reflexive Newton polygons listed in Fig. 3.1. It is beneficial if we can rearrange the Hamiltonian operator so that it is invariant under the reflection p → −p. We call such an operator p-parity even. From the point of view of perturbative calculation via the BWDifference function, the wave-functions of a p-parity odd Hamiltonian operator are complex, and the computation is significantly slowed down compared to the cases where wave-functions are real. This problem can be circumvented by turning on the option Imaginary->True in the BWDifference function, which then separates the real and the imaginary parts of complex wave-functions explicitly to cure the slowdown. From the point of view of numerical calculation, when working in the coordinate representation, the p operator is −i ∂/∂ x . As a consequence, if the Hamiltonian operator is p-parity even, the Hamiltonian matrix with entries n|H|m would be real symmetric instead of complex Hermitian, and thus the matrix diagonalisation would be faster.
Here is an appropriate place to recall the method of numerical calculation of spectrum (see for example [64]). We choose the basis of wave-functions in the domain of H to consist of the eigenfunctions of the quantum harmonic oscillator with both mass and frequency set to 1, i.e.
Here H n (x) are Hermite polynomials, and they obey the following orthogonality conditionŝ where L α n (z) are Laguerre polynomials. Then for the operator e rx+sp , we have n 1 |e rx+sp |n 2 = n 1 !n 2 ! e |z| 2 2 z n 1z n 2 min(n 1 ,n 2 ) Clearly the Hamiltonian matrix n 1 |H|n 2 is real and symmetric if and only if every monomial e rx+sp is paired with e rx−sp , in other words, the Hamiltonian operator is p-parity even. Among the 16 reflexive Newton polygons, the Hamiltonians of all but four geometries, namely F 6 , F 8 , F 10 , F 11 , can be put via a canonical transformation to a form that is p-parity even for appropriately chosen values of mass parameters. This is the form of the Hamiltonians presented in Tab. 3.1.
When mass parameters are non-negative, the Hamiltonian functions for the operators in Tab. 3.1 have a unique minimum, as is shown in the Appendix B 8 for real values of x, y, which is taken to be the classical ground state. The uniqueness of the classical ground state also indicates the absence of real instantons, and could be related to the Borel summability of the spectrum that we find here.
From the point of view of perturbative solutions, the BWDifference function expands around a minimum of the Hamiltonian function which it assumes to be (x, y) = (0, 0). Therefore when the actual minimum (x, y) = (x 0 , y 0 ) is not at the origin, we have to shift the coordinates x, y by hand before feeding the Hamiltonian function into the BWDifference function. Furthermore, the BWDifference function runs much faster if the Hamiltonian function after the shift of coordinates has no irrational coefficients. We can always achieve this by taking appropriate values of mass parameters.
As we have seen, in order to most efficiently use the BWDifference function, we would like to choose rational values of mass parameters such that geometry mass parameters geometry mass parameters   • the Hamiltonian operator is p-parity even (not applicable to F 6 , F 8 , F 10 , F 11 ); • the coordinates (e x 0 , e y 0 ) of the minimum of the Hamiltonian function are rational numbers.
We choose one set of mass parameters for each geometry satisfying these conditions, and list them in Tab. 3.2 (F 1 has no mass parameter). Let us focus for the moment on the polygon F 2 , which represents the Calabi-Yau threefold called the canonical bundle over the Hirzebruch surface F 0 or local F 0 , and we set the mass parameter m 1 = 1, as indicated in Tab. 3.2. As already mentioned in Section 3.1, we can compute the perturbative series of the ground state energy up to order 100 with relative ease. Now given the asymptotic series E (ν) ( ), we can compute the Borel transform which is a convergent series. The Borel transform may have poles in the ζ-plane, also known as the Borel plane, and the locations of the poles are the actions of the instantons of the relevant quantum mechanical system. If no pole lies on the positive real axis, we can perform the Laplace transformation on the Borel transform n is expected to be dictated by the saddles of the phase-space functional associated with the partition function of the difference operator. Indeed preliminary studies of the model of local F 0 indicate that E (0) n (−1) n ∼ n!/(2|S|), where S is the action of a complex instanton tunneling from the minimum at x = p = 0 to one of the closest complex minima (say p = 0, x = 2πi). On the other hand, generic cases are complicated by the fact that the instanton actions are complex (in the local F 0 model the leading instanton action is real and negative). We leave detailed studies of this kind for the future.
We proceed to compute the Borel-Padé sums of the perturbative ground state energy, evaluate them at = π, 2π, and 11π/7, and compare with the numerical results. As seen in Tab. 3.3, both sides agree extremely well: the column of = π agrees to 26 identical digits when 100 orders of are taken. To better illustrate the success of the Borel-Padé resummation, we define the matching degree between two order = π = 2π = 11π/7 which roughly speaking gives the number of identical digits between the two. We plot in Fig. 3.3b the matching degree between the Borel-Padé sum and the numerical result against the truncation order of the perturbative series. It is very satisfactory to see that the matching degree grows up consistently with the perturbation order up to a very high value. We perform the same analysis for the other 15 models listed in Tab. 3.2. We find that in all 15 models, there are no stable poles along the positive real axis in the Borel plane, and we find agreement between the Borel-Padé sums of the perturbative ground state energy and the numerical results, the degree of which improves consistently with increasing truncation order of the perturbative series. The plots of matching degrees for all 15 models are given in Figs. 3.3, 3.4. Finally, we give in Tab. 3.4 for all models the digits of the Borel-Padé sums which are both stabilised and identical with the numerical results.
We mention in passing that the underlying reason for the Borel summability of the spectrum is likely a consequence of the fact that no real-positive action instanton solutions exist in the limit of → 0 9 . To show this one would need to carefully study the stokes phenomena as the phase of is varied. We leave it as an open problem for the future.

Conclusions and future prospects
In this paper, we are interested in solutions to the eigenvalue problem of Hamiltonian operators which are difference operators, and in particular polynomials of e x , e p with the commutation relation [x, p] = i . We developed Bender-Wu like recursion relations that solve eigenenergies and wave-functions perturbatively in small , and implemented the algorithm for Mathematica in a function called BWDifference in the updated BenderWu package, originally developed in [49]. Our algorithm is very efficient, capable of computing more than one hundred orders of perturbative solutions for a typical Hamiltonian difference operator in a reasonable amount of time.
Typical Hamiltonian difference operators appear in the quantisation of mirror curves in topological string theory on toric Calabi-Yau threefolds. We studied all sixteen toric fano Calabi-Yau threefolds, whose associated Hamiltonian operators are unique, and computed the perturbative ground state energies for some choice of mass parameters. We find strong evidence that the perturbative eigenenergies are Borel summable, and the Borel sums are exact. Although we only studied and presented explicitly ground states, one can easily check the situation is the same in excited states. A possible reason for Borel summability is that the Hamiltonian difference operators arising in mirror curve quantization all have a unique real minimum as a function of x and p, so that that classical equations of motion do not allow for real-positiveaction instantons. However Borel summability (or even convergence) does not always mean that the re-summation gives the correct result, and non-perturbative corrections may still arise (see e.g. [9,10,93]). Nevertheless the perturbation theory is factorially growing, and is likely dictated by the complex instanton (or ghost instanton) solutions which are generically present in such systems. Our BWDifference function can be used to address such features in detail.
In addition, the mirror curves associated to toric fano Calabi-Yau threefolds are all genus-1 curves. In quantum mechanical systems with genus-1 curves obey a remarkable relation -theÁlvarez-Casares relation ( [12][13][14][15], more examples are later found in [9,10,16,[94][95][96][97]) -between the trivial perturbation theory and perturbation theory around instantons. Similar relation should exist for the ghost instantons of the quantum mirror curves. However in this case, as we argued, such nonperturbative objects do not contribute in the trans-series expansion, but rather dictate the   asymptotic growth. Because of this, we expect a form of self-resurgence [9] to hold, if the analogousÁlvarez-Casares relation holds in quantum mirror curves.
In this paper, we restrict ourselves to quantum mirror curves of fano Calabi-Yau threefolds. It would also be interesting to look at more generic toric Calabi-Yau threefolds, whose mirror curves have genera g great than one. There would be two different but related quantum mechanical problems. The first is again the quantisation of mirror curves. The associated Hamiltonian operators are no longer unique [65], but each of them still gives rise to a one-dimensional quantum mechanical system. Alternatively, one could also look at the cluster integral systems [54] associated to the toric Calabi-Yau threefolds. When the latter are not fano, the cluster integral systems are higher dimensional, involving g mutually commutative Hamiltonian operators, all of which are difference operators of the exponential-polynomial type. In order to study perturbatively the cluster integral systems, we would need to generalise our algorithm to treat multivariable systems. In addition, it would be desirable to further generalise the Bender-Wu algorithm for more generic difference operators, not necessarily of the exponential-difference type. Another interesting question is whether our Bender-Wu solutions to the Hamiltonian difference operators, assuming wave-functions can be expanded in terms of the wave-functions of harmonic oscillators, exhaust all possible wave-functions in the domain D of the Hamiltonian difference operators.

A The Bender-Wu recursion relations
We derive in detail here the solution to the following eigenvalue problem with Bender-Wu type recursion relations, so that to lowest orders the Hamiltonian operator becomes a simple harmonic oscillator with unit mass and frequency const. + 2 Hence we can define a reduced Hamiltonian Now we wish to solve the eigenvalue equation 10 h √ x, The eigenvalues and wave-functions are related to those of H √ x, √ p by We also comment here that although h( √ x, √ p) is a difference operator, at any finite order in √ expansion it is a polynomial inx,p and thus a finite order differential operator.
In order to solve the eigenvalue problem (A.10), it is convenient to write the coordinatex and momentump operators in terms of the creation and the annihilation operatorsx r,s β n 1β n 2 |β| 2n 3 n 1 !n 2 !n 3 ! g n 1 +n 2 +2n 3 −2 1 2 n 3 (a † ) n 1 a n 2 (A. 17) where n 1 , n 2 , n 3 run from 0 to infinity. Now we make a formal expansion of the wave-function and eigenvalue where ψ k (x) are eigenfunctions of the simple harmonic oscillator with unit mass frequency (i.e. ψ ν (x) are solutions of the leading order spectral problem). Using the fact that 11 (a † ) n 1 a n 2 ψ l = k!(k + n 1 − n 2 )! (k − n 2 )! ψ k−n 2 +n 1 , (A. 19) we get n 1 ,n 2 ,n 3 r,mã r,s β n 1β n 2 |β| 2n 3 n 1 !n 2 !n 3 ! (A.20) By equationg powers of g and coefficients of ψ k on both sides, we have n 1 ,n 2 ,n 3 r,mã r,s β n 1β n 2 |β| 2n 3 n 1 !n 2 !n 3 !2 n 3 Notice that we can formally assume that n 1,2,3 run from −∞ to +∞, noting that the factorials have poles at negative integer values, and that A k l vanishes for negative k or l. Then we can freely shift n 1 → n 1 + n 2 without worrying about the limits of the sum, and get n 1 ,n 2 ,n 3 r,mã r,s β n 1 |β| 2(n 2 +n 3 ) (n 1 + n 2 )!n 2 !n 3 !
where F (a, b; c; z) = 2 F 1 (a, b; c; z) is the hypergeometric function. Relabeling n 3 by q, n 1 by n or −n if n 1 is negative, we have that n≥0,q r,sã r,s β n |β| 2q n!q! (A.25) Finally we make the shift q → q + n in the second summation to make the expression above in a nicer form. where on the l.h.s. the n = 0 term has been singled out. This identity is valid for any k, l ≥ 0. 12 This simply follows from the definition of the hypergeometric function where (a) s = (a)(a + 1) . . . (a + s − 1) = Γ(a+s) Γ(a) . If a, b are negative integers −n, −m then we can use the Gamma-function reflection formula to get that (−n) s = (−1) s n! (n−s)! . Further if we take that c = q + 1 with q ∈ N 0 , we have that (q + 1) s = (q+s)! q! so that Next we consider (A.26) when l ≥ 2. Note that the the summand of the first summation when q = 0 always cancels with the term proportional to −2 on the r.h.s., and that the summand of the second summation when (n, q) = (1, 0), (2, 0) vanish due to the identities (A.15). Therefore (A.26) becomes Here we have shifted the index l → l + 2 on both sides, n → n + 2 on the r.h.s., and then singled out the terms proportional to A k l . Now notice that the sums on both the left and right hand side contain only coefficients A k l withl < l. So by inserting l = 0 all that remains is (2k Since not all A k 0 vanish, this identity can only be true if for some nonnegative integer ν 0 = 2ν + 1, A ν 0 = γ = 0 and A k 0 = 0 , ∀k = ν , (A. 32) where γ is an arbitrary nonvanishing constant. ν serves as the level of the eigenvalue/wavefunction solution. Fixing the level ν, we can normalise the wave-function so that A ν l = 0 , ∀l > 0 . (A.33) Following the normalisation of wave-function above, (A.30) gives us two recursion relations that solve A k l and l respectively. Assuming that A k l and l are known for all l < l, the expansion coefficients A k l and l can be solved from However, another source of irrational factors can be ξ which appears in the definition of β. However from (A.13), we can see that by definingβ to be β = ξβ and appropriately rescaling of the coupling, the difference equations can be converted to involve only ξ 2 in the imaginary part ofβ.
A little thought reveals that the difference equation can be setup in such a way that the real part ofÃ-coefficients depend only on the even power of ξ 2 , while the imaginary part always contains an odd power of ξ 2 . This means that the real part contains no irrational factors, for a choice of rational choice of all a r,s , and, if ξ 2 is irrational, the imaginary part ofÃ will be always proportional to ξ 2 , multiplying a rational number.
By splitting the difference equation forÃ into its real and imaginary parts, the irrational coefficients appear in a predictable manner, and such treatment of the difference equations whenever the imaginary part ofÃ is non-vanishing, speeds up the Mathematica algorithm of the package significantly. To activate this feature one needs to call the option Imaginary->True in the BWDifference function, as described in the text.

B Proof of uniqueness of minima
We prove here that the Hamiltonian of the form (1.2) with positive a r,s ≥ 0 has a unique real minimum or no minimum as a function of x and p. Equivalently we can also prove that as a function of X = e √ x and P = e √ p there is a unique minimum such that P > 0, X > 0.
The proof goes as follows. The minimum satisfies the equation Since B r > 0, we have two options. If r min ≥ 0 then the above polynomial has only positive coefficients. If r min < 0 then the coefficients of X k with k < −r min have a negative sign, while the rest are positive. We now invoke the rule of Decartes which says that the number n p of positive real roots of a polynomial is less then or equal to the number of the monomial sign variations of the coefficient. Since this sign variation is 0 or 1, we must have at most one solution.
The same argument can be invoked to show that the equation ∂ P H = 0 has only one solution in P for any X > 0. This concludes our proof.