Rigorous Asymptotics for the Lam\'e and Mathieu Functions and their Respective Eigenvalues with a Large Parameter

By application of the theory for second-order linear differential equations with two turning points developed in [Olver F.W.J., Philos. Trans. Roy. Soc. London Ser. A 278 (1975), 137-174], uniform asymptotic approximations are obtained in the first part of this paper for the Lam\'e and Mathieu functions with a large real parameter. These approximations are expressed in terms of parabolic cylinder functions, and are uniformly valid in their respective real open intervals. In all cases explicit bounds are supplied for the error terms associated with the approximations. Approximations are also obtained for the large order behaviour for the respective eigenvalues. We restrict ourselves to a two term uniform approximation. Theoretically more terms in these approximations could be computed, but the coefficients would be very complicated. In the second part of this paper we use a simplified method to obtain uniform asymptotic expansions for these functions. The coefficients are just polynomials and satisfy simple recurrence relations. The price to pay is that these asymptotic expansions hold only in a shrinking interval as their respective parameters become large; this interval however encapsulates all the interesting oscillatory behaviour of the functions. This simplified method also gives many terms in asymptotic expansions for these eigenvalues, derived simultaneously with the coefficients in the function expansions. We provide rigorous realistic error bounds for the function expansions when truncated and order estimates for the error when the eigenvalue expansions are truncated. With this paper we confirm that many of the formal results in the literature are correct.


Introduction
The Lamé equation is Lamé's equation first appeared in a paper by Gabriel Lamé in 1837 [14]. It appears in the method of separation of variables applied to the Laplace equation in elliptic coordinates. Lamé functions have applications in antenna research, occur when studying bifurcations in chaotic Hamiltonian systems, and in the theory of Bose-Einstein condensates, to name a few (see [24, § 29.19]). Mathieu's equation is where λ and h are real parameters. We consider the interval z ∈ (0, π). When λ assumes the special values a m or b m+1 for m = 0, 1, 2, . . . , Mathieu's equation admits even or odd periodic solutions denoted ce m (h, z) or se m+1 (h, z) respectively. These functions first arose in physical applications in 1868 inÉmile Mathieu's study of vibrations in an elliptic drum [15]. Since they have appeared in problems pertaining to vibrational systems, electrical and thermal diffusion, electromagnetic wave guides, elliptical cylinders in viscous fluids, and diffraction of sound and electromagnetic waves, to name a few. In general, they appear when studying solutions of differential equations that are separable in elliptic cylindrical coordinates. For an insight to see how Mathieu functions appear in physical applications see [16].
We wish to obtain uniform asymptotic approximations for the Lamé and Mathieu functions, and asymptotic expansions for their respective eigenvalues, as the parameters ν in Lamé's equation and h in Mathieu's equation become large. We denote for the moment the parameter h in Lamé's equation to be h L , to avoid confusion with the parameter h in Mathieu's equation. We have in the limit k → 0 + from [24, § 22.5(ii)] that lim k→0 + sn(z, k) = sin z, and lim k→0 + thus if ν → ∞ in such a way that as k → 0 + , ν(ν + 1)k = 2h for some constant h, then we can rewrite the limit of Lamé's equation in the form Thus for ν = − 1 2 + 1 4 + 2h k 2 we have lim k→0 + Ec m ν z, k 2 = ce m h, and With this in mind, we derive rigorous results for the Lamé functions and their eigenvalues, and deduce analogous results for Mathieu's equation using this limiting relation. For a general overview of the Lamé and Mathieu equations, see [2,5,29]. For a more detailed study of Mathieu's equation see also [17]. Whilst the results for asymptotic expansions of the Lamé functions and their respective eigenvalues for parameter ν large are not so abundant, analogous problems for the Mathieu functions and their eigenvalues for parameter h large have been studied extensively. The main results in the Lamé case can be found in [11,12,18,19,20]. These results are all formal, meaning they are not accompanied with any error analysis. None of the results about the functions have been published in [24, § 29.7], and only limited formal results about the corresponding eigenvalues can be found there. In [24, § 29.7(ii)] it is stated that one could derive from the results of [28] asymptotic approximations for the Lamé functions. However in that paper the results are given without much justification and the error bounds given for the approximations do not make sense in the intervals where the approximant is exponentially small. Here the results we give for the Lamé functions are new; we give a two term uniform asymptotic approximation for the Lamé functions in terms of parabolic cylinder functions for ν large complete with error bounds, and we show this holds uniformly in the interval z ∈ [0, K]. We also make a rigorous statement about the corresponding eigenvalues. We also give asymptotic expansions for the functions and eigenvalues in a shrinking neighbourhood of the origin, which correspond with the few formal results in the literature.
The main results in the Mathieu case can be found in [3,4,6,7,8,9,10,13,17,25,26,27,28]. In [17,28], error estimates are written down for a one term asymptotic approximation of the Mathieu functions but these are given without any justification, and the error does not make sense in the intervals where the approximant is exponentially small. In [13], similar results are given for the functions with similar issues for the error estimates; the results about the eigenvalues however seem reasonable, but methods of obtaining terms in the eigenvalue expansions seem cumbersome. The most satisfactory work thus far is contained in [4]. Here Dunster derives uniform asymptotic approximations for all complex values z when −2h 2 ≤ λ ≤ (2 − d)h 2 , d > 0, with error bounds either included or available for all approximations. These approximations involve both elementary functions and Whittaker functions. He also includes rigorous statements related to the eigenvalues a m and b m+1 . The remaining papers stated include only formal results, without any satisfactory error analysis. Here we consider only the real interval z ∈ (0, π) as many physical applications are restricted to real variables. The results we give are uniform asymptotic approximations complete with error bounds in the interval z ∈ 0, π 2 , given in the most natural form for the case we consider. Since we restrict ourselves to real variable analysis we can make stronger statements about the error bounds for the functions and their respective eigenvalues than given in [4]. We also give asymptotic expansions for the functions and eigenvalues in a shrinking neighbourhood of π 2 , which correspond with the formal results in the literature.

