Spectra of elliptic potentials and supersymmetric gauge theories

We study a relation between asymptotic spectra of the quantum mechanics problem with a four component elliptic function potential, the Darboux-Treibich-Verdier (DTV) potential, and the Omega background deformed N=2 supersymmetric SU(2) QCD models with four massive flavors in the Nekrasov-Shatashvili limit. The weak coupling spectral solution of the DTV potential is related to the instanton partition function of supersymmetric QCD with surface operator. There are two strong coupling spectral solutions of the DTV potential, they are related to strong coupling expansions of gauge theory prepotential at the magnetic and dyonic points in the moduli space. A set of duality transformations relate the two strong coupling expansions for spectral solution, and for gauge theory prepotential.


Introduction
In this paper we report results about a relation between N=2 SU (2) supersymmetric QCD with N f =4 hypermultiplets and asymptotic spectral solutions of Schrödinger equation with the Darboux-Treibich-Verdier (DTV) potential [1,2,3,4]. This is an example of the connection between Omega deformed N=2 supersymmetric gauge theories and quantum integrable models discovered by Nekrasov and Shatashvili [5]. The theoretical models discussed here are simple enough to dig up many details, meanwhile comprise of various nontrivial ingredients to demonstrate some aspects of this connection.
The Schrödinger equation with DTV potential is the elliptic form of Heun differential equation [6,7]. The Heun equation arises in a variety of physics problems, usually appears in different forms. There are some interest in obtaining the asymptotic spectral solution of this differential equation, it is difficult to obtain the solutions presented here by conventional perturbation method. In this paper we clarify the structure of asymptotic solutions of a particular kind, and carry out perturbative computations, by methods stem from its connection to the N f =4 super QCD model.
The N f =4 super QCD is one of the examples used by Seiberg and Witten to demonstrate the low energy physics of N=2 supersymmetric gauge theory [8,9]. We show in this paper that the spectral solutions of the DTV potential precisely produce the gauge theory prepotential expressed by elliptic modular functions, expanded both in the weak coupling region and in the strong coupling region.
The material in this paper is based on partial results obtained in Ref. [10], here we give a complete and consistent analysis for this problem. The spectral problem and the N f =4 super QCD model are more general than the examples we have studied before, the results obtained here confirm some facts and methods observed in other examples, provide more persuasive argument for some puzzling points encountered earlier, including the monodromy relations of strong coupling solutions [11,12,13].
The paper is organised as follows. In Section 2 we explain the structure of asymptotic solutions of the DTV potential which have a connection to gauge theory, this generalise a fact that was verified for a few simpler elliptic potentials in previous papers. In Section 3 we compute the weak coupling perturbative solution for the DTV potential, it is related to the low energy effective theory of super QCD model with surface operator under partial Ω deformation. In Section 4 we compute the strong coupling perturbative spectra of the DTV potential, discuss the role of two elliptic Hill potentials with infinitely many terms, given in (57) and (88). The two elliptic potentials, and their spectral solutions, are dual to each other in the same sense of monopole-dyon duality in gauge theory. In Section 5 we compute the strong coupling expansions of prepotential for massive N f =4 super QCD at the magnetic and dyonic singularities in the moduli space. In Section 6 we discuss various choices of coupling constants of the DTV potential by which the DTV potential reduces to the Lamé potential using Landen transformations of elliptic functions, correspondingly prepotential of the N f =4 super QCD model is connected to prepotential of the N=2 * super QCD model by similar reductions.

Notice on notations
This paper develops some results presented in a few earlier papers [10,11,12,13]. Various forms of elliptic functions are frequently encountered, in earlier papers we have used different notations for the elliptic modulus/nome according to content discussed there. To avoid confusion, we clarify the convention of notations used in each paper.
Refs modulus/nome gauge theory Jacobi sn/ cn/ dn Weierstrass ℘ E and K this paper q q p x 2 [10] q q p k 2 [11][12] [13] q k 2 q k 2 When the complementary elliptic modulus is used, it is 1 − q in Ref. [10] and in this paper, and it is k ′ 2 = 1−k 2 in Refs. [11,12,13]. In particular, notice that the instanton parameter q enters elliptic potentials differently for N=2 * theory and N f =4 theory. For the N=2 * theory, q appears as the nome of Weierstrass elliptic functions; for the N f =4 theory, q appears as the modulus of Jacobian elliptic functions.