Overview
This paper is written in two parts: In Part I we derive two term uniform asymptotic approximations for the Lamé functions Ec m ν z, k 2 and Es m+1 ν z, k 2 which hold for z ∈ [0, K], and rigorous approximations for the large order behaviour of their respective eigenvalues a m ν and b m+1 ν , m = 0, 1, . . . , as κ → ∞, where κ = ν(ν + 1)k. Treating the Mathieu functions and their respective eigenvalues as a special case of those in the Lamé case, we obtain simply uniform asymptotic approximations for the Mathieu functions ce m (h, z) and se m+1 (h, z) which hold for z ∈ 0, π 2 , and rigorous approximations for the large order behaviour of their respective eigenvalues a m and b m+1 , as h → ∞.
Part II uses a simplified method to derive asymptotic expansions for both the Lamé and Mathieu functions and their respective eigenvalues. We can compute as many terms as we like in these expansions. The price to pay is that these function expansions only hold for z = O(κ −1/2 ) and z = π 2 + O(h −1/2 ), as κ → ∞ and h → ∞ respectively. These intervals at least encapsulate all of the interesting oscillatory behaviour of the functions. We give rigorous and realistic error bounds for the function expansions once truncated, along with order estimates for the error when the eigenvalue expansions are truncated.
First, we will summarise the relevant properties of the Lamé and Mathieu functions and their respective eigenvalues in the upcoming section. and interlace such that The normalisations in [24, § 28.2] can be rewritten as and to complete their definitions, the signs are determined by continuity from We summarise their properties and give boundary conditions in Table 2.
From formal asymptotic expansions given in [24, § 29.7] we deduce that s → 0 as κ → ∞, hence in this limit (4.2) has two coalescing turning points. The turning points of our equation are at x = ±s and thus we apply the theory of Case I in [21]. In this case uniform asymptotic approximations are in terms of the parabolic cylinder functions U − 1 2 κσ 2 , √ 2κζ and U − 1 2 κσ 2 , √ 2κζ . For the standard notation see [24, § 12.2].
Following Olver [21], new variables relating {x, w} to {ζ, W } are introduced by the appropriate Liouville transformation given by the dot signifying differentiation with respect to ζ, where σ is defined by From this we denote that 0 < s < 1 corresponds to 0 < σ < σ * , where σ * = 2 arcsin(k) πk .
Since ζ = ±σ corresponds to x = ±s, integration of the second of (4.4) yields −s These equations define ζ as a real analytic function of x. There is a one-to-one correspondence between the variables x and ζ, where ζ is an increasing function of x, and we denote ζ = −ζ * , 0, ζ * to correspond to x = −1, 0, 1. It follows that x(ζ, σ) is analytic both in ζ and σ for ζ ∈ (−ζ * , ζ * ) and σ ∈ (−σ * , σ * ). Alsoẋ is non-zero in these intervals. Performing the substitution t = sτ in (4.5) we expand the integral and obtain and by reversion In the critical case s = σ = 0 we have from the third of (4.6) which gives Thus we deduce that as s, σ → 0 ζ * → 2 arctanh(k) k .
The transformed differential equation is now of the form where ψ ζ, k 2 =ẋ 2 φ x, k 2 +ẋ 1/2 d 2 By construction, the apparent singularities in the above function at ζ = ±σ, corresponding to x = ±s, cancel each other out so that ψ ζ, k 2 is well-behaved there. To check this one could expand ψ ζ, k 2 around this point. One could also note that in (4.2), φ x, k 2 has singularities at x = ±1, where as ψ ζ, k 2 does not blow up there. However ζ and therefore ψ x, k 2 will have branch point singularities there. For our advantage in the error analysis, we write our differential equation in the form where ψ ζ, k 2 = ψ ζ, k 2 − ψ 0, k 2 = ψ ζ, k 2 + σ 4 − s 4 2σ 2 s 4 , and correspondingly where This gives ψ 0, k 2 = 0. (4.12) On inspection it follows that since s, σ and the variables x and ζ are all bounded as κ → ∞, and the apparent singularities at x = ±s and ζ = ±σ cancel each other, we have for ζ ∈ [−ζ * , ζ * ] ψ ζ, k 2 = O(1), uniformly in this limit. On applying Theorem I of [21, § 6] with u = κ, α = σ and ζ 2 = ζ * , we obtain the following solutions of (4.9) valid when ζ ∈ [0, ζ * ). Choosing as a consequence of (4.12) and thus defining the variational operator as the bounds obtained are .
The functions E, M and later N are defined in [21, § 5.8]. It follows from (4.12) and the evenness of ψ ζ, k 2 that thus clearly l 1 − 1 2 κ σ 2 is bounded as κ → ∞. Hence Applying similar analysis to the above, again from Theorem 1 in [21, § 6], we obtain We can extend this analysis to include the point ζ * since ζ and ψ ζ, k 2 are bounded there.
where C m ν , S m+1 ν , η m ν,c and η m+1 ν,s are constants, and s m ν , σ m ν , W m ν,1 and W m ν,2 are defined with respect to either h = a m ν or h = b m+1 ν in the previous section. Let's first consider m odd. In correspondence with the boundary conditions we require The requirement at z = K in (5.1) gives thus from (4.16), (4.15) and (4.13), and since as κ → ∞ (see [24, § 12.9]) necessarily η m ν,c is exponentially small in this limit. Note that here we have assumed that κ σ 2 = O(1) as κ → ∞, and we show later, in (5.3), that this assumption is justified.
Similarly the requirement at z = K in (5.2) gives is exponentially small as κ → ∞. Hence for both these cases by considering the requirements at z = 0 we have the condition and to satisfy this we require where j is an odd integer. Let's now consider m even. In correspondence with the boundary conditions we require By similar reason to the m odd case, when considering the boundary condition at z = K both η m ν,c and η m+1 ν,s are necessarily exponentially small as κ → ∞. Hence by using the boundary condition at z = 0 we have the condition To satisfy this we require that where j is an even integer.
√ 2κζ tend to the zeros of D j √ 2κζ , the parabolic cylinder functions with exactly j zeros in it's oscillatory interval (see [24, § 12.11]). Thus in correspondence with the number of zeros of the Lamé functions and those of D j √ 2κζ , we deduce that j = m and as such from (4.7) and (4.10) we have that We deduce from (4.7) and (4.3) that 6 Approximations in terms of parabolic cylinder D functions The Lamé functions decay exponentially on either side of the oscillatory interval in (−K, K). Our approximant in Section 4, U − 1 2 κ σ 2 , √ 2κζ , displays the appropriate exponentially decreasing behaviour when ζ is large and positive, but when the variable is large and negative it becomes exponentially large. If our argument − 1 2 κ σ 2 had been exactly a negative half-integer, which it is not, then the approximant would have exhibited this wanted exponentially decaying behaviour when the variable is both large and positive and large and negative.
With respect to (5.3) we define thus it makes sense to split up (4.8) so that From (6.1) and since s m ν , σ m ν and the variables x and ζ are all bounded as κ → ∞, we have for in this limit uniformly. Now our approximant will have the desired property of decaying exponentially on either side of the oscillatory interval for large positive and large negative ζ. We obtain the solutions for (6.2) We have that In Part II, Section 9, we will show that We combine this with (4.11) and obtain This time in the error analysis we choose it so that and, hence, take as the variational operator The bounds for the errors are .
We define for convenience From (4.14) clearly l m,1 is a bounded constant as κ → ∞. Hence we have Applying similar analysis to the above, again from Theorem 1 in [21, § 6], we obtain Thus from (6.3), (4.4) and (4.1) we obtain are O e −κζ 2 * κ m+1/2 as κ → ∞. Only when ζ nears the endpoint the contribution from the D m √ 2κζ + m ν,2 ζ, k 2 term is comparable to D m √ 2κζ + m ν,1 ζ, k 2 . These functions will be made unique by their normalisation.
A second term in the approximation in terms of parabolic cylinder D functions In [21, § 11] uniform asymptotic expansions are included. The two term version is of the form in which we have taken A 0 = 1 and for B 0 (ζ, κ) we have dt, whereσ = (2m + 1)/κ. This coefficient depends on κ and its dominant part is To simplify the integrand we take the final equation in (4.6) and let κ → ∞. The result is This relation can be inverted We use this in and are able to evaluate the integral in (6.6) and obtain 32ζB 0 (ζ) = k 2 + 1 ln 1 where C(ζ, k) = 2k coth kζ 2 /2 − k 2 − 1. Thus our two term approximation is as κ → ∞ and similarly for Es m+1 ν z, k 2 .

Normalisation constant
We now want an approximation for C m ν and S m+1 ν . Due to the complicated nature of the mapping between x and ζ, it is not simple to express one in terms of the other. With respect to (3.1) we consider the integral Using only this first term in an expansion for the solution, we can only obtain coefficients for up to 1/h terms. Any further terms in an expansion for the solution would contribute to further terms in the normalisation constant expansion, which we will not consider here. Since the first terms in their respective function uniform approximations are the same for large κ, it will be the same for both C m ν and S m+1 ν . First we perform the transformation x = sn(z, k) to obtain then another transformation to the ζ variable to obtain As the oscillatory behaviour of D m √ 2κζ happens in a shrinking region of the origin as κ → ∞, we seek to approximate the integral around this point to get an approximation for C m ν . Thus we consider an expansion of the form as ζ → 0, and substituting this into the relation From this we can determine, by matching (6.7) and (6.8), the c k terms, and obtain as ζ → 0. From this we obtain as ζ → 0. With respect to this, we consider as a first approximation for C m ν the integral and letting t = √ 2κζ, we have the approximation In correspondence with (3.2) and that for large z, both η c m and η s m+1 are O e −2hζ 2 * h m+1/2 as h → ∞, the relationship between z and x is defined by x = cos z and by evaluating the integrals (4.6), the relationship between x and ζ is defined by and λ m corresponds to either a m or b m+1 depending on the solution we are considering, and We have that Two term approximations are given by where for B 0 (ζ) we have the relation