Critical points and asymptotic solutions
The connection between N f =4 super QCD and Schrödinger equation relies on the fact that instanton partition function of Ω-background deformed N=2 gauge theory with surface operator satisfies a partial differential equation, in the Nekrasov-Shatashvili (NS) limit the equation becomes an ordinary differential equation [14,15]. In the case of SU(2) N f =4 gauge theory, the ordinary differential equation is the Heun equation in normal form [16].
Another way to see this connection is through their relations to conformal field theory (CFT). As a special case of the correspondence discovered by Alday, Gaiotto and Tachikawa [17], the N f =4 supersymmetric QCD model is related to the 4-point sphere conformal block of Liouville CFT. An extension of this relation connects supersymmetric QCD model with surface operator to conformal block with a null vertex operator inserted. A conformal block including a null vertex operator satisfies a partial differential equation, the Belavin-Polyakov-Zamolodchikov equation [18]. The quantum mechanical models are obtained by taking a further limit, it is the large central charge limit for CFT. For the case of 4-point sphere conformal block with a null vertex operator, it also leads to the Heun equation in normal form [19]. Transforming the Heun equation to elliptic form, it is a Schrödinger equation with the DTV potential. In this paper we use the elliptic form of the equation to study the spectral problem. Using Weierstrass elliptic function, the Schrödinger equation is with potential The periods 2ω 1 , 2ω 2 are generators of the period lattice, the other two periods are 2ω 0 = 0 and 2ω 3 = 2ω 1 + 2ω 2 . We use the shifted Weierstrass elliptic function ℘(x; 2ω 1 , 2ω 2 ) = ℘(x; 2ω 1 , 2ω 2 ) + ζ 1 , because wave function for the shifted potential is more convenient to compare with instanton partition function of gauge theory. The constant ζ 1 is given by Weierstrass zeta function ζ(x; 2ω 1 , 2ω 2 ) at x = ω 1 , ζ 1 = ζ(ω 1 ; 2ω 1 , 2ω 2 )/ω 1 . The nome of Weierstrass elliptic function is p = exp 2πi ω 2 ω 1 . The potential can be represented in terms of elliptic theta functions as using the relation of elliptic function ℘(x, p) and theta function ϑ 1 ( πx 2ω 1 , p). Using Jacobian elliptic function the equation is transformed to with potential u(η) = b 0 q sn 2 η + b 1 q cn 2 η dn 2 η + b 2 1 sn 2 η + b 3 dn 2 η cn 2 η .
The modulus of Jacobian elliptic functions is q. The periods of potential u(x) are 2K, 2iK ′ and 2K + 2iK ′ , they are half periods of Jacobian elliptic functions. In the rest of this paper, periods/half-periods always refers to that for the potential. In later usage, the variables associated with Eq. (1) are called Weierstrass variables, the variables associated with Eq. (4) are called Jacobian variables. The elliptic nome and modulus are related by q = ϑ 4 2 (p)/ϑ 4 3 (p), other quantities are related by The Floquet exponent for Eqs. (1) and (4), which appear in later sections, denoted by ν and µ respectively, are related by νω 1 = µK. The two forms of equation and its solutions are equivalent, however, studying the solutions needs various manipulations of elliptic functions, one of the forms would be preferred when it is easier to perform computations.
The Heun equation appears in a wide range of physics and mathematics topics, various methods are used to find different kind of solutions. The solutions relevant to the N f =4 super QCD model we seek in this paper are of a particular kind, they are asymptotic solutions with coupling constants b s taking generic values, and the ratio of energy λ (or Λ) and coupling constants b s (or their geometric average) determines the perturbation solution be weak coupling or strong coupling. The DTV potential is a four components generalization of the Lamé potential, which in turn is the elliptic generalization of the Mathieu potential. Asymptotic solutions of the Mathieu potential has been known for many years, asymptotic solutions of the elliptic generalized Lamé potential are also known now [20, 21,12,13]. However, asymptotic solutions of the same nature for Heun equation become difficult to compute. For what is known about Heun equation up to recent time, see the books [22] and [7]. Among recent attempts, a study using WKB analysis can be found in Ref. [23], solutions are determined by conditions in algebraic form represented by equations of Bethe-Ansatz type in Ref. [24]. In this paper we study the spectral problem from a different point of view, by considering its relation with 4-dimensional N=2 supersymmetric gauge theory [5]. We would use accessible methods from both sides to study various aspects of this connection. For example, in the rest of this section we explain a fact about the global structure of the asymptotic solutions of the DTV potential. The elliptic potential has multiple critical points, every critical point is associated with an asymptotic spectral solution. This phenomenon has been verified for spectral problems of some simpler elliptic potentials, it seems a generic property for a class of elliptic potentials. It is related to the electro-magnetic-dyonic duality of the low energy effective theory of N f =4 super QCD model [8,9]. It is difficult to see this property by conventional approach in spectral theory. Now we explain in detail the relations of critical points of the DTV potential, weak and strong coupling asymptotic solutions, and Seiberg-Witten duality of the N f =4 super QCD.
This fact was studied in Ref. [10], for consistency we recall the results here which would be used later in this paper. We have shown there are six critical points for the DTV potential, whose locations can be computed recursively as follows. The critical points of the potential, in Darboux form (5), are given by solutions of the equation ∂ η u(η) = 0. It is a 12th order polynomial equation for the function snη, using the variable z defined by q sn 2 η = z, the equation can be rewritten as Q 6 (z) = 0 with Q 6 = C 6 z 6 + C 5 z 5 + C 4 z 4 + C 3 z 3 + C 2 z 2 + C 1 z + C 0 , where On the gauge theory side the UV theory is assumed weakly coupled, therefore q ≪ 1, we have this condition for all perturbative computations carried out later. Then the equation Q 6 (z) = 0 can be solved perturbatively by treating q as the expansion parameter. There are six distinct solutions, the leading order coefficients of the solutions are given explicitly bellow, other coefficients can be obtained by an iterative process as shown in Ref. [10].
The first pair of critical points are sn 2 η * 1,2 = b where u(η * 1,2 ) = (b where u(η * 3,4 ) = (b The expressions (10) and (12) indicate the first four critical points are on the same footing, they are converted to each other by permutations of b s , or by permutations of physical masses µ i according to relations given in Eq. (27) . As shown in Ref. [10], there is a single expression that represents perturbative spectral solutions around η * 1,2,3,4 , it is the weak coupling solution of the DTV potential, related to super QCD model in the weak coupling region of low energy moduli space.
There remains the last pair of critical points, where u(η * 5, In fact, coefficients s 5;n and s 6;n are the same, therefore sn 2 η * 5 and sn 2 η * 6 are related to each other by a simple sign inversion q The coefficients s 1;n , s 2;n , s 3;n , · · · s 6;n are rational functions of b 1 2 s , in later computations of gauge theory quantities they can be written as rational functions of mass parameters µ i . In the mainbody of the paper we study models with generic coupling constants b s (equivalently generic masses µ i ), so that the special values b 0 = b 1 or/and b 2 = b 3 are not applied to formulae such as (9)- (14). The analytical method applied here relies on the series expansion of various quantities which are not valid for these values. These values are of course allowed, they are related to some degenerate cases of the models separately discussed in Section 6.
A fact that would be useful later is the distance between strong coupling critical points.
When b s take generic values, which means no coefficients s 5,6;n in the expansion (13) take extremal small or large value, sn 2 η * 5,6 are dominated by the leading order terms which are As variables in elliptic potential and super QCD take complex values, it is their magnitudes which are compared. It indicates the distance η * 6 − η * 5 is very close to the half period K. In fact, it can be shown However, when the coupling constants take non-generic values such as approaching each other as b 0 → b 1 or/and b 2 → b 3 , the coefficients of later terms in Eq. (13) might blow up, the argument above fails, and the distance η * 6 − η * 5 could take any value in the period lattice. See subsection 6.1 for discussions about cases when the coupling constants take special values related to Landen transformations. We would use "weak coupling" and "strong coupling" for both the Schrödinger equation and the effective theory of super QCD model. While their couplings are characterised by very different variables, it happens that the weak/strong coupling solutions of spectral problem is related to the weak/strong coupling expansions of effective gauge theory. So in this paper terms like "weak/strong coupling solution/expansion" should be understood within the content where results about elliptic potential/super-QCD models are stated.

Weak coupling spectrum and N f =4 super QCD
In this section we compute the weak coupling (or large energy) spectrum of the DTV potential. In the end of this section we discuss their relation with the deformed prepotential of N f =4 super QCD model with surface operator.
By "weak coupling" solution of the spectral problem, it means the kinetic energy characterised by the quasimomentum ν is much larger than the potential energy characterised by the coupling constants b s , ν ≫ b 1 2 s . The quasimomentum is proportional to the scalar v.e.v < φ >= a in N=2 super QCD, and the coupling constants given in Eq. (27) are linear combinations of the square of flavor masses µ i , therefore, the weak coupling spectral solution corresponds to super QCD in the electric region where a ≫ µ i .

The spectrum at weak coupling
We use equation in the Weierstrass form (1) and (2) to compute the weak coupling solution.
The perturbation method for weak coupling solution is closely related to the procedure of generating infinity many conserved charges of KdV hierarchy [25,26]. Substituting the wave function ψ(x) = exp(´v(x)dx) into Eq. (1) leads to an equation for the integrand v(x), The weak coupling expansion for v(x) is equivalent to the large energy expansion, it takes the form where the coefficients are given by KdV Hamiltonian densities v ℓ (u, ∂ x u, ∂ 2 x u, · · · ). Periodicity of the potential requires the wave function to satisfy monodromy relation ψ(x + 2ω 1 ) = exp(±i2ω 1 ν)ψ(x), where the Floquet exponent ν is a constant depending on parameters λ, b s and q. As we have confirmed for other examples of elliptic potentials [11,12,13], the period associated to the weak coupling solution is 2ω 1 , this fact is also true for the DTV potential. So the relation of Floquet exponent and integral of v(x) over periodic 2ω 1 is given In connection with gauge theory, the integral over 2ω 1 is related to the α-cycle integration of Seiberg-Witten form on the curve of N f =4 super QCD which gives the scalar v.e.v a [8,9]. The spectrum of the differential equation always includes two sectors labeled by "±" signs, they correspond to performing the contour integrals clockwise/counter-clockwise in gauge theory.
The weak coupling eigenvalue and eigenfunction are obtained from definite and indefinite integrals of v(x). For the DTV potential, some relations of the shifted Weierstrass elliptic functions are needed to simplify expressions of v ℓ , they are given in Appendix A.
The eigenvalue is computed from the monodromy relation (17), expressions of v ℓ suitable for the integration are given in Ref. [12]. As v 2ℓ are total derivatives of periodic functionals, their integrals over a period vanish. For integrals of v 2ℓ−1 , the total derivative part of the integrand does not contribute, the remaining part contributes non-vanishing results, denoted by ε ℓ are expressed in terms of b s and modular functions e 1 , e 2 , e 3 , ζ 1 . The eigenvalue expansion is obtained by reversing the relation (17), the coefficients λ l are polynomials of ε ℓ . The expressions of λ l quickly become complicated, the first two of them are In the above summation the indices i, j, k take values in {1, 2, 3}. The summation restricted by i < j < k in fact requires i = 1, j = 2, k = 3. The summation restricted by i|j < k means when the index i takes a value in {1, 2, 3} the indices j, k take values in the remaining set and satisfy j < k.
The corresponding unnormalized wave functions are obtained by performing indefinite integral´v(x)dx for the integrand (16) and substituting λ 1 2 (ν) derived from the dispersion relation (19). They can be written in the form where the shifted Weierstrass zeta function ζ(x) is defined by ζ(x) = ζ(x) − ζ 1 x, it is a 2ω 1 -periodic function. The coefficients c l,s,d have three indices, among them l refers to the exponent of ν, s is in accordance with ω s , and d refers to the order of derivative of ζ(x). The coefficients appearing in the above are , Notice that the wave function (21) is already in Floquet form  [8,9]. Besides local operators, there are finer probes for gauge theories such as line operators and surface operators which detect more structure [27,28,29].

Instanton partition function of
For N=2 supersymmetric gauge theory there is another refinement, the Omega-background field is introduced to deform gauge theory, the moduli space of instanton is tamed so that the localization technic can be used to compute instanton partition function [30,31,32,33]. When the surface operator preserves supersymmetry, such as a plane surface in the x 1 − x 2 direction of spacetime, the instanton partition function with surface operator can be computed as well [34,35,36,37,38,39,14]. Correspondingly, the instanton counting parameter q splits up into multiple counting parameters, for SU(2) gauge group there are two, x 1 , x 2 (not the spacetime coordinates) subject to a relation x 1 x 2 = q. The SU(2) N f =4 super QCD is a prototype, relevant properties of instanton partition function with surface operator are discussed in Refs. [36,37], using the equivariant character developed in Ref. [40]. In the context of Gauge/Bethe correspondence, instanton partition function with surface operator in the NS limit contains complete spectral information of the associated quantum mechanical problem [5,36,19,37,16].
The expressions of instanton partition function for N f =4 super QCD are pretty lengthy, but a few explicit expansions are already enough to explain the points. The more useful quantity here is the logarithm of instanton partition function, the prepotential. So writing in exponential form and manifesting the pole structure of deformation parameters ǫ 1 , ǫ 2 , the instanton partition function takes the general form The function F inst is a polynomial of ǫ 2 , therefore the pole structure alone does not fully fix F inst and G inst . Functions F inst and G inst are both expanded as double series of parameters x 1 , x 2 , a further requirement about x 1 , x 2 determines them uniquely. The first piece is required to contain only terms with the same exponent for x 1 , x 2 , i.e. it is a q-series. In order to compare with the weak couplings spectral solution which is written as large ν series, we expanded the prepotential as large a series. The first piece is given by with F l,k rational functions of µ i and ǫ 1 , ǫ 2 . When ǫ 1 , ǫ 2 are assigned to mass dimension one, F l,k has mass dimension 2l + 2. The second piece is required to contain terms with different exponents for x 1 , x 2 , it is expanded as with G l,k 1 ,k 2 also rational functions of µ i and ǫ 1 , ǫ 2 , has mass dimension 2l + 1.
In fact, besides the instanton contribution, the full gauge theory partition function also contains classical and perturbative contributions, Z = Z cl Z pert Z inst . Correspondingly, the full prepotential is given by F = F cl + F pert + F inst and a similar formula for G. In the rest part of this subsection we show that, in the NS limit ǫ 2 → ǫ, ǫ 2 → 0, the function F is related to the eigenvalue (19) and the function G is related to the unnormalized wave function (21). As it would be clear, the relation between spectrum of the DTV potential and prepotential of the N f =4 super QCD is not quite straight to identify because the weak coupling spectrum is expressed in Weierstrass variables while the instanton partition function is primarily expressed in Jacobian variables, this causes some extra twists compared to a simpler example discussed in Ref. [13].

Eigenvalue from instanton partition function
The variables easy to identify are coupling constants b s and flavor masses µ i . By deriving the differential equation through their connection to CFT, the relation of b s and µ i can be fixed, see e.g. Ref. [19], The coupling constants can be written in the form b s = n s (n s − 1), with n s related to µ i by For gauge theory quantities n s are more primary objects because a few coefficients in the prepotential, in formulae (31) and (41) given below, do not involve combinations n s (n s − 1). For some quantities discussed in this paper, their properties could be determined by the leading order of ǫ-expansion. Therefore, the leading order values of b s are also useful, they For the prepotential F , the classical part is given by F cl = −a 2 ln q. The perturbative part is a one-loop effect and can be computed as regularized universal denominator extracted from Nekrasov instanton partition [41], it is a q-independent function therefore irrelevant to the discussion here. Up to a q-independent piece, the prepotential F is From the point of view 2-dimensional field theory on the surface defect, F in the NS limit is the effective twisted superpotential [5,42,43]. Notice that there are two elliptic variables involved, gauge theory quantities use the modulus q, spectral data in Weierstrass variables use the nome p. In order to see their connection, we rewrite the prepotential as a series in p using the relation q = ϑ 4 2 (q)/ϑ 4 3 (q). Then the leading order coefficient of large a expansion of F is F −1 = − 1 2 a 2 (ln p − 8 ln 2). For N f = 4 gauge theory, the instanton contribution F inst is an infinite series −a 2 ( 1 2 q + 13 64 q 2 + 23 192 q 3 + · · · ) at the leading order, combined with an extra piece from the classical contribution −a 2 ln q, it leads to F −1 in a single term because 1 2 ln p = 4 ln 2 + ln q + 1 2 q + 13 64 q 2 + 23 192 q 3 + · · · . As the leading order of large ν expansion of λ has a single term −ν 2 , therefore, when the a and ν is identified as πa = ǫω 1 ν, the relation at this order is simply 2π 2 p ∂ ∂p F −1 = −ǫ 2 ω 2 1 ν 2 . The nome p has a role in gauge theory. The massive N f =4 super QCD can be constructed starting from the massless/superconformal theory and adding mass deformation, ln p is the complex coupling of Coulomb branch effective theory for the massless theory [9]. The low energy effective coupling is renormalized by instanton effects, encoded in the relation of p and q [44,17].
In the next order, there is an a-independent piece in the prepotential F , there is no counterpart in eigenvalue. This piece is summarised in the function Φ given below. From order a −2 in the expansion of prepotential, when q is substituted by the series relation q = 16p 1 2 −128p+704p 3 2 +· · · , up to a scaling factor p ∂ ∂p F l precisely matches the corresponding order coefficient λ l in the expansion of eigenvalue.
In summary, the relation between the eigenvalue λ and the prepotential F is with In Eq. (30) the left side is expressed in terms of (quasi)modular functions e i , ζ 1 with nome p, while the right side is expressed as a q-series. Therefore, instead of using derivative p ∂ ∂p on the right side, we use

Wave function from instanton partition function
First we need to identify relation between the instanton parameters x 1 , x 2 and the coordinate of wave function x. The condition x 1 x 2 = q does not give us enough information. The next hint comes from the leading order coefficient of large a expansion of G inst in the NS limit. The function G inst −1 is an infinite series in x 1 , x 2 expanded as − 1 2 x 1 + 1 2 x 2 + 3 16 x 2 1 − 3 16 x 2 2 + · · · , but the leading order term in the (logarithm of) wave function contains only one term −iνx. The situation is similar to the case of matching F and the eigenvalue λ. Inspecting the series expansion of G inst −1 , we conclude that an extra piece is also needed to rewrite G inst −1 into a single term. The extra piece can be interpreted as the classical contribution of surface defect, the only reasonable form is 1 2 a ln(x 2 /x 1 ), up to terms independent of x 1 , x 2 . Therefore, the prepotential G is Then the leading order coefficient is x 4 1 x 2 2 + 35 2048 x 2 1 x 4 2 + 7 512 x 1 x 5 2 + 77 2048 Although the numerical coefficients of G −1 carry sufficient information to determine the functional relations x 1 (x) and x 2 (x), it still looks quite unlike the expected single term −iνx.
We directly provide the result, the details are explained in Appendix B. Define normalized shifted coordinates x ± = π 2ω 1 (x ± ω 2 2 ), then x 1 and x 2 are related to x through elliptic theta function, The property x 1 x 2 = q follows from transformation of theta function under translation The series expansion of them are Substituting the expansion into G −1 leads to G −1 = iπx/ω 1 . Therefore, the leading order piece of prepotential G is precisely the first piece in the wave function. The next order piece G 0 /a 0 is a-independent, there is no counterpart in the wave function.
It is given by Substituting (35) and using the relation ϑ 2 Then the next leading order of prepotential is up to constant terms which do not affect comparison with the unnormalized wave function. Or written in terms of the same shifted coordinate it is The third order coefficient, after arranging every term in the right place, is a series expanded as Nicely enough, each piece in brackets is summed up in a theta function, in sequence they are iω 1 2π ∂ x ln ϑ r (x + , p) for r = 1, 2, 4, 3. In fact, when substituting (36) into G 1 there are some constant terms appear, again they are neglected for comparison with the unnormalized wave function. So up to constant terms we have The shifted zeta-function is give by theta function as ζ(x) = ∂ x ln ϑ 1 (πx/2ω 1 , p), therefore, the third order piece in G precisely gives terms in the corresponding order of wave function, .
The fourth order coefficient shares the same structure as the previous one, it is It is recognized that G 2 ∼ ∂ x G 1 , again up to constant terms, Therefore the fourth order contribution of G gives a further piece in the wave function, In summary, under the identification of parameters the wave functions ψ ± (x) are related to the prepotential G through the relation Their relation is denoted by "⇐" because when substituting variables x 1 , x 2 on the right side by relation (48) some constant terms are dropped out to match the left side.

Strong coupling spectra and gauge theory duality
In this section we compute perturbative expansion of the strong coupling spectrum for the DTV potential. We use the Jacobian form of the equation for this purpose. The strong coupling spectral solution should be generalization of the corresponding solution of the Lamé potential and ellipsoidal potential. However, the method of perturbative computation for strong coupling spectra used in Refs. [12,13] does not straightly apply for the DTV potential. Substituting wave function ψ(η) = exp(´v(η)dη) into Eq. (4) we get the relation for v(η), v η + v 2 = u + Λ. With the potential u(η) in its original form (5), it is difficult to find the strong coupling expansion for v(η) that solves the nonlinear relation. In fact, as discussed in Refs. [10,12], in this case the strong coupling expansion parameter is which does not naturally appear in any simple manipulation of u(η). Before going to details of how the problem is solved, there are a few properties of the strong coupling solutions we can expect based on results obtained earlier for other elliptic potentials.
Similar to other examples, we should use the Jacobian form of the spectral problem, Eqs. (4) and (5), to study the strong coupling solution, so that expressions of the eigenvalue and eigenfunction would appear in compact form. Accordingly, the Floquet exponent is denoted by µ, which is different from ν of the weak coupling solution. In connection with gauge theory, µ is proportional to the magnetic dual a D of the scalar v.e.v, so µ is dual to ν in the sense of electric-magnetic duality of N=2 super QCD [8,9]. Consequently, computation of the strong coupling spectral solution also provides necessary materials to derive the dual expansion of prepotential for N f =4 super QCD with generic masses, presented in Section 5.

Effective potential at magnetic critical points η * 5
The leading order (in the ǫ-expansion) eigenvalue was computed in Ref. [10] using the normal form of Heun equation. Result obtained there indicates the coupling constants b s , now assumed large, are not the correct expansion parameter. The correct perturbation parameter is a single quantity of order Meanwhile, the periodicity of elliptic potential is not affected by perturbative treatment. This fact suggests that in the strong coupling region the four components potential u(η) can be "averaged" so that its behavior is similar to a single component elliptic potential with periods 2K and 2iK ′ .
It is helpful to learn from some simpler examples studied before [12], to have a clue about how to "average" the potential u(η). For elliptic potentials having multiple critical points, the strong coupling spectral solutions are eigenstates with small energy (complex valued) compared to the strength of potential. These eigenstates are spatially localised near the corresponding critical points of potential. That means the critical points looks like local vacuum, and the potential near the critical points behaves like a binding potential. To see more clearly the shape of the potential near a critical point we should change to a new coordinate system that the critical point is at its origin, then expand the potential near the critical point, so that the potential is rewritten as an "effective potential" in the new coordinate. The Fig. 1 schematically illustrates the local coordinate systems near the strong coupling critical points η * 5 and η * 6 , using a real potential with multiple critical points. In the region near η * 5 , when rewriting the DTV potential in the local χ-coordinate, as we show in the below, the effective potential contains a single dominant term with the largest coupling constant, therefore perturbation can be performed. The situation is similar for the region near η * 6 when rewritten in the local χ ′ -coordinate. A fact should be stressed, the effective potential used to compute spectrum at the critical point η * 5 can not be used to compute spectrum at the critical point η * 6 , and vice versa. For the Lamé potential and the ellipsoidal potential discussed in Refs. [12,13], they are not generic enough to fully reveal this feature because their strong coupling critical points are at special values of coordinate so that their spectral solutions can be computed in a single coordinate system, see the discussion in Appendix F. Regarding the relation with N=2 super QCD, change coordinate systems for the spectral problem is analogical to change degrees of freedom for the corresponding gauge theory models.
In quantum field theory, from the viewpoint of effective theory, the proper degrees of freedom used to formulate the theory vary within the parameter space. In particular for some field theory models, particles and solitons are used to formulate the effective theory at weak and strong coupling regimes, respectively. The duality phenomenon of this kind, which the N=2 N f =4 super QCD model relies on [8,9], constitutes an important development of current understanding of quantum field theory.
In the rest of this section we compute the effective potentials at the critical points η * 5 , η * 6 and their spectral solutions. To expand the DTV potential near the critical points η * 5 , introduce the coordinate χ, and the corresponding energy above local vacuum δ by where Λ * 5 = −u(η * 5 ), so that δ is the small perturbation near η * 5 , it satisfies the condition δ ≪ Λ * 5 . Then Eq. (4) is rewritten as with the effective potential u(χ) given by In order to do perturbative computation, we need to simplify the right side of the above expression using the value of sn 2 η * 5 . To keep periodicity of the potential, the effective potential u(χ) must be an elliptic function made out of snχ, cnχ, dnχ, and it should not contain terms with pole at χ = 0, K, iK ′ , K + iK ′ because the potential u(η) is finite at these points.
As the perturbation results at η * 5 and η * 6 are related in simple way by a sign inversion q 1 2 → −q 1 2 , in this section we only need to discuss the perturbation at η * 5 . The perturbation at η * 6 is obtained by the same sign inversion. So in the following discussion we drop numerical subscriptions.
The periodicity of the potential play a crucial role in the process of deriving expression of u(χ). The periods of potentials u(η) and u(χ) are 2K, 2iK ′ , 2K + 2iK ′ , we need to keep every elliptic function expression with the right periodic behavior when changing their forms. For example, on the right side of Eq. (52), there are terms sn 2 (η * + χ) − sn 2 η * , to simplify it we need the addition formula for elliptic function, Then sn 2 (η * + χ) is expressed as product of two parts. The factor (1 − q sn 2 η * sn 2 χ) −2 can be expanded as a q-series, in every term the elliptic function is a non-negative integer power of sn 2 η * sn 2 χ. The factor ( snη * cnχ dnχ + snχ cnη * dnη * ) 2 contains three terms sn 2 η * cn 2 χ dn 2 χ, The first two terms can be simplified to expressions containing only non-negative integer powers of sn 2 η * and sn 2 χ, using functional relations of Jacobian elliptic functions. The third term need careful treatment. The value of snη * , cnη * and dnη * can be obtained from the value of sn 2 η * , taking either + or − sector for square roots. But we shall not substitute cnχ and dnχ simply by ±(1 − sn 2 χ) 1 2 and ±(1 − q sn 2 χ) 1 2 to rewrite snχ cnχ dnχ as because an inconsistence of periodicity arises by this substitution. The left side of (55) appears in a term of the potential u(χ) in the form (· · · ) snχ cnχ dnχ, where the factor in bracket is a constant coefficient determined by snη * , the factor snχ cnχ dnχ is invariant under translation of χ by periods 2K, 2iK ′ and 2K + 2iK ′ . But the substituted expression on the right side of (55) does not maintain periodicity for translation over 2K because sn(χ + 2K) = − snχ, and cnχ, dnχ are invariant, so the sign flips. The situation is similar when simplifying other terms in (52), such as cn 2 Therefore, whenever a term containing powers of cnχ and dnχ arises, the correct substitution rule preserving periodicity of the potential are the following, for integer m 1.
The effective potential obtained in this way contains infinitely many terms, where the effective coupling constants Ω 2m and Ω 2m+1 are functions of b s and q, determined by snη * . The coupling constants are q-series of mass dimension zero, by counting powers of q we have Ω 2 ≫ Ω 3 ∼ Ω 4 ≫ Ω 5 ∼ Ω 6 · · · , so among the infinitely many coupling constants there is a single largest one which can be used as the perturbation parameter. Moreover, in the region −K χ K, | snχ| 1, so that later terms in the effective potential with increasing m-index become less important. Therefor, the effective potential is indeed very close to the Lamé potential, u(χ) ∼ Ω 2 sn 2 χ, but slightly corrected by smaller terms. This fact explains why the structure of strong coupling solutions of the DTV potential, obtained in this section, are the same as solutions of the Lamé potential, with higher order perturbations Ω m 3 included accordingly. Also notice that Ω 2 and Λ * 5 have the same magnitude, Ω 2 ∼ 2Λ * 5 .
For purpose of computing the dual prepotential of N f =4 super QCD, in Appendix C we give explicit expression for the first few coupling constants Ω m in the leading order ǫexpansion, in terms of physical masses using the relation of b s and µ i given in Eq. (27).

The asymptotic eigenvalue
Now we can proceed to find the strong coupling solution for the potential (57) associated to its critical point χ * 5 = 0 where u(χ * 5 ) = 0. The existence of a single dominant coupling constant of u(χ) makes equation (51) in the same class of elliptic spectral problem we have studied previously, including the Lamé equation and ellipsoidal wave equation [12,13], method used there can be applied straightforwardly.
Writing wave function in exponential form ψ(χ) = exp(´v(χ)dχ), Eq. (51) leads to the So we are seeking a strong coupling expansion of v(χ), with Ω 1 2 2 as the large parameter. It has a solution in the form it includes two sectors depending on the choice of ± sign. During the process of solving v ℓ (χ) recursively, whenever a term involving cn n χ or dn n χ with n 2 arises, then apply the substitution rule (56). The first few of them are There are two noticeable properties for the expansion of v ℓ (χ). The coordinate dependence is represented solely by powers of function snχ, this is crucial for discussion in subsection 4.5.3 about the duality transformation of wave functions. Also notice that this expansion is valid near the point χ ∼ 0, so that v ℓ (χ) are represented as Laurent series of snχ, this is what needed to perform the contour integrals in the next subsection. The general form of even v 2ℓ can be written as total derivatives, In the summations, the coefficient of constant terms c ′ 2ℓ,0 and c ′′ 2ℓ,0 are set to zero, and c ′ The eigenvalue δ(µ) at the strong coupling point η * 5 is given by the monodromy of wave function over period 2iK ′ . However, as we have shown for other examples in Refs. [11,12], the definition of Floquet exponent µ does not follow a direct application of classical Floquet theory. The correct relation of µ and the monodromy of wave function is the same as for other elliptic potentials we have studied, An explanation of the relation is given in subsection 4.3 by rewriting the corresponding wave function in canonical coordinate. As the DTV potential is a more general elliptic potential that reduces to simpler potentials in proper limits, this relation also guarantees its spectral solutions reduce to solutions of those potentials. The integral over period 2iK ′ in the χ-plane is mapped to the integral along β-contour in the ξ-plane defined by ξ = sn 2 χ, as explained in Appendix B of Ref. [12]. It is related to the β-cycle integration of Seiberg-Witten form on the curve of N f =4 super QCD [8,9].
As shown in Eq. (60), v 2ℓ are derivative of series in snχ which are invariant under translation χ → χ + 2iK ′ , therefore we havê As shown in Eq. (61), each v 2ℓ+1 contains two parts. In the first part, the series are made out of sn m χ with m = ∓(2n + 1), so the corresponding integrals are given by According to the prescription of contour integration compatible with Floquet property, described in Refs. [11,12], for positive exponent I 2n+1 = 0, and for negative exponent I −(2n+1) = 0. In the second part, the series in bracket made out of sn 2n+1 χ are invariant under translation χ → χ + 2iK ′ , hence their integrals over 2iK ′ vanish. Therefore, in the series of v 2ℓ+1 given by (61), terms in the first part with 2n + 1 < 0 and entire of the second part do not contribute to integral, only terms in the first part with 2n+1 > 0 (i.e. the principle part of the Laurent series) contribute nonzero integrals. Dropping all terms of vanishing integral, the Floquet exponent (62) is computed by The integrals I −(2n+1) can be computed by recurrence relations, with a choice for integral contour in the ξ-plane as explained in Appendix B of Ref. [12]. The summation only involves odd integer powers of ±Ω 2 leads to the sign inversion for µ. The integral (64) gives the functional relation µ(δ) as a large Ω 2 expansion, reversing the series leads to the strong coupling dispersion relation. Up to order O(Ω −2 2 ), the expansion of + sector eigenvalue is The − sector eigenvalue is obtained by inversion Ω This spectral solution is considered as a strong coupling solution in the sense that the strength of potential is much larger than the eigenvalue, Ω 2 ≫ δ, or equivalently characterised by the condition Ω 1 2 2 ≫ µ. The eigenvalue (65) generalises the strong coupling solution of other elliptic potentials such as the Lamé potential and ellipsoidal potential.
can be used as the strong coupling expansion parameter. In the spectral problem Ω 2 is a more convenient expansion parameter for eigenvalue (65) and eigenfunction (66), but in the related problem about dual expansion of prepotential for N f =4 super QCD, discussed in Section 5, the proper expansion parameter is

The asymptotic eigenfunction
The wave function ψ(χ) = exp(´v(χ)dχ) associated with the eigenvalue (65) The wave function provides a more transparent explanation for the monodromy relation (62). According to the prescription of integral over period 2K explained in Appendix B of Ref. [12], the integral trajectory in the χ-plane crosses the branch cut of ln[( dnχ+ cnχ)/ snχ], giving a monodromy iπ, but it avoids the branch cut of ln( dnχ − q 1 2 cnχ). The remaining terms in the wave function, such as ln snχ, sn m χ, cn 2m χ and cnχ dnχ, with m ∈ Z, are invariant under translation χ → χ+2iK ′ . Therefore, the monodromy of wave function comes from a single term ±iµ ln[( dnχ+ cnχ)/ snχ], the remaining part of wave function is periodic. This property indicates ψ ± (χ) given by (66) are Floquet wave functions. Considering it as the wave function of a quantum particle, µ is interpreted as the (quasi)momentum, so we could define the canonical coordinate ρ conjugating to µ by The factor (1 − q) 1 2 is included so that all terms of elliptic functions in the wave function can be expressed in terms of hyperbolic functions sinh ρ and cosh ρ, it causes a coordinate independent constant that does not change the unnormalized wave functions. Then the wave functions take standard Floquet form . See more details in Appendix E. The wave functions written in canonical coordinate explain a puzzling point in previous study [11,12,13]. There the monodromy relation (62) for elliptic potentials was determined by requiring the solutions of elliptic potential has correct limit to the solution of trigonometric potentials when taking the limit q → 0 and fixing other coupling constants properly. This requirement refuses the naive relation between Floquet exponent and monodromy of wave function that follows from coordinate translation χ → χ + 2iK ′ , This is because µ is really conjugate not to the original coordinate χ but to the canonical coordinate ρ, so the correct monodromy relation is It is equivalent to the relation in χ-coordinate which is the monodromy relation (62). A similar argument applies to strong coupling solution at the dyonic region obtained in subsection 4.5, whose wave functions are given in Appendix D.
In the strong coupling eigenvalue (65) and eigenfunction (66), higher order coupling constants Ω m 5 would enter the expressions in later terms. Setting Ω 2 , Ω 4 nonzero and all other coupling constants vanishing, we get the corresponding solution of ellipsoidal wave equation. Appendix F contains a discussion on the duality relation of its two strong coupling solutions.
At this point it is clear that the elliptic Hill potential (57) and its solutions generalise familiar examples of the Lamé equation and ellipsoidal wave equation, the procedure of solving v(χ) and preforming contour integral follows the prescription given in Ref. [12]. However, as would be shown in subsection 4.5, there are some particular features of the DTV potential that do not manifest for other elliptic potentials.

Leading order epsilon expansion
Besides Ω 1 2 2 , there is another parameter ǫ can be used as the expansion parameter. It is the Omega background deformation of N=2 super QCD that survives in the Nekrasov-Shatashvili limit. In the related quantum mechanical problem, it is interpreted as the Plank constant [5]. In the previous treatment ǫ is hidden in coupling constants b s , so it does not explicitly appear in various quantities of spectral solution. Counting the power of ǫ, the coupling constants b s , Ω m and the eigenvalues λ, Λ, δ are of order ǫ −2 , the exponents ν and µ are of By computing the strong coupling solution of elliptic potential u(χ), it leads to a way to compute the strong coupling expansion of prepotential for N f =4 super QCD. The magnetic and dyonic dual prepotentials of super QCD are computed in Section 5, restricted to the leading order of ǫ-expansion. For that purpose here we compute the leading order spectral solution in the strong coupling region which would be used later.
Expanding v(χ) with respect to ǫ, the wave function is The power of ǫ has been extracted, so p (n) are functions of order O(1). At the leading order, the Floquet exponent µ is computed from To perform the integral, we expand the integrand as a large Ω The left side is of order ǫ −1 , and Ω 2 is of order ǫ −2 , so the coefficients P β (0);ℓ (χ) are of order ǫ −(ℓ+1) . To write P β (0);ℓ (χ) as a Laurent series of snχ, powers cn n χ and dn n χ with n 2 are substituted according to the rule (56). It turns out that all even coefficients vanish, P β (0);2ℓ (χ) = 0. The odd coefficients P β (0);2ℓ+1 (χ) are given by The Floquet exponent at the leading order is given by integral over period 2iK ′ , Only I −(2n+1) with n 0 contribute nonzero integrals. Therefore up to order O(δ 2 ) for the where "⇒" indicates terms with vanishing integral are dropped. The integration leads to relation Reversing the relation we get the eigenvalue at the leading order of ǫ expansion, The − sector is obtained by inversion µ → −µ.
The corresponding wave function at the leading order is This computation gives a faster way to access the later terms of spectral solution in the leading order ǫ-expansion, compared to extracting them from the full spectrum (65) and (66) by taking the limit ǫ → 0. For example, to get the leading order spectrum (79)   The second strong coupling solution is associated to the critical point η * 6 of the DTV potential (5), in the η-coordinate. In the χ-coordinate it is associated to the critical point χ * 6 = η * 6 − η * 5 of the effective potential (57). The effective potential (57) has infinity many terms, it is valid only near the critical point η * 5 where χ ∼ 0, and can not be directly used to compute the spectrum near η * 6 where χ ∼ χ * 6 . Although the magnitude of coupling constants Ω m decreases with increasing m-index when q ≪ 1 is assumed, the positive power of function snχ grows with increasing m-index when snχ 1, therefore the potential (57) is not guaranteed to make sense in the entire range of period lattice because snχ is not bounded. As a result, the potential is valid only when restricted to the region −K χ K where snχ 1. A simple fact can be used to illustrate this point, the critical point of the potential (57) are given by solutions of ∂ χ u(χ) = 0 where The condition ∂ χ u(χ) = 0 is an equation with infinite many terms involving cnχ, dnχ and sn m χ with all m ∈ Z + . It is easy to see χ * 5 = 0 is a solution, but it is impossible to determine other five solutions. In particular, the critical point χ * 6 , which is symmetrically related to χ * 5 in the η-coordinate by a simple sign inversion q which is also only locally valid. From experience with other elliptic potentials [12,13], the second strong coupling spectrum is related to the monodromy of wave functions over period 2K + 2iK ′ , so the effective potential at η * 6 should be a polynomial with infinitely many terms made of Jacobi cn-function. However, it is not correct to get the effective potential at η * 6 by simply substituting sn 2m χ → (1 − cn 2 χ) m in the potential (57) at η * 5 as and expanding it as a polynomial in cnχ. This is because the sum in potential (57) is valid when snχ ≪ 1 hence cnχ ∼ 1 is already assumed, so the expansion of cnχ obtained by the substitution (82) does not make sense as an asymptotic series. In the region near χ ∼ K the condition cnχ ≪ 1 is satisfied. Therefore, we should write the effective potential at η * 6 using another coordinate system, denoted by κ, in which the critical point η * 6 is at κ * 6 = K.
Then an asymptotic series of cnκ does make sense near κ * 6 because the condition cnκ ≪ 1 is satisfied. The argument above is caused by the fact that the effective potentials at η * 5 and η * 6 are infinite series of sn-function and cn-function respectively, they are valid only in restricted regions of period lattice, one near χ ∼ 0 and the other near χ ∼ K. The substitutions sn 2m χ → (1 − cn 2 χ) m and cn 2m χ → (1 − sn 2 χ) m , however, do not move the effective potential to another region. Therefore, the substitution (82) can not be used to transform the effective potential at η * 5 to the effective potential at η * 6 . The situation is different for elliptic potentials with finite number of terms, such as the ellipsoidal potential, where the substitution becomes valid. See the discussion in Appendix F.
The effective potential at η * 6 can be directly computed following a procedure similar to that in subsection 4.1. There is another approach, however, which is simpler and display more clearly the relation between effective potentials at η * 5 and η * 6 . Firstly, define a new coordinate by where Λ * 6 = −u(η * 6 ), and σ is the small perturbation near the critical point η * 6 , so that σ ≪ Λ * 6 is satisfied. The equation (4) can be rewritten as the spectral problem for an effective potential u(χ ′ ) with eigenvalue σ, where u(χ ′ ) = u(η * 6 + χ ′ ) − u(η * 6 ). The effective potential u(χ ′ ) can be simplified to an expression involving snχ ′ , cnχ ′ and dnχ ′ using the solution of η * 6 , taking into account the substitution rule given in (56). Because solutions of the critical points η * 5 and η * 6 only differs by a sign inversion q 1 2 ↔ −q 1 2 , therefore, the effective potential u(χ ′ ) takes the same form as the effective potential u(χ) given in (57) with the argument of effective coupling constants So near the critical point η * 6 the effective potential is where Ω ′ m (q 2 ), and the hierarchical structure of coupling constant persists. But it is not correct to continue the subsequent procedure in subsection 4.2 to compute the eigenvalue at η * 6 . In the approach adopted here two elements are crucial to compute asymptotic solutions of elliptic potentials, one is deriving the integrand under proper large parameter assumption, the other is choosing the correct integral contour, the integrand and the integral contour must fit into each other so that the procedure works. As shown in Ref. [12], the spectrum at η * 6 is obtained by performing integrals over 2K + 2iK ′ in the χ ′coordinate. The integrand obtained in subsection 4.2, expressed in sn-function, is associated to integral over 2iK ′ , but it is incompatible with integral over 2K + 2iK ′ . We have to work in the κ-coordinate promised before, then the integrands are represented by the elliptic cn-function and becomes compatible with integral over 2K + 2iK ′ . Figure 2: The κ-coordinate replaces the χ ′ -coordinate.
The η-coordinate and κ-coordinate are related by η = η * 6 + κ − K, so from (83) the χ ′ -coordinate and κ-coordinate are related by χ ′ = κ − K. The region near χ ′ ∼ 0 becomes the region near κ ∼ K, as desired. In order to rewrite the potential (85) in the κ-coordinate, the following relation with lower sign are used, Substituting these relations into the potential (85) leads to an expression involving negative integer powers of dnκ, the following substitution rule converts it to an expression in cnκ, Expanding the obtained expression with respect to cnκ, we get the effective potential valid near the point κ ∼ K as an infinite series of cnκ, The effective coupling constants ✵ m are combinations of the coupling constants Ω ′ m , for m n 1, where C m n are the binomial coefficients. The effective coupling constants also have hierarchial values ✵ 2 ≫ ✵ 3 ∼ ✵ 4 ≫ ✵ 5 ∼ ✵ 6 · · · , so the dominant one ✵ 2 is the large expansion parameter. In the region −K κ − K K, cnκ 1, so that terms in the effective potential become less important as the m-index increases. As the magnitude of The critical points of the potential (88) are given by solutions of ∂ κ u(κ) = 0 where In this case it has an obvious solution κ * 6 = K, but it is blind to other five critical points.
This confirms the expectation that the potential (88) is only locally valid near κ ∼ K.

Perturbative spectrum
The process of computing strong coupling solution for the potential (88) is similar to the case in subsections 4.1 and 4.2. There are + and − sector solutions for the relation v 2 (κ)+v κ (κ) = u(κ) + σ, given by large with v −1 (κ) = cnκ, v 0 = − 1 2 ∂ κ ln cnκ, etc. The second strong coupling spectral solution is related to the monodromy of wave function over period 2K + 2iK ′ , so the Floquet exponent is given by This is the same relation holding for a few other elliptic potential spectral problems [11,12,13]. It also does not follow from direct application of classical Floquet theory. The reason to include a factor (1 − q) 1 2 is explained further in Appendix D by rewriting the corresponding wave function in canonical coordinate. The integral over period 2K + 2iK ′ in the κ-plane is mapped to the integral along γ-contour in the ξ-plane defined by ξ = sn 2 κ. It is related to the γ-cycle integration of Seiberg-Witten form on the curve of N f =4 super QCD [8,9].
For computation of integrals in (93), the integrands v ℓ (κ) can be simplified. Details of computation and explicit expressions of v ℓ (κ) are given in Appendix D. The even integrand v 2ℓ can be written as derivatives ∂ κ (· · · ) + ∂ 2 κ (· · · ) where the expression in brackets are series made of ln cnκ, cn 2n κ with both positive and negative integer n. These series have vanishing monodromy under translation of coordinate κ → κ + 2K + 2iK ′ , so that integrals of v 2ℓ vanish,ˆκ The odd integrand v 2ℓ+1 also contains two parts. The first part contains series made out of cn m κ with m = ∓(2n + 1), and the corresponding integrals are given by J m = κ 0 +2K+2iK ′ κ 0 cn m κdκ. By the same contour integral argument used in Ref. [12], similar to the case of strong coupling spectrum at η * 5 , for positive exponent J 2n+1 = 0, and for negative exponent J −(2n+1) = 0. The second part can be written as a derivative ∂ κ (· · · ), the series in bracket has vanishing monodromy under translation of coordinate κ → κ + 2K + 2iK ′ , therefore the second part does not contribute to integral (93). Taken together, the nonvanishing contribution for integrals of v(χ) only come from terms in the first part of v 2ℓ+1 with negative power of cnκ. The monodromy relation (93) is simplified to wherec 2ℓ+1,2n+1 are another set of coefficients, they are polynomials of σ, ✵ m 3 and q. The integrals J −(2n+1) are computed from recurrence relations, as explained in Appendix B of Ref. [12].
The definite integral (95) gives the functional relation µ(σ), from which the eigenvalue expansion σ(µ) can be derived. The corresponding wave function is obtained from indefinite integral´v(κ)dκ.

Duality for strong coupling spectra
We have obtained two strong coupling spectral solutions, they are directly computed by perturbation method when the corresponding effective potentials are known. It turns out that the two sets of spectral data are related by simple transformations, which would be shown in the following. Regarding the connection to gauge theory, these transformations are the monopole-dyon duality of N f =4 super QCD.

Duality transformation for eigenvalues
Eigenvalue σ turns out to be related to eigenvalue δ by a transformation given by (96).
To explain this, now we treat parameters µ, q and Ω m in (65) as independent variables. That is to say, as the duality transformations discussed here are limited to spectral data, the coupling constants Ω m are regarded as free parameters. The explicit form of Ω m given in Appendix C, which are functions of q, is a particular case for Ω m determined by the connection to N f =4 gauge theory. Then the expansion of eigenvalue σ(µ, ✵ 2m , ✵ 2m+1 , q) is given by the expansion of eigenvalue δ(µ, Ω 2m , Ω 2m+1 , q) with the arguments transformed as where the coupling constants ✵ m for m 1 are given by the relations (89). Notice the disparity of the transformation rules for the even coupling constants Ω 2m and odd coupling constants Ω 2m+1 . This kind of relation between spectra at the magnetic and dyonic expansion regions is also verified for other examples of elliptic potential/super gauge theory models [11,12]. But in those simpler cases the transformation rule does not involve transformation for coupling constants Ω m .

Duality transformation for integrands of eigenfunctions
There should be a duality transformation relates the wave functions near the critical points η * 5 and η * 6 . The wave functions at the magnetic point is given by (66), the wave functions at the dyonic point is given by (D.5) in Appendix D. In this subsection we use notations ψ D ± (χ) and ψ T ± (κ) to distinguish the wave functions, use notations v D (χ) and v T (κ) to distinguish the corresponding integrands, at the magnetic and dyonic points respectively.
where ⌈ * ⌉ is the ceiling function, the coupling constants ✵ m are related to Ω m by (89). The integrands v D ℓ (or v T ℓ ) only contain Ω m (or ✵ m ) with m 3, the dominant one Ω 2 (or ✵ 2 ) is the expansion parameter for the total integrand v D (or v T ). The transformation (97) combined with Ω 2 → ✵ 2 transforms the total integrands v D → v T . In particular, the duality transformation for elliptic function is a simple substitution snχ → cnκ, this simplicity becomes obscured after performing integration.
In Appendix F we discuss a truncated version of the transformation that relates strong coupling solutions of the ellipsoidal wave equation.

Inverse duality transformations
The duality discussed above are transformations mapping the magnetic solution to the dyonic solution. Because the magnetic and dyonic points are really symmetric to each other, we could have presented the strong coupling solutions in a different order, started from the potential and its spectrum at the dyonic point and found the transformations mapping them to the potential and spectrum at the magnetic point.
We first derive the inverse map of the relation (89) between effective coupling constants.
Starting from the effective potential at the critical point η * 6 which is written in the κcoordinate as given in Eq. (88). Following an argument similar to that in subsection 4.5.1, the effective potential at the critical point η * 5 (or κ * 5 ) can be expressed in the same form when another local coordinate κ ′ is used, where η = η * 5 + κ ′ − K, . Near the point η * 5 we have κ ′ ∼ K, therefore cnκ ′ ≪ 1 is satisfied, the potential is an asymptotic series in cnκ ′ .
To rewrite the potential (98) in the χ-coordinate, notice that by definition η = η * 5 + χ, so that κ ′ = χ + K. Using elliptic function relations given in (86) with upper sign the potential (98) becomes an expression involving negative integer powers of dnχ. Then the following substitution rule is used to convert it to an expression in snχ, The region near η * 5 is the same region near χ ∼ 0 in the χ-coordinate where snχ ≪ 1. The potential, when expanded in snχ, is an infinite series given by (57). The effective coupling constants Ω m are related to ✵ ′ m by the inverse transformations , In fact, they can be obtained from transformation (89) by substitutions Ω ℓ (±q 2 ) and q → − q 1−q . The inverse duality transformations mapping the dyonic spectral data to the magnetic spectral data are similar to the duality transformations just presented. The inverse transformation between eigenvalues δ(µ, Ω 2m , Ω 2m+1 , q) and σ(µ, ✵ 2m , ✵ 2m+1 , q) is The inverse transformation between integrands The coupling constants Ω m in (101) and (102) should be substituted by the relations (100).
The transformation (102) should combine ✵ 2 → Ω 2 to transform the total integrand as v T → v D .

Strong coupling expansion of gauge theory prepotential
The spectral problem studied so far is related to super QCD with ǫ-deformation, all quantities obtained can be expanded in ǫ. In this section we compute the dual prepotential of N f = 4 super QCD without ǫ-deformation. The computation can be extended to deformed super QCD by including higher order ǫ-expansion, as explained in Appendix G. A simpler example to test the strong coupling expansion of gauge theory is the N= 2 * super QCD.

Magnetic expansion of prepotential for N f =4 super QCD
The magnetic dual prepotential F D (a D ) is computed by the relation between scalar v.e.v a and its magnetic dual a D , where a and a D are given by integral of Seiberg-Witten differential form along conjugate cycles α and β on the curve of N f =4 super QCD [8,9]. With the connection to spectral problem, they can be computed from spectral data of the DTV potential. As the undeformed super QCD is the ǫ → 0 limit of the deformed theory, the leading order part in the ǫ-expansion of v(χ) is needed to compute a and a D . According to the discussion in subsection 4.4 they are given by and In the square root eigenvalue δ and coupling constants Ω m are the leading order part of their ǫ-expansion, the subscriptions are dropped for simplicity. The dual integral a D in (105) is related to the Floquet exponent µ of the first strong coupling spectral solution (76) by So the integral in (105) is the classical limit of the monodromy relation (62). The corresponding classical limit of eigenvalue δ (0) (µ) is given in (79), substituting the coefficients Ω m expressed in terms of flavor masses of super QCD, it is This expression reproduces the result computed from the normal form of Heun equation presented in Ref. [10], where the eigenvalue is characterized by another quantity ∆. The variableâ D in Ref. [10] is computed by integral along a cycle on the elliptic curve associated to the normal form of Heun equation, it is twice of the variableâ D used in this paper. Taking into account the factor two discrepancy forâ D , then δ is related to ∆ by Then we need to perform the integral in (104) to get the functional relation a(δ), substitute the eigenvalue δ(µ) to get the functional relation a(a D ), then integrate the relation (103) to obtain the magnetic dual expansion of prepotential. In this section the integral in (104) is computed in the region of strong effective coupling, so we should not bring here the result obtained in Section 3 which is the weak coupling solution. Neither we should use Ω 1 2 2 as a large expansion parameter, because in the expansions (74) and (75) the series P β (0);2ℓ+1 (χ) is used for integral over period 2iK ′ , they can not be used for integral over period 2K.
In the effective elliptic potential, terms with coefficients Ω m 3 are treated as perturbation to the dominant term Ω 2 sn 2 χ. This consideration leads to the correct expansion In the numerator of right side, power functions of cnχ and dnχ are simplified using the relations in (56). The even coefficient series actually vanish, P α (0);2ℓ (κ) = 0. The reason we keep them in the expansion is that when generalising the computation to include ǫ-deformation, the even coefficients are non-vanishing at higher order ǫ-expansion, see Appendix G. The odd coefficient series P α (0);2ℓ+1 (χ) are expanded in snχ, In the integrand given by (109) and (110) the elliptic functions appearing in the numerator P α (0);2ℓ+1 (χ) are either sn 2m+1 χ cnχ dnχ or sn 2m χ, they lead to two classes of integrals over period 2K. The first class integrals actually vanish, with 2m + 1 3(ℓ + 1) 0. The second class integrals give non-vanishing results, with 2m 3(ℓ + 1) 0. An easier way to see the vanishing or non-vanishing of these integrals is by rewriting the integrals over periods as contour integrals in the ξ-coordinate defined by ξ = sn 2 χ. The details are explained in Appendix H.
Therefore, for practical purpose terms with sn 2m+1 χ cnχ dnχ in the series P α (0);2ℓ+1 (χ) can be dropped. The resulting series, denoted by P α (0);2ℓ+1 (χ), contain only terms with sn 2m χ. Rewriting elliptic functions in terms of amplitude ϕ = amχ, the non-vanishing integrals in (109) becomê where x 2 = −Ω 2 /δ, with magnitude ||x 2 || ≫ 1. The numerator in (113) can be expanded as a series in sin 2 ϕ. Then computation of the integral in (104) boils down to computation of the integrals in (113), with indices satisfying 2m 3(ℓ + 1). Their values are given by complete elliptic integrals with elliptic modulus x 2 . The details are explained in Appendix I where recurrence relations of M 2ℓ+1,2m are given. Following the procedure from (109) to (113), the function a (0) (δ) is expanded in variable x as The coefficients CE 2ℓ+1 (Ω m 2 , q) and CK 2ℓ+1 (Ω m 2 , q) can be expanded as large Ω 2 -series, a few of them are given in Appendix G. Using the large-x expansion of complete elliptic integrals E(x 2 ) and K(x 2 ), substituting the relation δ(a D(0) ) given by (79), we get the relation a D(0) ). The magnetic dual expansion of prepotential is determined from the relation (103) by an integration with respect to a D(0) , in terms of coupling constants Ω m it is and ℓ i referring to the exponent of higher order coupling constants Ω m i 3 . See Appendix G for the first few of them needed to produce (117) given below. Notice that there is only one logarithm term in the prepotential, this fact depends on precise cancellations between many logarithm terms produced by the large-x expansion of the complete elliptic integrals.
Substituting the expressions of Ω m given in Appendix C, the first few orders of prepotential are The expansions of eigenvalue (107) and magnetic dual prepotential (117) are related by the relation taking into account the relationâ D = 1 2 µǫ. Or in terms of the quantity ∆ used in Ref.
[10] the relation is This is the magnetic dual version of the relation in N=2 gauge theories connecting moduli and prepotential [45].

Dyonic expansion of prepotential for N f =4 super QCD
The dyonic dual prepotential is defined by the relation where a T is the dyonic dual of scalar v.e.v, it is given by integral of Seiberg-Witten differential form along the γ-cycle on the curve of N f = 4 gauge theory. In its connection to spectral solution of the DTV potential, a and a T are given by integrals of v(κ) associated to elliptic potential at the critical point η * 6 . At the leading order of ǫ-expansion they are given by and where u(κ) is the effective potential given in (88). In the square root eigenvalue σ and coupling constants ✵ m are the leading order part of their ǫ-expansion.
The dual integral a T is related to the Floquet exponent of the second strong coupling spectral solution by So the integral in (122) is the classical limit of the monodromy relation (93). The classical limit of the corresponding eigenvalue σ (0) can be computed following a similar procedure as in subsection 4.4.
There is a simpler ways to obtain σ (0) , by applying the duality transformation (96) to the classical eigenvalue δ (0) given in Eq. (79). Up to order O(µ 2 ) the classical limit of the eigenvalue σ is Using the effective coupling constants ✵ m expressed in terms of flavor masses of super QCD, computed from Ω m using the relation (89), it is The computation of integral in (121) follows a similar process as in the magnetic dual case. In the dyonic expansion region the effective coupling constant ✵ 2 is the large parameter, ✵ 1 2 2 ≫ σ. However, simply performing the large ✵ 2 expansion for the integrand σ + u(κ) leads to an expression incompatible with integral over period 2K. Similar to the magnetic case, the correct expansion form is The even coefficient series vanish, Q α (0);2ℓ (κ) = 0. The odd coefficient series Q α (0);2ℓ+1 (κ) are expanded in cnκ, they are counterpart of the series P α (0);2ℓ+1 (χ). The first few are Because elliptic functions of Q α (0);2ℓ+1 (κ) are in the forms cn 2m+1 κ snκ dnκ and cn 2m κ, there are also two classes of integrals over 2K follow from (126) and (127). The first class integrals give vanishing value, with 2m + 1 3(ℓ + 1) 0. The second class integrals are nonzero, with 2m 3(ℓ + 1) 0. This fact is also explained in Appendix H using contour integrals in the ξ-coordinate defined by ξ = sn 2 κ. So terms involving cn 2m+1 κ snκ dnκ in the expansions of Q α (0);2ℓ+1 (κ) can be dropped, the resulting series are denoted by Q α (0);2ℓ+1 (κ) which contain only terms with cn 2m κ. To perform the remaining integrals, rewriting elliptic functions in terms of the amplitude ϕ = amκ, the non-vanishing integrals in (126) becomê where in this case x 2 = ✵ 2 /(✵ 2 +σ). As Q α (0);2ℓ+1 becomes a series in sin 2 ϕ by the substitution cn 2m κ = (1 − sin ϕ 2 ) m , the numerator in (130) is an infinite series in sin 2 ϕ. Therefore, in order to compute a (0) we need the same integrals as in (114), now with index satisfying m 0, and the elliptic modulus satisfying ||x 2 || ∼ 1. See more details in Appendix I. The procedure from (127) to (130) leads to the function a (0) (σ) which is expanded in variable x as The coefficients CE ℓ (✵ m 2 , q) and CK ℓ+1 (✵ m 2 , q) can be expanded as large ✵ 2 -series, see Appendix G. Substituting the expansion of complete elliptic integrals at x ∼ 1 and the series of σ given in (124), we get the relation a (0) (a T (0) ). The dyonic dual expansion of prepotential is determined from the relation (120) by an integration with respect to a T (0) , in terms of coupling constants ✵ m it is The coefficients D n (q). The first few of them are given in Appendix G. The complete elliptic integrals expanded at x ∼ 1 produce many logarithm terms with ln(x − 1), they almost cancel when computing the prepotential, leaving a single logarithm term.
The coupling constants ✵ m in terms of super QCD masses can be obtained from the transformation (89) and the expansions of Ω m , substituting them into (133) we get There is also a relation between the expansions of eigenvalue and dyonic dual prepotential, On the left side µ is substituted by the relation (123), on the right side a T is treated as a variable independent of q.
And the inverse transformation for the dual expansions of prepotential is Taking various mass decoupling limits for the magnetic and dyonic expansions of prepotentials, they lead to corresponding dual expansions for prepotentials of super QCD with less hypermultiplets and pure gauge theory. There is another class of limits, discussed in Section 6, which do not decouple any hypermultiplet, provide a connection between the N f =4 super QCD and the N=2 * super QCD.

Limits of gauge theory and Landen transformations
In this section we discuss a particular class of limits for the N f =4 super QCD and the related spectral problem. In gauge theory these limits are achieved by setting the flavor masses to special values, the partition function of N f =4 super QCD model reduces to the partition function of N=2 * super QCD model. In the corresponding spectral problem, these limits are obtained by setting the coupling constants b s to special values, the DTV potential reduces to the Lamé potential by Landen transformation of elliptic functions.
There are three Landen transformations, the duplication formula is also included as it is similar to Landen transformations. Transformations for Jacobian elliptic function can be found in Section 22.7 in Ref.
[20], transformations for the shifted Weierstrass elliptic function are given in Appendix A.

Descending Landen transformation
Making a special choice for the coupling constants of the DTV potential, b 0 = b 1 = b, b 2 = b 3 = 0, the potential (5) reduces to the Lamé potential, by the descending Landen transformation The hat quantities are related to the un-hat ones bŷ The hat half periods are For the Lamé potential (139), the critical points are atη * e = K + i K ′ ,η * m = 0 andη * d = K, the corresponding values of the potential are u(η * e ) = 2 1 + (1 − q) 1 2 b, u(η * m ) = qb and u(η * d ) = 2 1 − (1 − q) 1 2 b, respectively [12]. Notice that when rewriting derivative ∂ 2 η to ∂ 2 η in the Schrödinger equation, the equation is rescaled by an overall factor 1 + (1 −q) 1 2 2 . The potential appearing in the Lamé equation This comment applies to all forms of Lamé potential obtained in this section by taking limits associated to Landen transformations.
We would like to know how the six critical points η * 1−6 evolve in the limit of Landen transformations. The resulting Lamé potential has three critical points. There are a few possible fates for η * 1−6 when b s are turned to limiting values. They might persist as critical points, and move to merge with each other to form a critical point of the Lamé potential; or they move position but do not merge; some critical points might disappear, and new critical points emerge.
In the case of descending Landen transformation, however, it is not straightforward to see how the critical points evolve. The series solution of η * 1 has definite limit for By the relations (141) and (142), the position 1 2 K + iK ′ in the η-coordinate is mapped to the position K + i K ′ in theη-coordinate, it means η * 1 moves toη * e . For the points η * 2−5 , in their series solutions (9)-(14) the denominators contain powers of (b , taking the limit directly would leads to indefinite results. The problem is caused by the fact that we have used q-series solutions for the critical points η * 1−6 , each of them has valid regions in the parameter space. The series solution of a critical point does not necessary have definite limit even the point has. The descending limit pushes mass parameters to a region which is still valid for series solution of η * 1 , but is outside the valid regions for series solutions of η * 2−5 . There are two possibilities for η * 2−5 in the descending limit, they might have a limit position but their series solutions fail to capture the limit, or they might have no limit because they disappear when parameters are approaching the limiting value. The gauge theory side provides some help to trace their evolution. In N f =4 super QCD model, up to the leading order of ǫ-expansion flavor masses µ i are determined by b s using relations b 0 = ǫ −2 (µ 1 − µ 2 ) 2 etc, the descending limit forces three flavor masses to vanish, such as µ 1 = b 1 2 ǫ and µ 2,3,4 = 0. As shown in Ref. [10], the critical points η * 1−4 are in one-to-one associated with flavor masses µ 1−4 , and the critical points η * 5,6 are associated respectively with strong coupling mass scales ±4(µ 1 µ 2 µ 3 µ 4 ) 1 4 q 1 4 . When a mass µ i vanishes, the strong coupling mass scales also vanish. On the other hand, for the reduced potential (139) which is related to the N=2 * super QCD, the critical pointη * e is associated with the non-vanishing flavor mass µ 1 , the critical pointsη * m ,η * d are associated with the strong coupling mass scales ±ǫ(bq) 1 2 ∼ ±µ 1 q. Comparing the nature of initial and final critical points, we conjecture that in the descending limit b 0 = b 1 = b, b 2 = b 3 = 0 the behavior of the critical points is, η * 1 →η * e , and η * 2−6 disappear, two new critical pointsη * m andη * d emerge.
When the coupling constants take values b 0 = b 1 = 0, b 2 = b 3 = b, notice that translation of coordinate η → η + iK ′ for the DTV potential leads to interchange of coupling constants b 0 ↔ b 2 , b 1 ↔ b 3 , therefore, this limit follows the same discussion of the previous limit for potential u(η + iK ′ ). The DTV potential reduces to the Lamé potential with critical points atη * e = K,η * m = i K ′ andη * d = K + i K ′ . In this case one critical point has definite limit, so now it is η * 2 that moves toη * e . Similar to the previous case, taking the limit directly for series solutions of η * 1 and η * 3−6 given in Eqs. (9)-(14) leads to indefinite results. We conjecture that in the descending limit b 0 = b 1 = 0, b 2 = b 3 = b, the critical points evolve as η * 2 →η * e , and η * 1,3−6 disappear, two new critical pointsη * m andη * d emerge.

Ascending Landen transformation
It is also called Gauss' transformation in Ref. [46]. It corresponds to making a choice for the coupling constant using the ascending Landen transformation q sn 2 η + 1 The tilde quantities are related to the un-tilde ones bỹ η = (1 + q 1 2 )η,q = 4q The tilde half periods are There are three critical points for the Lamé potential (147) atη * e = K,η * m = i K ′ and η * d = K + i K ′ where the values of the potential are u(η * e ) = (1 + q)b, u(η * m ) = −2q The series solutions of critical points (9)- (14) are valid in the slice of ascending limit in the parameter space, so the limits of all critical points have definite value. Their limits are In the η-coordinate the critical points move as η * 1,2 → K + iK ′ , η * 3,4 → K, η * 5 → K + 1 2 iK ′ and η * 6 → 1 2 iK ′ . In theη-coordinate these limit positions are located at K+2i K ′ , K, K+i K ′ and i K ′ , respectively. The first two positions are actually the same up to the period 2i K ′ .

Irrational Landen transformation
It is also called generalized transformation in Ref. [20]. In this case the choice for the by the irrational Landen transformation The bar quantities are related to the un-bar ones bȳ The bar half periods are The critical points of the potential (154) are atη * e = 0,η * m = K + iK ′ andη * d = iK ′ , and the corresponding values of the potential are u(η * e ) = b, u(η * m ) = −2iq The series expansions of critical points are also valid in the slice of irrational limit in the parameter space, hence have definite limits, In the η-coordinate the critical points move as η * 1,2 → K − iK ′ , η * 3,4 → 0, η * 5 → 1 2 (K + iK ′ ) and η * 6 → 1 2 (K − iK ′ ). In theη-coordinate these limit positions are at −2iK ′ , 0, iK and K − iK ′ , respectively. The first two positions are the same, the fourth one equals toη * m , up to the period 2iK ′ . We conclude in this limit the critical points evolve as η * 1,2,3,4 →η * e , η * 5 →η * d , There is also another equivalent choice for the coupling constants b 0 = b 3 = 0, b 1 = b 2 = b, using the same translation of coordinate η → η + K to interchange coupling constants as with critical points atη * e = K,η * m = iK ′ andη * d = K + iK ′ . In the η-coordinate the critical points move as η * 1,2 → iK ′ , η * 3,4 → K, η * 5 → 1 2 (K − iK ′ ) and η * 6 → 1 2 (K + iK ′ ). In thē η-coordinate these limit positions are at −K + 2iK ′ , K, K − iK and iK ′ , respectively. The first two positions are the same, up to the period 2K − 2iK ′ ; and the third one equals tō η * d , up to the period 2iK ′ . The evolution of the critical points in this limit turns out the be the same as that in (159).
Inspecting the limit behavior of η * 1−6 in the limit of Landen transformations discussed above, there is a general phenomenon. For the weak coupling critical points η * 1,2,3,4 , when a flavor mass of the N f =4 super QCD is kept nonvanishing then the associated critical point persists, might moves to merge with other critical points of the same nature, while for those vanishing masses the associated critical points disappear. For the strong coupling critical points η * 5,6 , only when all four flavor masses are kept nonvanishing then they persist, might exchange the magnetic/dyonic nature with each other, otherwise they disappear and two new strong coupling critical points of the N=2 * super QCD emerge elsewhere.
Because the coordinate η is scaled by a factor 2, the half periods becomes 1 2 K, 1 2 iK ′ and 1 2 K + 1 2 iK ′ . The critical points of the potential (161) are at η * e = 1 2 K, η * m = 1 2 iK ′ and η * d = 1 2 K + 1 2 iK ′ where u(η * e ) = 4b, u(η * m ) = 0 and u(η * d ) = 4bq, respectively. In the duplication limit, there are two series expansions for η * 1 and η * 3 having definite limits, Because iK ′ is a period of the scaled coordinate 2η, so points η * 1 and η * 3 correspond to the same point η * e . The series solutions of η * 2 and η * 4,5,6 do not have definite limit. For the N f =4 super QCD model, in the duplication limit two flavor masses vanish while the other two flavor masses remain nonzero, so the strong coupling mass scales vanish accordingly.
Therefore, similar to the case of descending limit, we conjecture in the duplication limit two critical points η * 1 and η * 3 merge to a new critical point η * e , other critical points η 2,4−6 disappear, and two new critical points η * m and η * d emerge.

Landen transformations of effective potential u(χ)
The Landen transformations also apply to the effective potential u(χ) at the critical point η * 5 , but it manifests in a different way from discussion in the previous subsection.

Descending Landen transformation
We need to determine the limit value of Ω m . As discussed in subsection 6.1.1, the critical point η * 5 has no definite limit for descending transformation, consequently neither the coupling constants Ω m have. This is also explained using the expansion of Ω m given in Appendix C. The descending Landen transformation requires three flavors to have vanishing masses.
This is in conflict with the premise that the quantity (µ 1 µ 2 µ 3 µ 4 ) 1 2 q 1 2 is sufficiently large so the effective coupling constant Ω m have the right hierarchical structure to ensure validity of u(χ) near the critical point η * 5 . Therefore, the potential u(χ) has no descending Landen transformation.

Ascending Landen transformation
As discussed in subsection 6.1.2, the critical point η * 5 has a definite limit for ascending transformation. For the case b 0 = b 2 = b, b 1 = b 3 = 0, the limit of critical point is η * 5 → Then Ω m are determined by expanding the DTV potential nearη * 5 , following a similar procedure as in subsection 4.1, the effective potential is computed by For the case b 0 = b 2 = 0, b 1 = b 3 = b, the limit of critical point is η * 5 →η ′ * 5 = 1 2 iK ′ , the effective potential is computed by The two cases lead to the same limit for effective coupling constants, for m 1. The effective potential reduces to where the tilde quantities areχ (168)

Duplication formula
The It is again in conflict with the condition that (µ 1 µ 2 µ 3 µ 4 ) 1 2 q 1 2 is a large parameter for the potential u(χ). However, in this case a trick can be used to obtain definite limit. We first deform the limit to b 0 = b 2 = b + ε, b 1 = b 3 = b, then take the limit ε → 0. We get sn 2 η * 5 → q − 1 2 , the limit position of η * 5 is the same as the case of ascending transformation, η * 5 → K + 1 2 iK ′ . Or we may deform the limit for masses to µ 1 = µ 3 = b 1 2 ǫ, µ 2 = µ 4 = ε, so that Ω m given in Appendix C can be used. In this way, the limit of coupling constants are given by for m 1. The effective potential reduces to u(χ) = 4bq sn 2 (2χ, q). (175)

Landen transformations of effective potential u(κ)
The Landen transformations of the potential u(κ) at the critical point η * 6 are very close to the situation for the potential u(χ). We summarize the results.
There is no descending transformation for u(κ) because the limit leads to (µ 1 µ 2 µ 3 µ 4 ) 1 4 q 1 4 → 0, incompatible with the requirement that (µ 1 µ 2 µ 3 µ 4 ) 1 4 q 1 4 is sufficiently large. In the limit of ascending transformation, the limit of ✵ m are given by the limit of Ω m (−q 1 2 ) using the relations (89), and the limit of Ω m (−q 1 2 ) are obtained from (166) by a sign inversion q for m n 1. The effective potential becomes where the tilde quantities arẽ In the limit for irrational transformation, the limit of Ω m (−q for m l n 1. Then the effective potential is given by a triple sum where the bar quantities arē In the limit for duplication formula, the same trick as in the previous subsection is used to obtain definite limit for the coupling constants. The limit of Ω m (−q 1 2 ) are the same as the limit of Ω m (q 1 2 ) because the sign inversion q 1 2 → −q 1 2 does not affect integer powers of q. By the relation (89), the coupling constants ✵ m are given by where ⌊ * ⌋ is the floor function. The effective potential reduces to u(κ) = 4bq sn 2 (2κ, q).
According to the discussion in subsection 4.5.1 about the relations of coordinate systems, χ ′ = κ − K. Therefore, in these limit the potential u(κ) indeed reduces to a Lamé potential at the critical point η * 6 , as shown in Fig. 2.

Resolving a sign puzzle
The Landen transformations change the coordinates, elliptic nome/modulus and half periods of potentials u(η), u(χ) and u(κ) almost in the same way, except a sign mismatch. For example, in the case of ascending transformation, the tilde quantities given by (149) for u(η) and by (178) for u(κ) involve +q 1 2 , but the tilde quantities given by (168) for u(χ) involves −q 1 2 . The situation is similar for the irrational transformation. The sign mismatch for potentials u(χ) and u(κ) has a simple explanation. Taking the limits of ascending and irrational transformations, u(χ) and u(κ) reduce to a standard Lamé potential at points η * 5 and η * 6 , respectively. Because solutions of η * 5 and η * 6 are mapped to each other by the sign inversion q 1 2 → −q 1 2 , the resulting Lamé potentials are also mapped to each other by the sign inversion combined with χ → χ ′ for local coordinates.
However, the sign mismatch for potentials u(η) and u(χ) has a different reason. If we have chosen a different assignment for solutions of η * 5 and η * 6 from Eqs. (13), say the − sector solution for η * 5 and the + sector solution for η * 6 , then the sign mismatch for u(η) and u(χ) can be avoided, however, a sign mismatch for u(η) and u(κ) would have arisen.
Therefore, the sign inversion q 1 2 → −q 1 2 for effective potentials at η * 5 and η * 6 only exchanges where the sign mismatch shows up, it could not really eliminate it.
In the following we show the puzzle is resolved by certain delicate properties of elliptic functions. Taking the limit b 0 = b 2 = b, b 1 = b 3 = 0 for ascending Landen transformation as an example. The critical point η * 5 in this limit is determined by sn 2η * 5 = q − 1 2 , so it locates atη * 5 = K(q) + 1 2 iK ′ (q). Then the coordinates are related by η = χ + K(q) + 1 2 iK ′ (q). Now we use subscripts ± to distinguish the tilde quantities appearing in subsections 6.1.2 and 6.2.2 respectively. The tilde variables relevant for discussion in 6.1.2 are given by The tilde variables relevant for discussion in 6.2.2 are given by The variables χ,χ + do not appear in subsection 6.1.2 and the variables η,η − do not appear in subsection 6.2.2, nevertheless, they are needed for discussion in the following. The + and − sectors of the tilde variables are conjugate to each other in the following sense, and By definition of the + sector variables we haveη + =χ + + K + + i K ′ + . With the relations given above, it can be shown that the Lamé potential (147) and (153) are in fact the same as the Lamé potential (167), after removing a superficial sign difference.
In order to relate the potentials (147) and (167), we need the following relation, where in the last step a relation about Jacobian elliptic functions with conjugate coordinates and conjugate elliptic modulus is used. Therefore the potential (147) in fact is identical to the potential (167) up to a constant, In the limit b 0 = b 2 = 0, b 1 = b 3 = b, the critical point η * 5 moves toη ′ * 5 = 1 2 iK ′ , so the coordinates are related by η = χ + 1 2 iK ′ (q). The potentials (153) is also identical to the potential (167).
It can be shown in a similar way that for the case of irrational transformation the potentials (154) and (160) are actually the same as the potential (172).

Summary and discussion
We have examined the relations between a few SU(2) supersymmetric gauge theory models and Schrödinger equation with periodic potentials, the case studied in this paper is the richest one of this class. Various features of gauge theory have a corresponding interpretation in terms of spectral solutions of periodic potential. There is one aspect of this connection that is still missing, at the moment there are no gauge theory quantities that correspond to the strong coupling wave functions of the quantum mechanical models. In the weak coupling region, the instanton partition function with surface defect is related to the weak coupling wave function. Probably there are partition functions of effective gauge theory with defect at monopole and dyon singularities in the moduli space which are related to the strong coupling wave functions.
All the spectral solutions can be verified by plugging eigenvalue and wave function into the corresponding equation. Another way to verify is by checking a relation for second-order differential equations. For example, for the weak coupling solution obtained in subsection 3.1, taking the quotient function f (x) = ψ ± (x)/ψ ∓ (x) with either − or + sector wave function, then the following relation holds, where S(f (x)) is the Schwarzian derivative of f (x). The strong coupling spectra obtained in Section 4 also satisfy similar relations.
The effective potentials (57) and (88), derived by expanding the DTV potential near its magnetic and dyonic critical points, are very generic elliptic Hill potentials compatible with doubly periodic Floquet property, generalising other widely used periodic and elliptic potentials. The classical Floquet theory of periodic potential remains valid for elliptic potentials, however, for strong coupling solutions it only manifests clearly when the wave functions are written in canonical coordinates. A question it raises is, whether the DTV potential, or in the guise of effective potentials (57) and (88), are the most generic elliptic potentials of this kind.
The example studied in this paper is a non-relativistic two particles quantum mechanical problem with elliptic potential, having direct connection with supersymmetric gauge theory. It has a relativistic generalization, the corresponding spectral problem is a difference equation with coefficients given by elliptic functions [47,48,49]. Various relativistic or non-relativistic version of elliptic, trigonometric/hyperbolic, rational potential two particles quantum mechanical systems are obtained from limits of the relativistic-DTV potential. In the context of Gauge/Bethe correspondence, relativistic generalization of integrable system corresponds to growing a compact spacetime dimension for gauge theory [5,42,43]. Relativistic generalization of Toda class integrable models are related to 5-dimensional N=1 supersymmetric gauge theory compactificated on a circle, the partition function of gauge theory is computed by K-theoretic localization formula [33,41]. Some recent discussion of this connection are in Refs. [50,51,52,53]. It might be helpful to study the relativistic-DTV potential from the perspective of gauge theory, a version of relativistic Heun/van Diejen operator indeed shows up in supersymmetric field theory [54]. Nevertheless, we should be aware that making their connection precise might be not straightforward, variety of twists could arise as we have seen in the relation between non-relativistic modes with elliptic potentials and super QCD models (N f =4 and N=2 * ) studied so far.

A Useful relations of elliptic functions
The following formulae accompany some other relations listed in Ref. [12].

A.1 Simplifying Weierstrass elliptic function
In the course of computation weak coupling spectrum using potential in the Weierstrass form, the Hamiltonian densities of KdV hierarchy v ℓ (x) contain various products of shifted Weierstrass elliptic functions ℘(x + ω s ) and their derivatives ∂ n x ℘(x + ω s ), they are linear combination of two types of terms. The "P(ower)-type" terms are In the following procedure, the factors of coupling constants b ms+ns s are kept unchanged, therefore omitted in the discussion. The first step is to decrease d s and n s of D-type terms in (A.2). The higher order derivatives [∂ d x ℘(x)] n are simplified using relations such as In this section we briefly use c m ′ , etc, to denote constant coefficients. For terms with odd d, a derivative factor [∂ x ℘(x)] n still remains. Using the basic equation (A.3), an even exponent can be decreased to n = 0, an odd exponent can be decreased to n = 1, so a further simplification follows, Plugging the factor [∂ ds x ℘(x + ω s )] ns into the product in (A.2), it leaves multiple terms taking the general form with n ′ s = 0 or 1, order of poles of v ℓ (x) require the conditions 2m ′ s + 3n ′ s ℓ + 1. If in a monomial all n ′ s = 0, there is no derivative factor anymore, it is a P-type term. Otherwise, it is still a D-type term with a derivative factor which can be simplified further using with indices i, j, k ∈ {1, 2, 3} and i = j = k. For s n ′ s = 2 or 4, the derivative factor is totally dissolved; for s n ′ s = 1 or 3, a single derivative ∂ x ℘(x + ω s ) is left. Now the remaining D-type terms are in the form with t ∈ {0, 1, 2, 3}. Because the periods ω s are symmetric, there must be other three terms joining up to form a symmetric piece which is a total derivative of the form with 2m ′′ s ℓ. In summary, the first step simplification transforms part of D-type terms into P-type terms, the remaining part of the D-type terms are derivative of P-type terms.
The second step is to break products and decrease the exponent m s in P-type terms, including the P-type terms in (A.9), until no products of function ℘ is left. When there are products of functions with different period shifts ℘ ms (x + ω s ) ℘ mt (x + ω t ), repeatedly use the following relations to break products and decrease indices m s and m t until the smaller one reaches zero, (A.10) It leaves functions with the same period shift ℘ ms (x + ω s ) with m s 0. All power functions ℘ ms (x + ω s ) are traded by functions with total derivative terms, The second step simplification transforms P-type terms to linear combinations of constants, ℘(x+ ω s ) and even order derivative terms ∂ 2ds x ℘(x+ ω s ) with d s 1. Then the D-type terms, according to (A.9), are linear combinations of odd order derivative terms ∂ 2ds+1 x ℘(x + ω s ) with d s 0. The order of derivatives rise again, but there is no products anymore, this is exactly what needed for integration of v ℓ (x).
In the end, in term of zeta function, v ℓ (x) are simplified to a form ready for computing eigenvalue and eigenfunction, The upper bound of d comes from order of poles for v ℓ (x), the coefficients c ℓ,s,d are related, but not equal, to the coefficients c l,s,d in the wave function (21) despite denoted by similar symbols.
Finally, in terms of elliptic theta function Landen transformations are See Section 20.7 in Ref.
[20] and Ref. [55]. Notice that, the nome of ϑ r (z, − √ p) in the third transformation of (A.17) looks different from the nome used in Refs.

B Identifying coordinates of weak coupling wave function
Instead of directly searching the relation between x 1 , x 2 and x, we use an intermediary variable to analyse the series expansion G −1 . A trial choice is x 1 = q 1 2 e −i2φ , x 2 = q 1 2 e i2φ . Then G −1 given in (34) is rephrased to There is a relation in elliptic function theory getting close to, but not exactly as, this expression. If ϕ is the amplitude of elliptic function sn(η|q), then the variables η and ϕ are related by the elliptic integral of the first kind, η = F (ϕ, q). However, the q-series of F (ϕ, q) does not have the right coefficients to match G −1 . Not quite obvious, but it can be recognized that a further twist on F (ϕ, q) could get all coefficients correct. Using the tilde variables (η,φ,q) resulting from the ascending Landen transformation of un-tilde variables (η, ϕ, q), given in Section 6, then G −1 is precisely given by whereη is related toφ by sn(η|q) = sinφ. To verify the relation using the following integral, performing the q-series expansion of integrand before integrating, and K = K(q) = F ( π 2 ,q). It looks odd to use the tilde variables for prepotential G. We can transform back to the un-tilde variables. The un-tilde amplitude ϕ is related to the tilde amplitudeφ by relatioñ ϕ = arcsin (1 + q (B.5) Now it can be verified G −1 is given by the elliptic integral of the first kind with un-tilde elliptic variables, where sn(η|q) = sin ϕ and K = K(q). The q-series of η is computed from From the relation between coordinates η and x given in (6), we have G −1 = iπ(x − ω 2 )/ω 1 . We are finally arriving at the functional relations x 1 (x) and x 2 (x). The first guess x 1 (φ) and x 2 (φ) are q 1 2 (cos 2φ ∓ i sin 2φ) = q 1 2 ( cn 2η − sn 2η ∓ 2i snη cnη). (B.8) Using the ascending Landen transformation, they are transformed to un-tilde variables By the relation of sn-function and theta function, the relation of coordinates η and x given in (6), they equal −q . (B.10)

C Effective coupling constants Ω m in terms of super QCD masses
The effective coupling constants Ω m and ✵ m are rational functions of b 1 2 s and q, in turn they are rational functions of µ 1 2 i , ǫ and q. When ǫ is treated as a quantization parameter, the "classical" limit of b s are b 0 ∼ ǫ −2 (µ 1 − µ 2 ) 2 , etc. In Section 5 the strong coupling expansions of prepotential are computed in the leading order of ǫ-expansion, for that purpose in the following we give q-series of the effective coupling constants Ω m up to m = 6 in the leading order of ǫ-expansion.
In the following expressions indices i, j, k ∈ {1, 2, 3, 4}. The summations restricted by i < j or i < j < k are customarily used. The summations restricted by i|j < k means the index i takes a value in {1, 2, 3, 4}, meanwhile indices j, k take values in the remaining set and satisfy j < k, that is The vacuum energy at η * 5 is The effective coupling constants are expanded as q-series whose coefficients are functions of masses. They are given by The coupling constants ✵ m are obtained from Eq. (89). It can be shown that ✵ m and Ω m with the same m-index are of the same order by counting the leading power of q. D More about the spectrum at η * 6 The second strong coupling solution of the DTV potential is associated with the critical point η * 6 . The locally valid effective potential u(κ) is used to compute the spectral solution.
Performing the indefinite integral´v(κ)dκ and substituting the eigenvalue (D.1) for the variable σ, the corresponding unnormalized wave functions are obtained, For solution at the dyonic critical point, the trajectory of integral over 2K + 2iK ′ in the κ-plane is chosen to cross the branch cut of ln{[ dnκ+(1−q) 1 2 snκ]/ cnκ}, give a monodromy −iπ, but avoid the branch cut of ln( dnκ +iq 1 2 snκ). See the discussion in Appendix B of Ref. [12]. The remaining terms in the wave function contain ln cnκ, cn m κ, sn 2m κ and snκ dnκ, with m ∈ Z + , are invariant under translation κ → κ+2K+2iK ′ . Therefore, the monodromy of wave function comes from a single term ±(1 − q) − 1 2 µ ln{[ dnκ + (1 − q) 1 2 snκ]/ cnκ}. In this case we could define the canonical coordinate θ conjugating to µ by 2 so that all terms of elliptic functions in the wave function can be expressed in terms of sinh θ and cosh θ. Then the wave functions take the Floquet form where p ± 3 (µ, θ) are periodic functions invariant under translation θ → θ + iπ, related by p + 3 (µ, θ) = p − 3 (−µ, θ). See more details in the next Appendix E. Now the monodromy relation (93) can be explained by a reasoning similar to that in subsection 4.3. As µ is conjugate to the canonical coordinate θ, the correct monodromy relation is ± iµπ (D.8) 2 We could define the canonical coordinate by then the wave function is written in a more familiar form as ψ ± (κ) = exp(±iµθ)p ± 3 (µ, θ).
Or in terms of κ-coordinate it is ± iµπ Then elliptic functions of χ are expressed by hyperbolic functions of ρ as snχ = 1 The periodic functions p ± 2 in Eq. (68) become Some constants arise, they are absorbed into the normalization factor of wave functions.

E.2 Wave function at η * 6
For wave functions of the dyonic strong coupling solution, by definition of the canonical coordinate (D.6) we get relations dnκ + (1 − q) [20]), is another extensively studied example. Both its weak and strong coupling asymptotic spectra have been known [56,12]. But it was considered not directly related to N=2 gauge theory. In fact, the ellipsoidal potential can be viewed as a truncated version of the elliptic Hill potential (57) obtained in subsection 4.1, by keeping Ω 2 , Ω 4 nonzero and setting all other coupling constants zero.
Rewriting the ellipsoidal potential in terms of function cnχ, it becomes a truncated version of the elliptic Hill potential (88) obtained in subsection 4.5.1. This fact brings closer the ellipsoidal potential and N=2 supersymmetric gauge theory.
In the following, we show that there is a duality relation between the two strong coupling asymptotic solutions of ellipsoidal wave equation. It is a truncated version of the duality For the wave function, the integrands of dyonic solution v T ℓ ( cnχ, σ, Ω 4 , q) can be obtained from the integrands of magnetic solution v D ℓ ( snχ, δ, q, Ω 4 , q) by transformation The inverse transformations map dyonic solution to magnetic solution, noticing that Ω 2 = −✵ 2 − 2✵ 4 and Ω 4 = ✵ 4 , they are and G Some details of computing the dual prepotential

G.2 Going to higher orders of ǫ-expansion
To include the ǫ-deformation, we need to find the integrand v(χ), compatible with integration over period 2K, from the equation v 2 (χ)+∂ χ v(χ) = u(χ)+δ. The computation in the leading order indicates it is crucial to separate δ + Ω 2 sn 2 χ from higher order coupling terms, and use δ + Ω 2 sn 2 χ as the expansion quantity for v(χ). The following expansion for v(χ) is assumed, (G.6) The coefficient series P α L (χ) contain only higher order coupling constants Ω m 3 . Substituting into the equation, the derivative ∂ χ v(χ) produces terms like − Ω 2 LP α L (χ) snχ cnχ dnχ (δ + Ω 2 sn 2 χ) containing Ω 2 in the numerator. This contradicts with the assumed form of the numerator.
The coupling constant Ω 2 is traded by then a term as above becomes This substitution allows v(χ) to be solved order by order.
In accordance with the potential u(χ), any higher order powers of cnχ and dnχ in P α L (χ) are further simplified using the substitution rule (56). The procedure leads to a solution, the first few coefficients are Notice that P α L (χ) are divided up into different parts according to their order by counting the power of ǫ. The leading order parts in both P α 2ℓ−1 (χ) and P α 2ℓ (χ) with ℓ 0 are of order ǫ −2ℓ , terms in brackets {· · · } are of higher order ǫ −2ℓ+2 , ǫ −2ℓ+4 , etc. Now the integrals over period 2K are the same form as those in Eqs. (111) and (112) used in Section 5.1, but with different range for the indices. Denoting them by S L,2m+1 and S L,2m , the indices L and m take values according to: 1. L take value of all integers L −1; 2. When L = 2ℓ, there are S 2ℓ,2m+1 with m −ℓ − 1 and S 2ℓ,2m with m −ℓ + 2; 3. When L = 2ℓ + 1, there are S 2ℓ+1,2m+1 with m −ℓ + 1 and S 2ℓ+1,2m with m −ℓ − 1.
The inequalities for m-index are saturated when L 1.
The necessary integrals for computation of dual expansions of ǫ-deformed gauge theory are explained in Appendix H and Appendix I.

H More prescriptions of contour integrals
When computing the spectra of elliptic Hill potentials (57) and (88) in Section 4, we need integrals I −(2n+1) over periods 2iK ′ and J −(2n+1) over periods 2K + 2iK ′ . Integrals over periods in the χ-coordinate are transformed to contour integrals in the ξ-coordinate through a coordinate change sn 2 χ = ξ, integrals over periods in the κ-coordinate are transformed similarly through sn 2 κ = ξ. From the first few ones I −1 , I −3 and J −1 , J −3 , other integrals are computed by recurrence relations, as explained in Ref. [12].
When computing the strong coupling expansion of gauge theory prepotential, we have encountered some other integrals of Jacobian elliptic functions over period 2K. In the following, we explain all integrals needed to compute the ǫ-deformed gauge theory by transforming periodic integrals in the χ-coordinate (or the κ-coordinate) to contour integrals in the ξcoordinate. with magnitude || δ Ω 2 || ≪ 1. The α-contour is chosen to encircle points ξ = 0 and 1, but avoid the point ξ = − δ Ω 2 . There are four possibilities:

H.1 Contour integrals for magnetic dual expansion
3. When L = 2ℓ + 1 −1 and m 0, the integrand has a branch cut [− δ Ω 2 , ∞) in the ξ-plane. The contour can be arranged to avoid the branch cut so that the integral vanishes, S 2ℓ+1,2m+1 = 0, as shown in Figure 3 Points at ξ = 0, 1, 1 q are branch points; the point ξ = − δ Ω 2 is a pole when L is even, a branch point when L is odd. In this case there are two possibilities: 1. When L = 2ℓ 0, the branch cuts can be arranged as [0, ∞) and [1, 1 q ], the contour cross both branch cuts so the integral is nonzero. Rewriting elliptic functions in terms of amplitude ϕ as in Section 5.1, 2. When L = 2ℓ + 1 −1, the branch cuts can be arranged as [0, − δ Ω 2 ] and [1, 1 q ] as shown in Figure 4, the contour cross both branch cuts so the integral is also nonzero.