Uniform asymptotic expansions
In Part I, uniform asymptotic approximations were given for the Lamé and Mathieu functions when a parameter became large. A third term in an asymptotic expansion could not be computed due to the complicated nature of the transformation of the independent variable. In Section 8 of this part we employ a simpler transformation than in Part I such that we can construct formal asymptotic expansions for the Lamé functions and their corresponding eigenvalues in the forms where t = √ 2κ sn (z, k), the coefficients A s (t) and B s (t) are just polynomials, and the µ s are constants in terms of m and k. The function expansions will clearly only make sense when t = O(1) as κ → ∞. In Section 9, we analyse the asymptotic expansions for the eigenvalues and give an order estimate for the error corresponding to the truncated eigenvalue expansion. Then in Section 10 we use these results to give rigorous and realistic error bounds for the function expansions upon truncation. In Sections 11 and 12 we identify our expansions with the Lamé functions and then summarise our results. Finally in Section 13 we give analogous results for the Mathieu functions and their corresponding eigenvalues as a special limiting case of those given for the Lamé functions and their corresponding eigenvalues.

Uniform asymptotic expansions of the Lamé functions
In this section we construct formal asymptotic expansions for solutions of Lamé's equation and corresponding eigenvalues. We do this by performing a simpler transformation than employed in Part I. The oscillatory behaviour of the Lamé functions happens in a shrinking neighbourhood of the origin as κ → ∞, and it can be shown that around the origin ζ behaves approximately like sn(z, k). Thus the variable in the parabolic cylinder function around this point behaves approximately like √ 2κ sn(z, k). This motivates the next simpler transformation.

t-plane
Letting t = √ 2κ sn(z, k) in (1.1) we have where z ∈ (−K, K) corresponds to t ∈ − √ 2κ, √ 2κ . Supposing in accordance with (5.4) that where n is a positive integer and µ n can be re-expanded in a sensible manner. We can then write Lamé's equation in the form This equation is split in such a way that constructing a formal asymptotic expansion in terms of parabolic cylinder functions D m (t) in the form seems sensible. However one should observe that this splitting only makes sense when t = O(1) as h → ∞, hence this formal expansion is only sensible for this range of t. We denote this range as (−t * , t * ). One should note that whilst this appears to be a new ansatz, there are similar expansions given in the literature for the Mathieu functions. These expansions are given in terms of D m−4j (t) for j ∈ Z instead of in terms of D m (t) and D m (t), but using the recurrence relations for the parabolic cylinder functions it is obvious these expansions are equivalent up to normalisation. An expansion in the form we have given allows us to differentiate easily and hence seems the most natural in this case, also allowing us to perform rigorous error analysis. We seek solutions which are either even or odd respective to the parity of m. Since D m (t) is either even or odd respective to when m is either even or odd, we deduce that A s (t) and B s (t) must be even and odd respectively. Substituting (8.3) into (8.2) and equating powers of κ, we have the recurrence relations for A s (t) and B s (t) +t(m) tA s−2 (t) + 2tB s−2 (t) + 2B s−2 (t) = 0, (8.5) wheret(m) = 1 4 t 2 − m − 1 2 . Neither of these relations separately determine solutions for A s (t) or B s (t) from previous coefficients, thus we differentiate (8.4) to obtain an expression for A s (t) and substitute it into (8.5); this gives the third order inhomogeneous differential equation for B s (t) Once B s (t) is determined, we can use (8.4) to determine A s (t). There will be freedom in choosing the integration constants in the A s (t) terms, with identification of our solutions made unique by their normalisation.

General coef f icients B s (t) and A s (t)
Using variation of parameters we obtain the general solution for B s (t) By expanding the Wronskians we derive the relations Then without loss of generality we can rewrite the indefinite integrals in b i s (t) 3 i=1 as definite integrals from 0 to t, and since B s (t) is supposed to be an odd function we have to take c 1 = c 2 = 0. Thus Although we consider t to be O(1) as κ → ∞, we still want an expansion which exhibits the correct behaviour at the endpoints of the interval. Expanding the squared term in (8.8) and splitting it into three separate integrals, large variable asymptotics for the parabolic cylinder functions (see [24, § 12.9]) tells us that the terms involving D 2 m (t) or D m (t)D m (t) grow no faster than polynomials when t becomes large. Since to ensure the boundedness of our formal expansion as t becomes large, µ s is determined uniquely by the condition that This condition will also ensure boundedness as t → −∞. Note that in the next section we will derive an alternative method to compute the µ s terms which we use to compute the coefficients as it is simpler. Consider first s = 0. From (8.7) we see that b 0 (t) = 0, thus ensuring that B 0 (t) is odd we obtain the general solution and then Rearranging, on the κ 0 level of the asymptotic expansion (8.3) we have Thus having this term in B 0 (t) equates to B 0 (t) = 0 and an extra constant term in A 0 (t), and since we have freedom in the arbitrary constant terms in A s (t), we can take c 3 = 0 without loss of generality. For simplicity we take A 0 (t) = 1 and adopt the convention A s (t) = 0 for s ≥ 1.
With these choices we have and expressing the D m (t) in (8.9) via [24, formula (12.7.2)] in terms of Hermite polynomials we can evaluate the integral in (8.9) and obtain In the case s = 1, it can be shown using integration by parts that for µ 1 which satisfies (8.9) Thus clearly in this case, if the correct c 3 is chosen in (8.8) then B 1 (t) is exactly an odd polynomial. Using similar observations as in the s = 0 case, we also note that if this multiple of D m D m is included in the B 1 (t) term, the expansion can be rearranged so that this term is instead represented in the constant term of A 1 (t). Taking B 1 (t) to be exactly an odd polynomial, we get that A 1 (t) is an even polynomial. If one would go through the details to compute the representation for B 1 (t) as a polynomial plus this multiple of D m D m , one would see that it would appear for all s ≥ 1 that if the previous A s (t) and B s (t) terms are all polynomials, then the A s (t) and B s (t) terms can be represented as polynomials.

Polynomial coef f icients B s (t) and A s (t)
To obtain explicit expressions for A s (t) and B s (t) we try substituting in polynomial expansions with undetermined coefficients. Take A 0 (t) = 1 and B 0 (t) = 0, then for s ≥ 1 we consider A s (t) and B s (t) in the form Substituting these into (8.6) and (8.4) we obtain the recurrence relations for coefficients (2i + 3)(2i + 2)(2i + 1)b s,i+1 + (4m + 2)(2i + 1)b s,i − 2ib s,i−1 From these recurrence relations it is observed that only a finite number of the a s,i and b s,i are non-zero. The orders are displayed in Table 3.
In this manner the coefficients are determined recursively. We deduce by considering the difference of (8.10) and (8.11) in the case i = 0 that µ s = (2m + 1)b s,0 − 2a s,1 . (8.12) In the case that s is even then we start with i = 2s − 1 in (8.10) to determine b s,2s−2 . Then we take i = 2s − 2 and determine b s,2s−3 , and so on. Once every power of t is eliminated, we are left with a constant equation which we must make zero by specifying µ s ; in this manner the eigenvalue terms are determined uniquely. These µ s terms are the same as those specified by the condition (8.9), since it is required for solutions B s (t) which grow no faster than polynomials. Once the coefficients in the polynomial expression for B s (t) are determined, and the corresponding µ s , then the coefficients in the polynomial expression for A s (t) can be determined from (8.11). In the case that s is odd then we have to start with i = 2s in (8.10) to determine b s,2s−1 .
Returning to the z-plane: odd solutions In [24, § 28.8] expansions for the odd Mathieu functions are given which exhibit the correct odd behaviour, which an expansion in the form (8.3) would not when transformed back into the z-plane. We want something similar for the Lamé functions.
Once we transform our formal expansion back into the z-plane, over the whole real line they behave like the even Lamé functions Ec m ν z, k 2 . We want to construct a similar expansions which when transformed back into the z-plane behave like the odd Lamé functions Es m+1 ν z, k 2 . We need our expansion to have the property The Jacobi elliptic function cn(z, k) is even around the origin and odd around z = −K and z = K. Hence again letting t = √ 2κ sn(z, k) we consider a formal expansion for a solution of (8.2) of the form This has the correct behaviour at z = −K and z = K. Again we will require that P s (t) is even and Q s (t) is odd. By writing we can express the formal solution in the form Expanding this square root, we can rewrite this formal expansion as where the connection between P s (t) and A s (t), and Q s (t) and B s (t) is given by where a b is the generalised binomial coefficient. Thus we determine the P s (t) and Q s (t) terms by connection with the A s (t) and B s (t) derived previously. By connection with the previous case we get the same eigenvalue expansion to all orders for both a m ν and b m+1 ν .

An interlude: error analysis for the eigenvalue expansions
In this section we match our eigenvalue expansions with the eigenvalues a m ν and b m+1 ν , and give order estimates for the expansions on truncation. The theory of this subsection is included in [17], where it is employed for Mathieu's equation. For this we consider Sturm-Liouville theory; for a fuller treatment on Sturm-Liouville theory see [1]. We have the differential equation The regular SL theory states that these eigenvalues are real and can be ordered such that λ 0 < λ 1 < · · · < λ ν < · · · → ∞ and these eigenvalues λ ν correspond to unique eigenfunctions y ν (x) with exactly ν zeros in (a, b). Normalised eigenfunctions form an orthonormal basis such that The following theorem is Theorem 1 in [17, § 1.52] Theorem 9.1. Consider SL to be a self-adjoint differential operator defined on a subspace A of L 2 [[a, b], w(x)dx] containing all the twice differentiable functions which satisfy the boundary conditions (9.1) and corresponding λ ν such that Let y(x) ∈ A be such that ( y, y) = 1. Now define a remainder function R(x) such that Proof . We have that (R, y k ) = ((SL + λ) y, y k ) = ( y, SLy k ) + λ( y, y k ) = (λ − λ k )( y, y k ), and using this in and these eigenvalues are real and can be ordered such that

Eigenvalues of Lamé's equation
The functions Ec m ν z, k 2 have exactly m zeros in (−K, K). The normalised eigenfunctions w m ν z, k 2 form an orthonormal basis such that Let L ν be the operator Letting t = √ 2κ sn(z, k), we define the truncated expansions corresponding to w m ν z, k 2 as where the A s (t) and B s (t) terms were derived previously, and c m ν,n is defined to be a function of κ so that where the dash represents differentiation with respect to t. Thus clearly w m ν,n z, k 2 ∈ A since cn(−K, k) = cn(K, k) = 0. We also define the truncated eigenvalue expansions a m ν,n = (2m + 1)κ + 2 where the µ s were derived previously and define R m ν,n z, k 2 such that L ν + a m ν,n w m ν,n z, k 2 = R m ν,n z, k 2 .
We consider the operator L ν acting on D m (t) and derive L ν (D m (t)) = (2m + 1) 1 + k 2 2 t 2 − κ − k 2 (2m + 1) + κ 1 + k 2 t 4 4κ + k 2 t 6 8κ D m (t) Using the recurrence relations tD m (t) = D m+1 (t) + mD m−1 (t), D m (t) = mD m−1 (t) − 1 2 tD m (t), (9.5) we can rewrite w m ν,n z, k 2 given in (9.2) such that for s ∈ {0, . . . , n}, A s (t)D m (t) and B s (t)D m (t) are sums of varying orders of parabolic cylinder functions with constant coefficients dependent only on m and k; it then follows from (9.4) that with the solution rewritten in this form, we have R m ν,n z, k 2 = c m ν,n κ −n (· · · ) + κ −1 (· · · ) + · · · , where the terms inside the brackets are sums of varying orders of parabolic cylinder functions with constant coefficients dependent only on k and m. Note that this remainder term will be finite sum and clearly R m ν,n z, k 2 ∈ L 2 [[−K, K], dz]. As w m ν,n z, k 2 ∈ A we have that min k a m ν,n − a k ν 2 < R m ν,n , R m ν,n .
In Part I, Section 5, we proved that We need an order estimate for R m ν,n , R m ν,n . The integrals we must consider then are of the form where i, j ∈ N 0 . Performing the substitution t = √ 2κ sn(z, k) we obtain Since when t is large D m is exponentially small it follows that Considering the first term of the expansion for w m ν,n z, k 2 in (9.2) we have as κ → ∞. Hence it follows from (9.3) that c m,n = O(κ 1/4 ) as κ → ∞. Similar observations would give us that and so as κ → ∞. The error analysis for b m ν,n would in the same manner give the order estimate Hence we have from (8.1) that a m ν = (2m + 1)κ + 2 In the same manner we would deduce that b m+1 ν = (2m + 1)κ + 2 This result also follows from the difference between a m ν and b m+1 ν being exponentially small as κ → ∞ (see [24, § 28.8]).

Error analysis for the uniform asymptotic expansions of the Lamé functions
We now use the results of Section 9 to obtain strict and realistic error bounds for the functions expansions derived in Section 8. Define the differential operator and consider t ∈ (−t * , t * ). We have the truncated expansion corresponding to an even solution such that we have the exact solution w m ν t, k 2 = w m ν,n t, k 2 + m ν,n t, k 2 . (10.1) We define the remainder term R m ν,n t, k 2 such that L ν w m ν,n t, k 2 = R m ν,n t, k 2 and split the eigenvalue such that h 2κ = m + 1 2 + n s=1 µ s κ s + µ n+1 κ n+1 , and note we just proved that µ n+1 = O(1) as κ → ∞. Since the coefficients A s (t) and B s (t) satisfy (8.4) and (8.5) it follows that R m ν,n t, k 2 = O κ −n−1 , as κ → ∞. Applying L ν to (10.1) we obtain m ν,n and denoting the right hand side of this equation Ω m ν,n t, k 2 , by use of variation of parameters we have In accordance with [22, § 6.10.2] we define J(τ ) = 1, H(τ ) = 1 + k 2 − k 2 τ 2 2κ and Since we consider t ∈ (−t * , t * ) where t * = O(1) as κ → ∞, the error analysis is much simpler than the analysis Olver uses in [21] as we have |K (t, τ ) | ≤ k 0 , and |∂K (t, τ ) /∂t| ≤ k 1 , where k 0 and k 1 are O(1) as κ → ∞. Thus again in accordance with [22, § 6.10.2] we define P 0 (t) = k 0 , Q(τ ) = 1, P 1 (t) = k 1 (10.2) (we do not define P 2 (t) as we do not need to bound |∂ 2 K(t, τ )/∂t 2 | to carry out our analysis), and finally the constants κ = 1, κ 0 = k 0 , κ 1 = k 1 .
These eigenvalue coefficients match the formal results given in [24, § 29.7].