Asymptotics for the Lamé and Mathieu Functions

Rigorous Asymptotics for the Lamé and Mathieu
Functions and their Respective Eigenvalues
with a Large Parameter\star\star$\star$This paper is a contribution to the
Special Issue on Orthogonal Polynomials, Special Functions and Applications.
The full collection is available at
http://www.emis.de/journals/SIGMA/OPSFA2015.html

Karen OGILVIE and Adri B. OLDE DAALHUIS K. Ogilvie and A.B. Olde Daalhuis
Maxwell Institute and School of Mathematics, The University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh EH9 3FD, UK<EMAIL_ADDRESS><EMAIL_ADDRESS>
Received August 03, 2015, in final form November 20, 2015; Published online
November 24, 2015

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é 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.

Lamé functions; Mathieu functions; uniform asymptotic approximations;
coalescing turning points

33E10; 34E05; 34E20

###  1 Introduction

The Lamé equation is

$\displaystyle\frac{d^{2}w}{dz^{2}}+\left(h-\nu(\nu+1)k^{2}\operatorname{sn}^{2}(z,k)\right)w=0,$
(1.1)

where $h$, $k$ and $\nu$ are real parameters such that $0<k<1$ and
$\nu\geq-\frac{1}{2}$, and $\operatorname{sn}(z,k)$ is a Jacobian elliptic
function (see [24, § 22.2]). We consider the interval $z\in(-K,K)$, where
$K=K(k)$ is Legendre’s complete elliptic integral of the first kind (see [24,
§ 19.2(ii)]). When $h$ assumes the special values $a_{\nu}^{m}$ or
$b_{\nu}^{m+1}$ for $m=0,1,\dots$, Lamé’s equation admits even or odd periodic
solutions denoted $\operatorname{Ec}_{\nu}^{m}(z,k^{2})$ or
$\operatorname{Es}_{\nu}^{m+1}(z,k^{2})$ respectively.

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

$\displaystyle\frac{d^{2}w}{dz^{2}}+\left(\lambda-2h^{2}\cos 2z\right)w=0,$
(1.2)

where $\lambda$ and $h$ are real parameters. We consider the interval
$z\in(0,\pi)$. When $\lambda$ assumes the special values $a_{m}$ or $b_{m+1}$
for $m=0,1,2,\dots$, Mathieu’s equation admits even or odd periodic solutions
denoted $\operatorname{ce}_{m}(h,z)$ or $\operatorname{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 $\nu$ 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\to 0_{+}$ from [24, § 22.5(ii)] that

$\displaystyle\lim_{k\to
0_{+}}\operatorname{sn}(z,k)=\sin{z},\qquad\textrm{and}\qquad\lim_{k\to
0_{+}}K(k)=\frac{\pi}{2},$

thus if $\nu\to\infty$ in such a way that as $k\to 0_{+}$,
$\sqrt{\nu(\nu+1)}k=2h$ for some constant $h$, then we can rewrite the limit
of Lamé’s equation in the form

$\displaystyle\frac{d^{2}w}{dz^{2}}+\left(h_{L}-2h^{2}-2h^{2}\cos{2\left(\frac{\pi}{2}-z\right)}\right)w=0.$

Thus for $\nu=-\frac{1}{2}+\sqrt{\frac{1}{4}+\left(\frac{2h}{k}\right)^{2}}$
we have

$\displaystyle\lim_{k\to
0_{+}}\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=\operatorname{ce}_{m}\left(h,\frac{\pi}{2}-z\right),\qquad\lim_{k\to
0_{+}}\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)=\operatorname{se}_{m+1}\left(h,\frac{\pi}{2}-z\right),$
(1.3)

and

$\displaystyle a_{m}=\lim_{k\to 0_{+}}a_{\nu}^{m}-2h^{2},\qquad
b_{m+1}=\lim_{k\to 0_{+}}b_{\nu}^{m+1}-2h^{2}.$ (1.4)

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 $\nu$ 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 $\nu$ large complete with error bounds, and we show this holds
uniformly in the interval $z\in[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}\leq\lambda\leq(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\in(0,\pi)$ 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\in\big{[}0,\frac{\pi}{2}\big{]}$, 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 $\frac{\pi}{2}$, which correspond with the formal results in
the literature.

###  2 Overview

This paper is written in two parts: In Part I we derive two term uniform
asymptotic approximations for the Lamé functions
$\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)$ and
$\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)$ which hold for
$z\in[0,K]$, and rigorous approximations for the large order behaviour of
their respective eigenvalues $a_{\nu}^{m}$ and $b_{\nu}^{m+1}$, $m=0,1,\dots$,
as $\kappa\to\infty$, where $\kappa=\sqrt{\nu(\nu+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 $\operatorname{ce}_{m}(h,z)$ and $\operatorname{se}_{m+1}(h,z)$
which hold for $z\in\big{[}0,\frac{\pi}{2}\big{]}$, and rigorous
approximations for the large order behaviour of their respective eigenvalues
$a_{m}$ and $b_{m+1}$, as $h\to\infty$.

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=\mathcal{O}(\kappa^{-1/2})$ and
$z=\frac{\pi}{2}+\mathcal{O}(h^{-1/2})$, as $\kappa\to\infty$ and $h\to\infty$
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.

###  3 Properties of Lamé and Mathieu functions

#### Lamé functions

We will summarise their important properties here, for a full treatment see
[24, § 29]. These functions are either $2K$-periodic or $2K$-antiperiodic,
depending on the parity of $m$. The functions have exactly $m$ zeros in the
interval $(-K,K)$ and their eigenvalues are ordered such that

$\displaystyle a_{\nu}^{0}<a_{\nu}^{1}<a_{\nu}^{2}<\cdots,\qquad
b_{\nu}^{1}<b_{\nu}^{2}<b_{\nu}^{3}<\cdots,$

and interlace such that

$\displaystyle a_{\nu}^{m}<b_{\nu}^{m+1},\qquad b_{\nu}^{m}<a_{\nu}^{m+1}.$

The eigenvalues coalesce such that

$\displaystyle
a_{\nu}^{m}=b_{\nu}^{m},\qquad\textrm{when}\quad\nu=0,1,\dots,m-1.$

Since the Jacobian elliptic function $\operatorname{dn}(z,k)$ (see [24, § 22])
is even, we can rewrite the normalisations given in [24, § 29.3] as

$\displaystyle\int_{-K}^{K}\operatorname{dn}(z,k)\left\\{\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)\right\\}^{2}dz=\int_{-K}^{K}\operatorname{dn}(z,k)\left\\{\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)\right\\}^{2}dz=\frac{\pi}{2},$
(3.1)

To complete their definitions we have

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(K,k^{2}\right)>0,\qquad{\rm
and}\qquad\frac{d\operatorname{Es}_{\nu}^{m}\left(z,k^{2}\right)}{dz}\bigg{|}_{z=K}<0.$
(3.2)

They satisfy the orthogonality conditions for $m\neq n$, $(n=0,1,\dots)$

$\displaystyle\int_{-K}^{K}\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)\operatorname{Ec}_{\nu}^{n}\left(z,k^{2}\right)dz=0,\qquad{\rm
and}\qquad\int_{-K}^{K}\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)\operatorname{Es}_{\nu}^{n+1}\left(z,k^{2}\right)dz=0.$

We summarise their properties and give boundary conditions in Table 1.

Table 1: Properties and boundary conditions for Lamé functions. eigenfunctions | eigenvalues | periodicity | $\begin{array}[]{@{}c@{}}\text{parity}\\\ \text{at $z=K$}\end{array}$ | $\begin{array}[]{@{}c@{}}\text{parity}\\\ \text{at $z=0$}\end{array}$ | boundary conditions
---|---|---|---|---|---
$\operatorname{Ec}_{\nu}^{2m}\left(z,k^{2}\right)$ | $a_{\nu}^{2m}$ | period $2K$ | even | even |  $w^{\prime}(0)=w^{\prime}(K)=0$ 2pt1pt
$\operatorname{Ec}_{\nu}^{2m+1}\left(z,k^{2}\right)$ | $a_{\nu}^{2m+1}$ | antiperiod $2K$ | even | odd |  $w(0)=w^{\prime}(K)=0$ 1pt
$\operatorname{Es}_{\nu}^{2m+1}\left(z,k^{2}\right)$ | $b_{\nu}^{2m+1}$ | antiperiod $2K$ | odd | even |  $w^{\prime}(0)=w(K)=0$ 1pt
$\operatorname{Es}_{\nu}^{2m+2}\left(z,k^{2}\right)$ | $b_{\nu}^{2m+2}$ | period $2K$ | odd | odd |  $w(0)=w(K)=0$1pt

#### Mathieu functions

We will summarise their important properties here, for a full treatment see
[24, § 28]. These are either $\pi$-periodic or $\pi$-antiperiodic, depending
on the parity of $m$. Both functions have $m$ zeros in the interval $(0,\pi)$
and their eigenvalues are ordered such that

$\displaystyle a_{0}<a_{1}<\cdots\to\infty,\qquad
b_{1}<b_{2}<\cdots\to\infty,$

and interlace such that

$\displaystyle a_{0}<b_{1}<a_{1}<b_{2}<a_{2}<\cdots.$

The normalisations in [24, § 28.2] can be rewritten as

$\displaystyle\int_{0}^{\pi}\left\\{\operatorname{ce}_{m}(h,z)\right\\}^{2}dz=\int_{0}^{\pi}\left\\{\operatorname{se}_{m+1}(h,z)\right\\}^{2}dz=\frac{\pi}{2},$

and to complete their definitions, the signs are determined by continuity from

$\displaystyle\operatorname{ce}_{0}(0,z)=\frac{1}{\sqrt{2}},\qquad\operatorname{ce}_{m}(0,z)=\cos{mz},\qquad\operatorname{se}_{m}(0,z)=\sin{mz}.$

We summarise their properties and give boundary conditions in Table 2.

Table 2: Properties and boundary conditions for Mathieu functions. eigenfunctions | eigenvalues | periodicity | $\begin{array}[]{@{}c@{}}\text{parity of}\\\ \text{$w(h,z)$}\end{array}$ | $\begin{array}[]{@{}c@{}}\text{parity of}\\\ \text{$w\left(h,z+\frac{\pi}{2}\right)$}\end{array}$ | boundary conditions
---|---|---|---|---|---
$\operatorname{ce}_{2m}(h,z)$ | $a_{2m}$ | period $\pi$ | even | even |  $w^{\prime}(0)=w^{\prime}\left(\pi/2\right)=0$ 2pt1pt
$\operatorname{ce}_{2m+1}(h,z)$ | $a_{2m+1}$ | antiperiod $\pi$ | even | odd |  $w^{\prime}(0)=w\left(\pi/2\right)=0$ 1pt
$\operatorname{se}_{2m+1}(h,z)$ | $b_{2m+1}$ | antiperiod $\pi$ | odd | even |  $w(0)=w^{\prime}\left(\pi/2\right)=0$ 1pt
$\operatorname{se}_{2m+2}(h,z)$ | $b_{2m+2}$ | period $\pi$ | odd | odd |  $w(0)=w\left(\pi/2\right)=0$1pt

##  Part I Uniform asymptotic approximations

We wish to obtain uniform approximations for the Lamé functions
$\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)$ and
$\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)$, and rigorous
approximations for their respective eigenvalues $a_{\nu}^{m}$ and
$b_{\nu}^{m+1}$, for $m=0,1,\dots$. In [24, § 29.7], the first few terms are
given for formal asymptotic expansions of the eigenvalues $a_{\nu}^{m}$ and
$b_{\nu}^{m+1}$ as $\nu\to\infty$. These indicate that (1.1) will have two
turning points, one either side of the origin, which coalesce symmetrically
into the origin of the $z$-plane as $\nu\to\infty$. In Section 4 we will use
this fact and derive uniform approximations for solutions in terms of the
parabolic cylinder functions
$U\big{(}{-}\frac{1}{2}\kappa\sigma^{2},\sqrt{2\kappa}\zeta\big{)}$ and
$\overline{U}\big{(}{-}\frac{1}{2}\kappa\sigma^{2},\sqrt{2\kappa}\zeta\big{)}$
as $\nu\to\infty$, where $\kappa=\sqrt{\nu(\nu+1)}k$ and $\zeta$ arises from a
special complicated transformation of the variable $z$, and $\sigma$ is
related to the eigenvalue parameter $h$. Only when
$-\frac{1}{2}\kappa\sigma^{2}$ is exactly a negative half integer does this
function decay exponentially on both the positive and negative real axis when
the variable is large (see [24, § 12.9]). Respective of this, in Section 5 we
get a rigorous approximation for the eigenvalues $a_{\nu}^{m}$ and
$b_{\nu}^{m+1}$ and use this in Section 6 to give an approximation in terms of
the special parabolic cylinder functions
$D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}$, where
$D_{m}(z)=U\left(-m-\frac{1}{2},z\right)$, and these decay exponentially for
both large positive and large negative values of the variable, and have
exactly $m$ zeros in their oscillatory interval (again see [24, § 12.2]).
Finally in Section 7 we use these results for the Lamé equation to obtain in a
limiting form uniform approximations for the Mathieu functions
$\operatorname{ce}_{m}(h,z)$ and $\operatorname{se}_{m+1}(h,z)$, and rigorous
results for their eigenvalues $a_{m}$ and $b_{m+1}$, as $h\to\infty$.

###  4 Uniform approximations for the Lamé functions

The periodic coefficient in (1.1) is troublesome, thus we transform the
independent variable to obtain an algebraic equation, and then transform the
dependent variable to remove the subsequent term in the first derivative; this
is done by letting

$\displaystyle x=\operatorname{sn}(z,k),\qquad
w\left(z,k^{2}\right)=\frac{1}{\left(\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)\right)^{1/4}}\widetilde{w}\left(x,k^{2}\right),$
(4.1)

and denoting $\kappa=\sqrt{\nu(\nu+1)}k$, we obtain Lamé’s equation in the
form

$\displaystyle\frac{d^{2}\widetilde{w}\left(x,k^{2}\right)}{dx^{2}}=\left(\kappa^{2}\frac{x^{2}-s^{2}}{\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}+\phi\left(x,k^{2}\right)\right)\widetilde{w}\left(x,k^{2}\right),$
(4.2)

where

$\displaystyle s^{2}=\frac{h}{\kappa^{2}}$ (4.3)

and

$\displaystyle\phi\left(x,k^{2}\right)=-\frac{2k^{2}\left(k^{2}+1\right)x^{4}+\left(k^{4}-10k^{2}+1\right)x^{2}+2\left(1+k^{2}\right)}{4\left(1-x^{2}\right)^{2}\left(1-k^{2}x^{2}\right)^{2}}.$

We correspondingly consider the interval $x\in(-1,1)$, where $x=-1,0,1$
corresponds to $z=-K,0,K$. We now write $\kappa\to\infty$ to correspond to
$\nu\to\infty$.

From formal asymptotic expansions given in [24, § 29.7] we deduce that $s\to
0$ as $\kappa\to\infty$, hence in this limit (4.2) has two coalescing turning
points. The turning points of our equation are at $x=\pm s$ and

$\displaystyle\frac{x^{2}-s^{2}}{\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}<0,\qquad-s<x<s,$

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\big{(}{-}\frac{1}{2}\kappa\sigma^{2},\sqrt{2\kappa}\zeta\big{)}$ and
$\overline{U}\big{(}{-}\frac{1}{2}\kappa\sigma^{2},\sqrt{2\kappa}\zeta\big{)}$.
For the standard notation see [24, § 12.2].

Following Olver [21], new variables relating
$\left\\{x,\widetilde{w}\right\\}$ to $\left\\{\zeta,W\right\\}$ are
introduced by the appropriate Liouville transformation given by

$\displaystyle
W\left(\zeta,k^{2}\right)=\dot{x}^{-\tfrac{1}{2}}\widetilde{w}\left(x,k^{2}\right),\qquad\dot{x}^{2}\frac{x^{2}-s^{2}}{\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}=\zeta^{2}-\sigma^{2}$
(4.4)

the dot signifying differentiation with respect to $\zeta$, where $\sigma$ is
defined by

$\displaystyle\int_{-s}^{s}\sqrt{\frac{s^{2}-t^{2}}{\left(1-t^{2}\right)\left(1-k^{2}t^{2}\right)}}dt=\int_{-\sigma}^{\sigma}\sqrt{\sigma^{2}-\tau^{2}}d\tau=\frac{\pi}{2}\sigma^{2}.$
(4.5)

From this we denote that

$\displaystyle 0<s<1\quad\textrm{corresponds to}\quad
0<\sigma<\sigma_{*},\qquad\textrm{where}\quad\sigma_{*}=2\sqrt{\frac{\arcsin(k)}{\pi
k}}.$

Since $\zeta=\pm\sigma$ corresponds to $x=\pm s$, integration of the second of
(4.4) yields

$\displaystyle\int_{x}^{-s}\sqrt{\frac{t^{2}-s^{2}}{\left(1-t^{2}\right)\left(1-k^{2}t^{2}\right)}}dt=\int_{\zeta}^{-\sigma}\sqrt{\tau^{2}-\sigma^{2}}d\tau,\qquad$
$\displaystyle-1<x\leq-s,$
$\displaystyle\int_{-s}^{x}\sqrt{\frac{s^{2}-t^{2}}{\left(1-t^{2}\right)\left(1-k^{2}t^{2}\right)}}dt=\int_{-\sigma}^{\zeta}\sqrt{\sigma^{2}-\tau^{2}}d\tau,\qquad$
$\displaystyle-s\leq x\leq s,$
$\displaystyle\int_{s}^{x}\sqrt{\frac{t^{2}-s^{2}}{\left(1-t^{2}\right)\left(1-k^{2}t^{2}\right)}}dt=\int_{\sigma}^{\zeta}\sqrt{\tau^{2}-\sigma^{2}}d\tau,\qquad$
$\displaystyle s\leq x<1.$ (4.6)

These equations define $\zeta$ as a real analytic function of $x$. There is a
one-to-one correspondence between the variables $x$ and $\zeta$, where $\zeta$
is an increasing function of $x$, and we denote $\zeta=-\zeta_{*},0$,
$\zeta_{*}$ to correspond to $x=-1,0,1$. It follows that $x(\zeta,\sigma)$ is
analytic both in $\zeta$ and $\sigma$ for $\zeta\in(-\zeta_{*},\zeta_{*})$ and
$\sigma\in(-\sigma_{*},\sigma_{*})$. Also $\dot{x}$ is non-zero in these
intervals.

Performing the substitution $t=s\tau$ in (4.5) we expand the integral and
obtain

$\displaystyle\sigma^{2}=s^{2}+\frac{1+k^{2}}{8}s^{4}+\frac{3+2k^{2}+3k^{4}}{64}s^{6}+\mathcal{O}\left(s^{8}\right),\qquad
s\to 0,$

and by reversion

$\displaystyle
s^{2}=\sigma^{2}-\frac{1+k^{2}}{8}\sigma^{4}-\frac{\left(1-k^{2}\right)^{2}}{64}\sigma^{6}+\mathcal{O}\left(\sigma^{8}\right),\qquad\sigma\to
0.$ (4.7)

In the critical case $s=\sigma=0$ we have from the third of (4.6)

$\displaystyle\frac{\operatorname{arctanh}(k)}{k}=\int_{0}^{1}\frac{t}{\sqrt{(1-t^{2})(1-k^{2}t^{2})}}dt=\int_{0}^{\zeta_{*}}\tau
d\tau=\textstyle{{1}\over{2}}\displaystyle\zeta_{*}^{2},$

which gives

$\displaystyle\zeta_{*}=\sqrt{\frac{2\operatorname{arctanh}(k)}{k}}.$

Thus we deduce that as $s,\sigma\to 0$

$\displaystyle\zeta_{*}\to\sqrt{\frac{2\operatorname{arctanh}(k)}{k}}.$

The transformed differential equation is now of the form

$\displaystyle\frac{d^{2}W\left(\zeta,k^{2}\right)}{d\zeta^{2}}=\left(\kappa^{2}\left(\zeta^{2}-\sigma^{2}\right)+{\psi}\left(\zeta,k^{2}\right)\right)W\left(\zeta,k^{2}\right),$
(4.8)

where

$\displaystyle{\psi}\left(\zeta,k^{2}\right)=\dot{x}^{2}\phi\left(x,k^{2}\right)+\dot{x}^{1/2}\frac{d^{2}}{d\zeta^{2}}\big{(}\dot{x}^{-1/2}\big{)}$
$\displaystyle\hphantom{{\psi}\left(\zeta,k^{2}\right)}{}=\frac{2\sigma^{2}+3\zeta^{2}}{4\left(\sigma^{2}-\zeta^{2}\right)^{2}}+\frac{\sigma^{2}-\zeta^{2}}{4}\left(-k^{2}+\frac{1+k^{2}\left(1-3s^{2}\right)}{s^{2}-x^{2}}\right.$
$\displaystyle\left.\hphantom{{\psi}\left(\zeta,k^{2}\right)=}{}+\frac{3\left(1-2\left(1+k^{2}\right)s^{2}+3k^{2}s^{4}\right)}{\left(s^{2}-x^{2}\right)^{2}}+\frac{5s^{2}\left(s^{2}-1\right)\left(1-k^{2}s^{2}\right)}{\left(s^{2}-x^{2}\right)^{3}}\right).$

By construction, the apparent singularities in the above function at
$\zeta=\pm\sigma$, corresponding to $x=\pm s$, cancel each other out so that
${\psi}\left(\zeta,k^{2}\right)$ is well-behaved there. To check this one
could expand ${\psi}\left(\zeta,k^{2}\right)$ around this point. One could
also note that in (4.2), $\phi\left(x,k^{2}\right)$ has singularities at
$x=\pm 1$, where as ${\psi}\left(\zeta,k^{2}\right)$ does not blow up there.
However $\zeta$ and therefore $\psi\left(x,k^{2}\right)$ will have branch
point singularities there.

For our advantage in the error analysis, we write our differential equation in
the form

$\displaystyle\frac{d^{2}W\left(\zeta,k^{2}\right)}{d\zeta^{2}}=\big{(}\kappa^{2}\left(\zeta^{2}-\widetilde{\sigma}^{2}\right)+{\widetilde{\psi}}\left(\zeta,k^{2}\right)\big{)}W\left(\zeta,k^{2}\right),$
(4.9)

where

$\displaystyle\widetilde{\psi}\left(\zeta,k^{2}\right)={\psi}\left(\zeta,k^{2}\right)-{\psi}\left(0,k^{2}\right)={\psi}\left(\zeta,k^{2}\right)+\frac{\sigma^{4}-s^{4}}{2\sigma^{2}s^{4}},$

and correspondingly

$\displaystyle\widetilde{\sigma}^{2}={\sigma}^{2}+\frac{\sigma^{4}-s^{4}}{2\kappa^{2}\sigma^{2}s^{4}},$
(4.10)

where

$\displaystyle\frac{\sigma^{4}-s^{4}}{2\sigma^{2}s^{4}}=\frac{1+k^{2}}{8}+\mathcal{O}\left(s^{2}\right),\qquad
s\to 0.$ (4.11)

This gives

$\displaystyle\widetilde{\psi}\left(0,k^{2}\right)=0.$ (4.12)

On inspection it follows that since $s$, $\sigma$ and the variables $x$ and
$\zeta$ are all bounded as $\kappa\to\infty$, and the apparent singularities
at $x=\pm s$ and $\zeta=\pm\sigma$ cancel each other, we have for
$\zeta\in[-\zeta_{*},\zeta_{*}]$

$\displaystyle\widetilde{\psi}\left(\zeta,k^{2}\right)=\mathcal{O}(1),$

uniformly in this limit. On applying Theorem I of [21, § 6] with $u=\kappa$,
$\alpha=\widetilde{\sigma}$ and $\zeta_{2}=\zeta_{*}$, we obtain the following
solutions of (4.9)

$\displaystyle
W_{1}\left(\zeta,k^{2}\right)=U\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)+\epsilon_{1}\left(\zeta,k^{2}\right),$
$\displaystyle
W_{2}\left(\zeta,k^{2}\right)=\overline{U}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)+\epsilon_{2}\left(\zeta,k^{2}\right),$
(4.13)

valid when $\zeta\in[0,\zeta_{*})$. Choosing as a consequence of (4.12)

$\displaystyle\Omega(z)=z,$

and thus defining the variational operator as

$\displaystyle\mathcal{V}_{a,b}\big{(}\widetilde{\psi}\big{)}=\int_{a}^{b}\frac{\big{|}\widetilde{\psi}\left(t,k^{2}\right)\big{|}}{\sqrt{2\kappa}t}dt,$

the bounds obtained are

$\displaystyle\left|\epsilon_{1}\left(\zeta,k^{2}\right)\right|\leq\frac{\textbf{M}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}{\textbf{E}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}\left[\exp\left(\textstyle{{1}\over{2}}\displaystyle\sqrt{\pi/\kappa}l_{1}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2}\right)\mathcal{V}_{\zeta,\zeta_{*}}\big{(}\widetilde{\psi}\big{)}\right)-1\right],$
$\displaystyle\frac{\left|\epsilon_{2}\left(\zeta,k^{2}\right)\right|}{\textbf{E}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}\leq\textbf{M}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)\left[\exp\left(\textstyle{{1}\over{2}}\displaystyle\sqrt{\pi/\kappa}l_{1}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2}\right)\mathcal{V}_{0,\zeta}\big{(}\widetilde{\psi}\big{)}\right)-1\right],$

where

$\displaystyle
l_{1}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2}\right)=\sup_{z\in(0,\infty)}\left\\{\frac{\Omega(z)\textbf{M}^{2}\left(-\frac{1}{2}\kappa\widetilde{\sigma}^{2},z\right)}{\Gamma\left(\frac{1}{2}\left(1+\kappa\widetilde{\sigma}^{2}\right)\right)}\right\\}.$

The functions E, M and later N are defined in [21, § 5.8]. It follows from
(4.12) and the evenness of $\widetilde{\psi}\left(\zeta,k^{2}\right)$ that

$\displaystyle\mathcal{V}_{0,\zeta}\big{(}\widetilde{\psi}\big{)}=\mathcal{O}\big{(}\kappa^{-1/2}\big{)}\qquad\textrm{and}\qquad\mathcal{V}_{\zeta,\zeta_{*}}\big{(}\widetilde{\psi}\big{)}=\mathcal{O}\big{(}\kappa^{-1/2}\big{)}\qquad\textrm{as}\quad\kappa\to\infty.$

From [21, § 5.8] we have

$\displaystyle\textbf{M}^{2}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},z\right)=\mathcal{O}\left(z^{-1}\right)\qquad\textrm{as}\quad
z\to\infty,$ (4.14)

thus clearly $l_{1}\left(-\frac{1}{2}\kappa\widetilde{\sigma}^{2}\right)$ is
bounded as $\kappa\to\infty$. Hence

$\displaystyle\epsilon_{1}\left(\zeta,k^{2}\right)=\frac{\textbf{M}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}{\textbf{E}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}\mathcal{O}\left(\kappa^{-1}\right),$
$\displaystyle\epsilon_{2}\left(\zeta,k^{2}\right)=\textbf{E}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)\textbf{M}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)\mathcal{O}\left(\kappa^{-1}\right).$
(4.15)

Applying similar analysis to the above, again from Theorem 1 in [21, § 6], we
obtain

$\displaystyle\frac{\partial\epsilon_{1}\left(\zeta,k^{2}\right)}{\partial\zeta}=\frac{\textbf{N}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}{\textbf{E}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)}\mathcal{O}\big{(}\kappa^{-1/2}\big{)},$
$\displaystyle\frac{\partial\epsilon_{2}\left(\zeta,k^{2}\right)}{\partial\zeta}=\textbf{E}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)\textbf{N}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)\mathcal{O}\big{(}\kappa^{-1/2}\big{)}.$
(4.16)

We can extend this analysis to include the point $\zeta_{*}$ since $\zeta$ and
$\psi\left(\zeta,k^{2}\right)$ are bounded there.

###  5 An interlude: eigenvalues

Thus we obtain from (4.13), (4.4) and (4.1) that for $z\in[0,K]$

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=C_{\nu}^{m}\left(\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\right)^{1/4}\left(W_{\nu,1}^{m}\left(\zeta,k^{2}\right)+\eta_{\nu,c}^{m}W_{\nu,2}^{m}\left(\zeta,k^{2}\right)\right),$
$\displaystyle\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)=S_{\nu}^{m+1}\left(\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\right)^{1/4}\left(W_{\nu,1}^{m}\left(\zeta,k^{2}\right)+\eta_{\nu,s}^{m+1}W_{\nu,2}^{m}\left(\zeta,k^{2}\right)\right),$

where $C_{\nu}^{m}$, $S_{\nu}^{m+1}$, $\eta_{\nu,c}^{m}$ and
$\eta_{\nu,s}^{m+1}$ are constants, and $s_{\nu}^{m}$, $\sigma_{\nu}^{m}$,
$W_{\nu,1}^{m}$ and $W_{\nu,2}^{m}$ are defined with respect to either
$h=a_{\nu}^{m}$ or $h=b_{\nu}^{m+1}$ in the previous section.

Let’s first consider $m$ odd. In correspondence with the boundary conditions
we require

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(0,k^{2}\right)=\frac{d\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)}{dz}\bigg{|}_{z=K}=0,$
(5.1)
$\displaystyle\operatorname{Es}_{\nu}^{m+1}\left(0,k^{2}\right)=\operatorname{Es}_{\nu}^{m+1}\left(K,k^{2}\right)=0.$
(5.2)

The requirement at $z=K$ in (5.1) gives

$\displaystyle\textstyle{{1}\over{2}}\displaystyle\zeta_{*}W_{\nu,1}^{m}\left(\zeta_{*},k^{2}\right)+\big{(}\zeta_{*}^{2}-\left(\sigma_{\nu}^{m}\right)^{2}\big{)}\frac{dW_{\nu,1}^{m}\left(\zeta,k^{2}\right)}{d\zeta}\bigg{|}_{\zeta=\zeta_{*}}$
$\displaystyle\qquad{}+\eta_{\nu,c}^{m}\left(\textstyle{{1}\over{2}}\displaystyle\zeta_{*}W_{\nu,2}^{m}\left(\zeta_{*},k^{2}\right)+\big{(}\zeta_{*}^{2}-\left(\sigma_{\nu}^{m}\right)^{2}\big{)}\frac{dW_{\nu,2}^{m}\left(\zeta,k^{2}\right)}{d\zeta}\bigg{|}_{\zeta=\zeta_{*}}\right)=0,$

thus from (4.16), (4.15) and (4.13), and since as $\kappa\to\infty$ (see [24,
§ 12.9])

$\displaystyle
U\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta_{*}\right)\sim
e^{-\frac{\kappa\zeta_{*}^{2}}{2}}\left(\sqrt{2\kappa}\zeta_{*}\right)^{\textstyle{{1}\over{2}}\displaystyle\left(\kappa\widetilde{\sigma}^{2}-1\right)},$
$\displaystyle
U^{\prime}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta_{*}\right)\sim-\sqrt{\frac{\kappa}{2}}\zeta_{*}U\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta_{*}\right),$
$\displaystyle\overline{U}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta_{*}\right)\sim\sqrt{\frac{2}{\pi}}\Gamma\left(\textstyle{{1}\over{2}}\displaystyle\left(1+\kappa\widetilde{\sigma}^{2}\right)\right)e^{\frac{\kappa\zeta_{*}^{2}}{2}}\left(\sqrt{2\kappa}\zeta_{*}\right)^{-\textstyle{{1}\over{2}}\displaystyle\left(\kappa\widetilde{\sigma}^{2}+1\right)},$
$\displaystyle\overline{U}^{\prime}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta_{*}\right)\sim\sqrt{\frac{\kappa}{2}}\zeta_{*}\overline{U}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta_{*}\right),$

necessarily $\eta_{\nu,c}^{m}$ is exponentially small in this limit. Note that
here we have assumed that $\kappa\widetilde{\sigma}^{2}=\mathcal{O}(1)$ as
$\kappa\to\infty$, and we show later, in (5.3), that this assumption is
justified.

Similarly the requirement at $z=K$ in (5.2) gives

$\displaystyle
W_{\nu,1}^{m}\left(\zeta_{*},k^{2}\right)+\eta_{\nu,s}^{m+1}W_{\nu,2}^{m}\left(\zeta_{*},k^{2}\right)=0,$

then necessarily $\eta_{\nu,s}^{m+1}$ is exponentially small as
$\kappa\to\infty$.

Hence for both these cases by considering the requirements at $z=0$ we have
the condition

$\displaystyle
U\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2},0\right)+\mathcal{O}\left(\kappa^{-1}\right)=\frac{2^{\left(\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2}-1\right)/4}\sqrt{\pi}}{\Gamma\left(\frac{1}{4}\left(3-\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2}\right)\right)}+\mathcal{O}\left(\kappa^{-1}\right)=0.$

and to satisfy this we require

$\displaystyle\textstyle{{1}\over{2}}\displaystyle\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2}=j+\textstyle{{1}\over{2}}\displaystyle+\mathcal{O}\left(\kappa^{-1}\right),$

where $j$ is an odd integer.

Let’s now consider $m$ even. In correspondence with the boundary conditions we
require

$\displaystyle\frac{d\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)}{dz}\bigg{|}_{z=0}=\frac{d\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)}{dz}\bigg{|}_{z=K}=0,$
$\displaystyle\frac{d\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)}{dz}\bigg{|}_{z=0}=\operatorname{Es}_{\nu}^{m+1}\left(K,k^{2}\right)=0.$

By similar reason to the $m$ odd case, when considering the boundary condition
at $z=K$ both $\eta_{\nu,c}^{m}$ and $\eta_{\nu,s}^{m+1}$ are necessarily
exponentially small as $\kappa\to\infty$. Hence by using the boundary
condition at $z=0$ we have the condition

$\displaystyle
U^{\prime}\left(-\textstyle{{1}\over{2}}\displaystyle\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2},0\right)+\mathcal{O}\left(\kappa^{-1}\right)=-\frac{2^{\left(\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2}+1\right)/4}\sqrt{\pi}}{\Gamma\left(\frac{1}{4}\left(1-\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2}\right)\right)}+\mathcal{O}\left(\kappa^{-1}\right)=0.$

To satisfy this we require that

$\displaystyle\textstyle{{1}\over{2}}\displaystyle\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2}=j+\textstyle{{1}\over{2}}\displaystyle+\mathcal{O}\left(\kappa^{-1}\right),$

where $j$ is an even integer.

As $\kappa\to\infty$, the zeros of
$U\big{(}{-}\frac{1}{2}\kappa\left(\widetilde{\sigma}_{\nu}^{m}\right)^{2},\sqrt{2\kappa}\zeta\big{)}$
tend to the zeros of $D_{j}\big{(}\sqrt{2\kappa}\zeta\big{)}$, 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}\big{(}\sqrt{2\kappa}\zeta\big{)}$, we deduce
that $j=m$ and as such from (4.7) and (4.10) we have that

$\displaystyle\left(\sigma_{\nu}^{m}\right)^{2}=\frac{2m+1}{\kappa}+\mathcal{O}\left(\kappa^{-2}\right).$
(5.3)

We deduce from (4.7) and (4.3) that

$\displaystyle a_{\nu}^{m}=(2m+1)\kappa+\mathcal{O}(1),\qquad
b_{\nu}^{m+1}=(2m+1)\kappa+\mathcal{O}(1),\qquad\kappa\to\infty.$ (5.4)

###  6 Approximations in terms of parabolic cylinder $\boldsymbol{D}$
functions

The Lamé functions decay exponentially on either side of the oscillatory
interval in $(-K,K)$. Our approximant in Section 4,
$U\left(-\frac{1}{2}\kappa\widetilde{\sigma}^{2},\sqrt{2\kappa}\zeta\right)$,
displays the appropriate exponentially decreasing behaviour when $\zeta$ is
large and positive, but when the variable is large and negative it becomes
exponentially large. If our argument
$-\frac{1}{2}\kappa\widetilde{\sigma}^{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

$\displaystyle\left({\omega}_{\nu}^{m}\right)^{2}=\left({\sigma}_{\nu}^{m}\right)^{2}-\frac{2m+1}{\kappa}=\mathcal{O}\left(\kappa^{-2}\right)\qquad(\kappa\to\infty),$
(6.1)

thus it makes sense to split up (4.8) so that

$\displaystyle\frac{d^{2}W_{\nu}^{m}\left(\zeta,k^{2}\right)}{d\zeta^{2}}=\left(\kappa^{2}\left(\zeta^{2}-\frac{2m+1}{\kappa}\right)+\widehat{\psi}_{\nu}^{m}\left(\zeta,k^{2}\right)\right)W_{\nu}^{m}\left(\zeta,k^{2}\right),$
(6.2)

where

$\displaystyle\widehat{\psi}_{\nu}^{m}\left(\zeta,k^{2}\right)=\psi_{\nu}^{m}\left(\zeta,k^{2}\right)-\kappa^{2}\left({\omega}_{\nu}^{m}\right)^{2}.$

From (6.1) and since $s_{\nu}^{m}$, $\sigma_{\nu}^{m}$ and the variables $x$
and $\zeta$ are all bounded as $\kappa\to\infty$, we have for
$\zeta\in[-\zeta_{*},\zeta_{*}]$

$\displaystyle\widehat{\psi}_{\nu}^{m}\left(\zeta,k^{2}\right)=\mathcal{O}(1),$

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 $\zeta$. We obtain the solutions for (6.2)

$\displaystyle
W_{\nu,1}^{m}\left(\zeta,k^{2}\right)=D_{m}\left(\sqrt{2\kappa}\zeta\right)+\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right),$
$\displaystyle
W_{\nu,2}^{m}\left(\zeta,k^{2}\right)=\overline{D}_{m}\left(\sqrt{2\kappa}\zeta\right)+\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right),$
(6.3)

valid when $\zeta\in[0,\zeta_{*}]$, where we denoted for convenience

$\displaystyle\overline{D}_{m}(t)=\overline{U}\left(-m-\textstyle{{1}\over{2}}\displaystyle,t\right).$

We have that

$\displaystyle\widehat{\psi}_{\nu}^{m}\left(0,k^{2}\right)=-\kappa^{2}\left({\omega}_{\nu}^{m}\right)^{2}+\psi_{\nu}^{m}\left(0,k^{2}\right)=-\kappa^{2}\left({\sigma}_{\nu}^{m}\right)^{2}+(2m+1)\kappa+\frac{\left({s}_{\nu}^{m}\right)^{4}-\left({\sigma}_{\nu}^{m}\right)^{4}}{2\left({\sigma}_{\nu}^{m}\right)^{2}\left({s}_{\nu}^{m}\right)^{4}}$

In Part II, Section 9, we will show that

$\displaystyle\left({\sigma}_{\nu}^{m}\right)^{2}=\frac{2m+1}{\kappa}-\frac{1+k^{2}}{8\kappa^{2}}+\mathcal{O}\left(\kappa^{-3}\right).$

We combine this with (4.11) and obtain

$\displaystyle\widehat{\psi}_{\nu}^{m}\left(0,k^{2}\right)=\mathcal{O}\left(\kappa^{-1}\right)\qquad(\kappa\to\infty).$
(6.4)

This time in the error analysis we choose it so that

$\displaystyle\Omega(z)=1+z,$

and, hence, take as the variational operator

$\displaystyle\mathcal{V}_{a,b}\left(\widetilde{\psi}_{\nu}^{m}\right)=\int_{a}^{b}\frac{\big{|}\widetilde{\psi}_{\nu}^{m}\left(t,k^{2}\right)\big{|}}{1+\sqrt{2\kappa}t}dt.$

The bounds for the errors are

$\displaystyle\left|\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right)\right|\leq\frac{\textbf{M}_{m}\left(\sqrt{2\kappa}\zeta\right)}{\textbf{E}_{m}\left(\sqrt{2\kappa}\zeta\right)}\left[\exp\left(\textstyle{{1}\over{2}}\displaystyle\sqrt{\pi/\kappa}l_{m,1}\mathcal{V}_{\zeta,\zeta_{*}}\big{(}\widetilde{\psi}_{\nu}^{m}\big{)}\right)-1\right],$
$\displaystyle\left|\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right)\right|\leq\textbf{E}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}\textbf{M}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}\left[\exp\left(\textstyle{{1}\over{2}}\displaystyle\sqrt{\pi/\kappa}l_{m,1}\mathcal{V}_{0,\zeta}\big{(}\widetilde{\psi}_{\nu}^{m}\big{)}\right)-1\right],$

where

$\displaystyle
l_{m,1}=\sup_{z\in(0,\infty)}\left\\{\frac{\Omega(z)\textbf{M}^{2}_{m}(z)}{\Gamma(m+1)}\right\\}.$

We define for convenience

$\displaystyle\textbf{E}\left(-m-\textstyle{{1}\over{2}}\displaystyle,z\right)=\textbf{E}_{m}(z),\qquad\textbf{M}\left(-m-\textstyle{{1}\over{2}}\displaystyle,z\right)=\textbf{M}_{m}(z),\qquad\textbf{N}\left(-m-\textstyle{{1}\over{2}}\displaystyle,z\right)=\textbf{N}_{m}(z).$

We have the bounds

$\displaystyle\mathcal{V}_{a,b}\big{(}\widetilde{\psi}_{\nu}^{m}\big{)}=\int_{a}^{b}\frac{\big{|}\widetilde{\psi}_{\nu}^{m}\big{(}t,k^{2}\big{)}\big{|}}{1+\sqrt{2\kappa}t}dt\leq\int_{a}^{b}\frac{\big{|}\widetilde{\psi}_{\nu}^{m}\big{(}0,k^{2}\big{)}\big{|}}{1+\sqrt{2\kappa}t}dt+\int_{a}^{b}\frac{\big{|}\widetilde{\psi}_{\nu}^{m}\big{(}t,k^{2}\big{)}-\widetilde{\psi}_{\nu}^{m}\big{(}0,k^{2}\big{)}\big{|}}{\sqrt{2\kappa}t}dt,$

thus from (6.4) we have

$\displaystyle\mathcal{V}_{0,\zeta}\big{(}\widehat{\psi}_{\nu}^{m}\big{)}=\mathcal{O}\big{(}\kappa^{-1/2}\big{)}\qquad\textrm{and}\qquad\mathcal{V}_{\zeta,\zeta_{*}}\big{(}\widehat{\psi}_{\nu}^{m}\big{)}=\mathcal{O}\big{(}\kappa^{-1/2}\big{)}\qquad\textrm{as}\quad\kappa\to\infty.$

From (4.14) clearly $l_{m,1}$ is a bounded constant as $\kappa\to\infty$.
Hence we have

$\displaystyle\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right)=\frac{\textbf{M}_{m}\left(\sqrt{2\kappa}\zeta\right)}{\textbf{E}_{m}\left(\sqrt{2\kappa}\zeta\right)}\mathcal{O}\left(\kappa^{-1}\right),\qquad\frac{\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right)}{\textbf{E}_{m}\left(\sqrt{2\kappa}\zeta\right)}=\textbf{M}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}\mathcal{O}\left(\kappa^{-1}\right).$
(6.5)

Applying similar analysis to the above, again from Theorem 1 in [21, § 6], we
obtain

$\displaystyle\frac{\partial\epsilon_{1}\left(\zeta,k^{2}\right)}{\partial\zeta}=\frac{\textbf{N}_{m}\left(\sqrt{2\kappa}\zeta\right)}{\textbf{E}_{m}\left(\sqrt{2\kappa}\zeta\right)}\mathcal{O}\big{(}\kappa^{-1/2}\big{)},\qquad\frac{{\partial\epsilon_{2}\left(\zeta,k^{2}\right)}/{\partial\zeta}}{\textbf{E}_{m}\left(\sqrt{2\kappa}\zeta\right)}=\textbf{N}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}\mathcal{O}\big{(}\kappa^{-1/2}\big{)}.$

Thus from (6.3), (4.4) and (4.1) we obtain

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=C_{\nu}^{m}\left(\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\right)^{1/4}$
$\displaystyle\hphantom{\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=}{}\times\left(D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right)+\eta_{\nu,c}^{m}\left(\overline{D}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right)\right)\right),$
$\displaystyle\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)=S_{\nu}^{m+1}\left(\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\right)^{1/4}$
$\displaystyle\hphantom{\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)=}{}\times\left(D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right)+\eta_{\nu,s}^{m+1}\left(\overline{D}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right)\right)\right),$

valid as $\kappa\to\infty$ for $z\in[0,K]$, where
$\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right)$ and
$\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right)$ are given in (6.5), and
$\eta_{\nu,c}^{m}$ and $\eta_{\nu,s}^{m+1}$ are constants to be determined.
Since we have the boundary conditions

$\displaystyle\frac{d\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)}{dz}\bigg{|}_{z=K}=0,\qquad\operatorname{Es}_{\nu}^{m+1}\left(K,k^{2}\right)=0,$

we obtain that necessarily both $\eta_{\nu,c}^{m}$ and $\eta_{\nu,s}^{m+1}$
are $\mathcal{O}\big{(}e^{-\kappa\zeta_{*}^{2}}{\kappa}^{m+1/2}\big{)}$ as
$\kappa\to\infty$. Only when $\zeta$ nears the endpoint the contribution from
the
$\overline{D}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\epsilon_{\nu,2}^{m}\left(\zeta,k^{2}\right)$
term is comparable to
$D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\epsilon_{\nu,1}^{m}\left(\zeta,k^{2}\right)$.
These functions will be made unique by their normalisation.

#### A second term in the approximation in terms of parabolic cylinder
$\boldsymbol{D}$ functions

In [21, § 11] uniform asymptotic expansions are included. The two term version
is of the form

$\displaystyle W_{\nu,1}^{m}\left(\zeta,k^{2}\right)\sim
A_{0}D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\frac{B_{0}(\zeta,\kappa)}{\kappa^{2}}\frac{d}{d\zeta}D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)},$

in which we have taken $A_{0}=1$ and for $B_{0}(\zeta,\kappa)$ we have

$\displaystyle
B_{0}(\zeta,\kappa)=\int_{\tilde{\sigma}}^{\zeta}\frac{\widehat{\psi}_{\nu}^{m}\left(t,k^{2}\right)}{2\sqrt{\left(\zeta^{2}-{\tilde{\sigma}}^{2}\right)\left(t^{2}-{\tilde{\sigma}}^{2}\right)}}dt,$

where $\tilde{\sigma}=\sqrt{(2m+1)/\kappa}$. This coefficient depends on
$\kappa$ and its dominant part is

$\displaystyle
B_{0}(\zeta)=\lim_{\kappa\to\infty}B_{0}(\zeta,\kappa)=\frac{1}{2\zeta}\int_{0}^{\zeta}\frac{\lim\limits_{\kappa\to\infty}\widehat{\psi}_{\nu}^{m}\left(t,k^{2}\right)}{t}dt.$
(6.6)

To simplify the integrand we take the final equation in (4.6) and let
$\kappa\to\infty$. The result is

$\displaystyle\textstyle{{1}\over{2}}\displaystyle\zeta^{2}=\int_{0}^{x}\frac{tdt}{\sqrt{\left(1-t^{2}\right)\left(1-k^{2}t^{2}\right)}}=\frac{1}{k}\ln\left(\frac{1+k}{\sqrt{1-k^{2}x^{2}}+k\sqrt{1-x^{2}}}\right).$

This relation can be inverted

$\displaystyle
x^{2}=k^{-1}\sinh\left(k\zeta^{2}\right)-\left(1+k^{-2}\right)\sinh^{2}\left(k\zeta^{2}/2\right).$

We use this in

$\displaystyle\lim_{\kappa\to\infty}\widehat{\psi}_{\nu}^{m}\left(\zeta,k^{2}\right)=\frac{3}{4\zeta^{2}}+\frac{\zeta^{2}}{4}\left(\frac{k^{2}\left(x^{2}+x^{4}\right)+x^{2}-3}{x^{4}}\right)+\frac{1+k^{2}}{8},$

and are able to evaluate the integral in (6.6) and obtain

$\displaystyle 32\zeta
B_{0}(\zeta)=\left(k^{2}+1\right)\ln\left(\textstyle{{1}\over{4}}\displaystyle\zeta^{2}C(\zeta,k)\right)-\frac{3\left(k^{2}-1\right)^{2}}{2C(\zeta,k)}+3k\coth\left(k\zeta^{2}/2\right)+2k^{2}\zeta^{2}-\frac{6}{\zeta^{2}},$

where $C(\zeta,k)=2k\coth\left(k\zeta^{2}/2\right)-k^{2}-1$. Thus our two term
approximation is

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)\sim
C_{\nu}^{m}\left(\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\right)^{1/4}\left(D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}+\frac{B_{0}(\zeta)}{\kappa^{2}}\frac{d}{d\zeta}D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}\right),$

as $\kappa\to\infty$ and similarly for
$\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)$.

#### Normalisation constant

We now want an approximation for $C_{\nu}^{m}$ and $S_{\nu}^{m+1}$. Due to the
complicated nature of the mapping between $x$ and $\zeta$, it is not simple to
express one in terms of the other. With respect to (3.1) we consider the
integral

$\displaystyle\left(C_{\nu}^{m}\right)^{2}\int_{-K}^{K}\operatorname{dn}(z,k)\sqrt{\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}}D^{2}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}dz.$

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
$\kappa$, it will be the same for both $C_{\nu}^{m}$ and $S_{\nu}^{m+1}$.

First we perform the transformation $x=\operatorname{sn}(z,k)$ to obtain

$\displaystyle\left(C_{\nu}^{m}\right)^{2}\int_{-1}^{1}\sqrt{\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{\big{(}x^{2}-\left(s_{\nu}^{m}\right)^{2}\big{)}\left(1-x^{2}\right)}}D^{2}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}dx,$

then another transformation to the $\zeta$ variable to obtain

$\displaystyle\left(C_{\nu}^{m}\right)^{2}\int_{-\zeta_{*}}^{\zeta^{*}}\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\sqrt{1-k^{2}x^{2}}D^{2}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}d\zeta.$

As the oscillatory behaviour of $D_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}$
happens in a shrinking region of the origin as $\kappa\to\infty$, we seek to
approximate the integral around this point to get an approximation for
$C_{\nu}^{m}$. Thus we consider an expansion of the form

$\displaystyle x=\sum_{k=1}^{\infty}c_{k}\zeta^{k},$ (6.7)

as $\zeta\to 0$, and substituting this into the relation

$\displaystyle\frac{dx}{d\zeta}=\sqrt{\frac{\big{(}\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}\big{)}\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}}.$
(6.8)

From this we can determine, by matching (6.7) and (6.8), the $c_{k}$ terms,
and obtain

$\displaystyle
x=\frac{\sigma_{\nu}^{m}}{s_{\nu}^{m}}\zeta+\left(\frac{(\sigma_{\nu}^{m})^{4}-(s_{\nu}^{m})^{4}}{6(s_{\nu}^{m})^{5}\sigma_{\nu}^{m}}-\frac{\left(1+k^{2}\right)(\sigma_{\nu}^{m})^{3}}{6(s_{\nu}^{m})^{3}}\right)\zeta^{3}+\mathcal{O}\left(\zeta^{5}\right),$

as $\zeta\to 0$. From this we obtain

$\displaystyle\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\sqrt{1-k^{2}x^{2}}=\frac{\left(\sigma_{\nu}^{m}\right)^{2}}{\left(s_{\nu}^{m}\right)^{2}}+\frac{2\left(\sigma_{\nu}^{m}\right)^{4}-2\left(s_{\nu}^{m}\right)^{4}-k^{2}\left(s_{\nu}^{m}\right)^{2}\left(\sigma_{\nu}^{m}\right)^{4}}{2\left(s_{\nu}^{m}\right)^{6}}\zeta^{2}+\mathcal{O}\left(\zeta^{4}\right),$

as $\zeta\to 0$. With respect to this, we consider as a first approximation
for $C_{\nu}^{m}$ the integral

$\displaystyle\left(C_{\nu}^{m}\right)^{2}\int_{-\zeta_{*}}^{\zeta_{*}}\left(\frac{\left(\sigma_{\nu}^{m}\right)^{2}}{\left(s_{\nu}^{m}\right)^{2}}+\frac{\left(2\sigma_{\nu}^{m}\right)^{4}-2\left(s_{\nu}^{m}\right)^{4}-k^{2}\left(s_{\nu}^{m}\right)^{2}\left(\sigma_{\nu}^{m}\right)^{4}}{2\left(s_{\nu}^{m}\right)^{6}}\zeta^{2}\right)D^{2}_{m}\big{(}\sqrt{2\kappa}\zeta\big{)}d\zeta.$

and letting $t=\sqrt{2\kappa}\zeta$, we have the approximation

$\displaystyle\frac{\left(C_{\nu}^{m}\right)^{2}}{\sqrt{2\kappa}}\int_{-\infty}^{\infty}\left(1+\frac{\left(1+k^{2}\right)(2m+1)+t^{2}-k^{2}t^{2}}{8\kappa}\right)D^{2}_{m}(t)dt\sim\frac{\left(C_{\nu}^{m}\right)^{2}m!}{\sqrt{\kappa/\pi}}\left(1+\frac{2m+1}{4\kappa}\right).$

In correspondence with (3.2) and that for large $z$,

$\displaystyle D_{m}(z)>0,\qquad D_{m}^{\prime}(z)<0\qquad\textrm{(see
\cite[cite]{[\@@bibref{}{NIST:DLMF}{}{}, \S~{}12.9(i)]}},$

we deduce from (3.1) that as $\kappa\to\infty$

$\displaystyle\left.\begin{array}[]{@{}l@{}}C_{m}\\\ S_{m+1}\\\
\end{array}\right\\}\sim\frac{\left(\pi\kappa\right)^{{1}/{4}}}{\sqrt{2m!}}\left(1-\frac{2m+1}{8\kappa}\right).$

###  7 A special case: Mathieu functions

We obtain rigorous uniform results for Mathieu’s equation, using the limiting
arguments given in (1.3) and (1.4). Thus (reader can check the details, we
will just summarise) we have the uniform approximations as $h\to\infty$

$\displaystyle\operatorname{ce}_{m}(h,z)=C_{m}\left(\frac{\zeta^{2}-\sigma_{m}^{2}}{x^{2}-s_{m}^{2}}\right)^{1/4}$
$\displaystyle\hphantom{\operatorname{ce}_{m}(h,z)=}{}\times\left(D_{m}\big{(}2\sqrt{h}\zeta\big{)}+\epsilon_{m,1}(\zeta)+\eta_{m}^{c}\big{(}\overline{D}_{m}\big{(}\sqrt{2}h\zeta\big{)}+\epsilon_{m,2}(\zeta)\big{)}\right),$
$\displaystyle\operatorname{se}_{m+1}(h,z)=S_{m+1}\left(\frac{\zeta^{2}-\left(\sigma_{\nu}^{m}\right)^{2}}{x^{2}-\left(s_{\nu}^{m}\right)^{2}}\right)^{1/4}$
$\displaystyle\hphantom{\operatorname{se}_{m+1}(h,z)=}{}\times\left(D_{m}\big{(}2\sqrt{h}\zeta\big{)}+\epsilon_{m,1}(\zeta)+\eta_{m+1}^{s}\big{(}\overline{D}_{m}\big{(}\sqrt{2}h\zeta\big{)}+\epsilon_{m,2}(\zeta)\big{)}\right),$

valid for $\zeta\in[0,\zeta_{*}]$, i.e., $z\in[0,\frac{\pi}{2}]$, as
$h\to\infty$, where

$\displaystyle\epsilon_{m,1}(\zeta)=\frac{\textbf{M}_{m}\left(\sqrt{2}h\zeta\right)}{\textbf{E}_{m}\left(\sqrt{2}h\zeta\right)}\mathcal{O}\left(h^{-1}\right),\qquad\epsilon_{m,2}(\zeta)=\textbf{E}_{m}\big{(}\sqrt{2}h\zeta\big{)}\textbf{M}_{m}\big{(}\sqrt{2}h\zeta\big{)}\mathcal{O}\left(h^{-1}\right)$

both $\eta_{m}^{c}$ and $\eta_{m+1}^{s}$ are
$\mathcal{O}\big{(}e^{-2h\zeta_{*}^{2}}h^{m+1/2}\big{)}$ as $h\to\infty$, the
relationship between $z$ and $x$ is defined by $x=\cos{z}$ and by evaluating
the integrals (4.6), the relationship between $x$ and $\zeta$ is defined by

$\displaystyle\int_{x}^{-s_{m}}\sqrt{\frac{t^{2}-s_{m}^{2}}{(1-t^{2})\left(1-k^{2}t^{2}\right)}}dt=\int_{\zeta}^{-\sigma_{m}}\sqrt{\tau^{2}-\sigma_{m}^{2}}d\tau,\qquad$
$\displaystyle-1<x\leq-s_{m},$
$\displaystyle\int_{-s_{m}}^{x}\sqrt{\frac{s_{m}^{2}-t^{2}}{(1-t^{2})\left(1-k^{2}t^{2}\right)}}dt=\int_{-\sigma_{m}}^{\zeta}\sqrt{\sigma_{m}^{2}-\tau^{2}}d\tau,\qquad$
$\displaystyle-s_{m}\leq x\leq s_{m},$
$\displaystyle\int_{s_{m}}^{x}\sqrt{\frac{t^{2}-s_{m}^{2}}{(1-t^{2})\left(1-k^{2}t^{2}\right)}}dt=\int_{\sigma_{m}}^{\zeta}\sqrt{\tau^{2}-\sigma_{m}^{2}}d\tau,\qquad$
$\displaystyle s_{m}\leq x<1,$

where

$\displaystyle
s^{2}_{m}=\frac{\lambda_{m}+2h^{2}}{4h^{2}},\qquad\sigma_{m}^{2}=s^{2}_{m}\,{}_{2}F_{1}\left({\frac{1}{2},\frac{1}{2}\atop
2};s_{m}^{2}\right),$

and $\lambda_{m}$ corresponds to either $a_{m}$ or $b_{m+1}$ depending on the
solution we are considering, and

$\displaystyle\left.\begin{array}[]{@{}l@{}}C_{m}\\\
S_{m+1}\end{array}\right\\}\sim\left(\frac{\pi
h}{2m!^{2}}\right)^{1/4}\left(1-\frac{2m+1}{16h}\right).$

We have that

$\displaystyle
s_{m}^{2}=\frac{m+\frac{1}{2}}{h}+\mathcal{O}\left(h^{-2}\right)\qquad\textrm{as}\quad
h\to\infty,$

thus from (4.3) and (1.4) we obtain

$\displaystyle a_{m}=-2h^{2}+(4m+2)h+\mathcal{O}(1),\qquad
b_{m+1}=-2h^{2}+(4m+2)h+\mathcal{O}(1),$

as $h\to\infty$. One should note that this result about the eigenvalues was
given formally in [17, § 2.33, Satz 5 & Satz 6].

Two term approximations are given by

$\displaystyle\operatorname{ce}_{m}(h,z)\sim
C_{m}\left(D_{m}\big{(}2\sqrt{h}\zeta\big{)}+\frac{B_{0}(\zeta)}{h^{2}}\frac{d}{d\zeta}D_{m}\big{(}2\sqrt{h}\zeta\big{)}\right),$
$\displaystyle\operatorname{se}_{m+1}(h,z)\sim
S_{m+1}\left(D_{m}\big{(}2\sqrt{h}\zeta\big{)}+\frac{B_{0}(\zeta)}{h^{2}}\frac{d}{d\zeta}D_{m}\big{(}2\sqrt{h}\zeta\big{)}\right),$

where for $B_{0}(\zeta)$ we have the relation

$\displaystyle
256B_{0}(\zeta)=\frac{3\zeta}{4-\zeta^{2}}-\frac{2}{\zeta}\ln\left(1-\frac{\zeta^{2}}{4}\right).$

##  Part II 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

$\displaystyle w_{\nu}^{m}(t)\sim
D_{m}(t)\sum_{s=0}^{\infty}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{\infty}\frac{B_{s}(t)}{\kappa^{s}},\qquad
h\sim(2m+1)\kappa+\sum_{s=0}^{\infty}\frac{\mu_{s}}{\kappa^{s}},$

where $t=\sqrt{2\kappa}\operatorname{sn}{(z,k)}$, the coefficients $A_{s}(t)$
and $B_{s}(t)$ are just polynomials, and the $\mu_{s}$ are constants in terms
of $m$ and $k$. The function expansions will clearly only make sense when
$t=\mathcal{O}(1)$ as $\kappa\to\infty$. 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.

###  8 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
$\kappa\to\infty$, and it can be shown that around the origin $\zeta$ behaves
approximately like $\operatorname{sn}(z,k)$. Thus the variable in the
parabolic cylinder function around this point behaves approximately like
$\sqrt{2\kappa}\operatorname{sn}(z,k)$. This motivates the next simpler
transformation.

####  $\boldsymbol{t}$-plane

Letting $t=\sqrt{2\kappa}\operatorname{sn}(z,k)$ in (1.1) we have

$\displaystyle\frac{d^{2}w}{dt^{2}}+\frac{2h-t^{2}\kappa}{4\kappa}w-\frac{1}{2\kappa}\left(\left(\left(1+k^{2}\right)t^{2}-\frac{k^{2}t^{4}}{2\kappa}\right)\frac{d^{2}}{dt^{2}}+t\left(1+k^{2}-\frac{k^{2}t^{2}}{\kappa}\right)\frac{d}{dt}\right)w=0,$

where $z\in(-K,K)$ corresponds to
$t\in\big{(}{-}\sqrt{2\kappa},\sqrt{2\kappa}\big{)}$. Supposing in accordance
with (5.4) that

$\displaystyle\frac{h}{2\kappa}=m+\textstyle{{1}\over{2}}\displaystyle+\sum_{s=1}^{n}\frac{\mu_{s}}{\kappa^{s}}+\frac{\widetilde{\mu}_{n+1}}{\kappa^{n+1}},$
(8.1)

where $n$ is a positive integer and $\widetilde{\mu}_{n}$ can be re-expanded
in a sensible manner. We can then write Lamé’s equation in the form

$\displaystyle\frac{d^{2}w}{dt^{2}}+\left(m+\frac{1}{2}-\frac{t^{2}}{4}\right)w+\frac{1}{\kappa}\left\\{\left(-\frac{1}{2}\left(1+k^{2}\right)t^{2}+\frac{k^{2}t^{4}}{4\kappa}\right)\frac{d^{2}}{dt^{2}}\right.$
$\displaystyle\left.\qquad{}+\frac{t}{2}\left(-1-k^{2}+\frac{k^{2}t^{2}}{\kappa}\right)\frac{d}{dt}+\sum_{s=0}^{n}\frac{\mu_{s+1}}{\kappa^{s}}+\frac{\widetilde{\mu}_{n+1}}{\kappa^{n+1}}\right\\}w=0.$
(8.2)

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

$\displaystyle
w_{\nu}^{m}\left(t,k^{2}\right)=D_{m}(t)\sum_{s=0}^{\infty}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{\infty}\frac{B_{s}(t)}{\kappa^{s}}$
(8.3)

seems sensible. However one should observe that this splitting only makes
sense when $t=\mathcal{O}(1)$ as $h\to\infty$, 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\in\mathbb{Z}$ instead of in terms
of $D_{m}(t)$ and $D_{m}^{\prime}(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 $\kappa$, we have the
recurrence relations for $A_{s}(t)$ and $B_{s}(t)$

$\displaystyle
2A_{s}^{\prime}(t)+B_{s}^{\prime\prime}(t)-\frac{t+tk^{2}}{2}\left(2tA_{s-1}^{\prime}(t)+A_{s-1}(t)+tB_{s-1}^{\prime\prime}(t)+B_{s-1}^{\prime}(t)+\tilde{t}(m)tB_{s-1}(t)\right)$
$\displaystyle\qquad{}+\frac{k^{2}t^{3}}{4}\left(2tA_{s-2}^{\prime}(t)+2A_{s-2}(t)+tB_{s-2}^{\prime\prime}(t)+2B_{s-2}^{\prime}(t)+\tilde{t}(m)tB_{s-2}(t)\right)$
$\displaystyle\qquad{}+\sum_{j=1}^{s}\mu_{j}B_{s-j}(t)=0,$ (8.4)
$\displaystyle
A_{s}^{\prime\prime}(t)+2\tilde{t}(m)B_{s}^{\prime}(t)+\frac{t}{2}B_{s}(t)+\sum_{j=1}^{s}\mu_{j}A_{s-j}(t)-\frac{t+tk^{2}}{2}\bigl{(}tA_{s-1}^{\prime\prime}(t)+A_{s-1}^{\prime}(t)+\frac{t^{2}}{2}B_{s-1}(t)$
$\displaystyle\qquad{}+\tilde{t}(m)\left(tA_{s-1}(t)+2tB_{s-1}^{\prime}(t)+B_{s-1}(t)\right)\bigr{)}+\frac{k^{2}t^{3}}{4}\bigl{(}tA_{s-2}^{\prime\prime}(t)+2A_{s-2}^{\prime}(t)+\frac{t^{2}}{2}B_{s-2}(t)$
$\displaystyle\qquad{}+\tilde{t}(m)\left(tA_{s-2}(t)+2tB_{s-2}^{\prime}(t)+2B_{s-2}(t)\right)\bigr{)}=0,$
(8.5)

where $\tilde{t}(m)=\frac{1}{4}t^{2}-m-\frac{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^{\prime\prime}_{s}(t)$ and substitute it into (8.5); this gives the third
order inhomogeneous differential equation for $B_{s}(t)$

$\displaystyle
B_{s}^{\prime\prime\prime}(t)-\left(t^{2}-4m-2\right)B_{s}^{\prime}(t)-tB_{s}(t)=b_{s}(t),$
(8.6)

where

$\displaystyle
b_{s}(t)=\sum_{j=1}^{s}\mu_{j}\left(2A_{s-j}(t)-B_{s-j}^{\prime}(t)\right)$
$\displaystyle\hphantom{b_{s}(t)=}{}+\frac{1+k^{2}}{2}\left[3tA_{s-1}^{\prime}(t)+A_{s-1}(t)+t^{2}B_{s-1}^{\prime\prime\prime}(t)+3tB_{s-1}^{\prime\prime}(t)+B_{s-1}^{\prime}(t)\right.$
$\displaystyle\left.\hphantom{b_{s}(t)=}{}-\textstyle{{1}\over{2}}\displaystyle
t^{3}B_{s-1}(t)-\tilde{t}(m)t^{2}\left(2A_{s-1}(t)+3B_{s-1}^{\prime}(t)\right)\right]$
$\displaystyle\hphantom{b_{s}(t)=}{}-\frac{k^{2}t^{2}}{2}\left[3tA_{s-2}^{\prime}(t)+3A_{s-2}(t)+\textstyle{{1}\over{2}}\displaystyle
t^{2}B_{s-2}^{\prime\prime\prime}(t)+3tB_{s-2}^{\prime\prime}(t)+3B_{s-2}^{\prime}(t)\right.$
$\displaystyle\left.\hphantom{b_{s}(t)=}{}-\textstyle{{1}\over{4}}\displaystyle
t^{3}B_{s-2}(t)-\textstyle{{1}\over{2}}\displaystyle
t^{2}\tilde{t}(m)\left(2A_{s-2}(t)+3B_{s-2}^{\prime}(t)\right)\right].$ (8.7)

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 coefficients $\boldsymbol{B_{s}(t)}$ and $\boldsymbol{A_{s}(t)}$

Using variation of parameters we obtain the general solution for $B_{s}(t)$

$\displaystyle
B_{s}(t)=b^{1}_{s}(t)D^{2}_{m}(t)+b^{2}_{s}(t)\overline{D}^{2}_{m}(t)+b^{3}_{s}(t)D_{m}(t)\overline{D}_{m}(t),$

where $\overline{D}_{m}(t)=\overline{U}\left(-m-\frac{1}{2},t\right)$ (for
$\overline{U}$ see formula (12.2.21) in [24]), and

$\displaystyle
b^{1}_{s}(t)=\int\frac{b_{s}(t)\mathcal{W}\\{\overline{D}^{2},D\overline{D}\\}}{\mathcal{W}\\{D^{2},\overline{D}^{2},D\overline{D}\\}}dt+c_{1},\qquad
b^{2}_{s}(t)=-\int\frac{b_{s}(t)\mathcal{W}\\{D^{2},D\overline{D}\\}}{\mathcal{W}\\{D^{2},\overline{D}^{2},D\overline{D}\\}}dt+c_{2},$
$\displaystyle
b^{3}_{s}(t)=\int\frac{b_{s}(t)\mathcal{W}\\{D^{2},\overline{D}^{2}\\}}{\mathcal{W}\\{D^{2},\overline{D}^{2},D\overline{D}\\}}dt+c_{3}.$

By expanding the Wronskians we derive the relations

$\displaystyle\mathcal{W}\\{D^{2},D\overline{D}\\}=\mathcal{W}\\{D,\overline{D}\\}D^{2}_{m}(t),\qquad$
$\displaystyle\mathcal{W}\\{\overline{D}^{2},D\overline{D}\\}=-\mathcal{W}\\{D,\overline{D}\\}\overline{D}^{2}_{m}(t),$
$\displaystyle\mathcal{W}\\{D^{2},\overline{D}^{2}\\}=2\mathcal{W}\\{D,\overline{D}\\}D_{m}(t)\overline{D}_{m}(t),\qquad$
$\displaystyle\mathcal{W}\\{D^{2},\overline{D}^{2},D\overline{D}\\}=-2\mathcal{W}\\{D,\overline{D}\\}^{3}.$

Then without loss of generality we can rewrite the indefinite integrals in
$\left\\{b^{i}_{s}(t)\right\\}_{i=1}^{3}$ 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

$\displaystyle
B_{s}(t)=\frac{\pi}{4m!^{2}}\int_{0}^{t}b_{s}(\tau)\left(\overline{D}_{m}(\tau)D_{m}(t)-{D}_{m}(\tau)\overline{D}_{m}(t)\right)^{2}d\tau+c_{3}D_{m}(t)\overline{D}_{m}(t).$
(8.8)

Although we consider $t$ to be $\mathcal{O}(1)$ as $\kappa\to\infty$, 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_{m}^{2}(t)$
or ${D}_{m}(t)\overline{D}_{m}(t)$ grow no faster than polynomials when $t$
becomes large. Since

$\displaystyle\overline{D}_{m}(t)^{2}\sim\frac{2m!^{2}}{\pi}e^{t^{2}/2}t^{-2m-2},\qquad
t\to\infty,$

to ensure the boundedness of our formal expansion as $t$ becomes large,
$\mu_{s}$ is determined uniquely by the condition that

$\displaystyle\int_{0}^{\infty}b_{s}(\tau){D}_{m}^{2}(\tau)\,d\tau=0.$ (8.9)

This condition will also ensure boundedness as $t\to-\infty$. Note that in the
next section we will derive an alternative method to compute the $\mu_{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

$\displaystyle B_{0}(t)=c_{3}D_{m}(t)\overline{D}_{m}(t),$

and then

$\displaystyle A_{0}(t)=-\textstyle{{1}\over{2}}\displaystyle
c_{3}\big{(}D_{m}^{\prime}(t)\overline{D}_{m}(t)+D_{m}(t)\overline{D}_{m}^{\prime}(t)\big{)}+c.$

Rearranging, on the $\kappa^{0}$ level of the asymptotic expansion (8.3) we
have

$\displaystyle\left(\tfrac{1}{2}c_{3}\mathcal{W}(D_{m},\overline{D}_{m})+c\right)D_{m}(t).$

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\geq 1$. With these choices we have

$\displaystyle
b_{1}(t)=2\mu_{1}+\frac{1+k^{2}}{4}\left(2-t^{2}\left(t^{2}-4m-2\right)\right),$

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

$\displaystyle\mu_{1}=-\textstyle{{1}\over{8}}\displaystyle\left(1+k^{2}\right)\left(1+2m+2m^{2}\right).$

In the case $s=1$, it can be shown using integration by parts that for
$\mu_{1}$ which satisfies (8.9)

$\displaystyle\frac{\pi}{4m!^{2}}\int_{0}^{t}b_{1}(\tau)\big{(}\overline{D}_{m}(\tau)D_{m}(t)-{D}_{m}(\tau)\overline{D}_{m}(t)\big{)}^{2}d\tau$
$\displaystyle\qquad{}=\frac{1+k^{2}}{16}\left(t^{3}-(2m+1)t+(-1)^{m}\frac{2m+1}{m!}\sqrt{\frac{\pi}{2}}{D}_{m}(t)\overline{D}_{m}(t)\right).$

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}\overline{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}\overline{D}_{m}$, one would see that it would appear for all $s\geq
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 coefficients $\boldsymbol{B_{s}(t)}$ and
$\boldsymbol{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\geq 1$ we consider $A_{s}(t)$ and
$B_{s}(t)$ in the form

$\displaystyle
A_{s}(t)=\sum_{i=1}^{\infty}a_{s,i}t^{2i}\qquad\textrm{and}\qquad
B_{s}(t)=\sum_{i=0}^{\infty}b_{s,i}t^{2i+1}.$

Substituting these into (8.6) and (8.4) we obtain the recurrence relations for
coefficients

$\displaystyle(2i+3)(2i+2)(2i+1)b_{s,i+1}+(4m+2)(2i+1)b_{s,i}-2ib_{s,i-1}$
$\displaystyle\qquad{}+2\sum_{j=0}^{s-1}\mu_{s-j}\left(\left(i+\textstyle{{1}\over{2}}\displaystyle\right)b_{j,i}-a_{j,i}\right)-\left(1+k^{2}\right)\Bigl{[}\left(3i+\textstyle{{1}\over{2}}\displaystyle\right)a_{s-1,i}+\textstyle{{1}\over{2}}\displaystyle\left(2i+1\right)^{3}b_{s-1,i}$
$\displaystyle\qquad{}+\left(m+\textstyle{{1}\over{2}}\displaystyle\right)\left(a_{s-1,i-1}+\left(3i-\textstyle{{3}\over{2}}\displaystyle\right)b_{s-1,i-1}\right)-\textstyle{{1}\over{4}}\displaystyle
a_{s-1,i-2}-\left(\textstyle{{3}\over{4}}\displaystyle
i-\textstyle{{7}\over{8}}\displaystyle\right)b_{s-1,i-2}\Bigr{]}$
$\displaystyle\qquad{}+k^{2}\Bigl{[}\left(3i-\textstyle{{3}\over{2}}\displaystyle\right)a_{s-2,i-1}\\!+2i\left(i^{2}\\!-\textstyle{{1}\over{4}}\displaystyle\right)b_{s-2,i-1}\\!+(2m+1)\left(\textstyle{{1}\over{4}}\displaystyle
a_{s-2,i-2}\\!+\left(\textstyle{{3}\over{4}}\displaystyle
i-\textstyle{{9}\over{8}}\displaystyle\right)b_{s-2,i-2}\right)$
$\displaystyle\qquad{}-\textstyle{{1}\over{8}}\displaystyle
a_{s-2,i-3}-\left(\textstyle{{3}\over{8}}\displaystyle
i-\textstyle{{13}\over{16}}\displaystyle\right)b_{s-2,i-3}\Bigr{]}=0,$ (8.10)
$\displaystyle
2(2i+2)a_{s,i+1}+(2i+3)(2i+2)b_{s,i+1}+\sum_{j=0}^{s-1}\mu_{s-j}b_{j,i}$
$\displaystyle\qquad{}-\textstyle{{1}\over{2}}\displaystyle\left(1+k^{2}\right)\left((4i+1)a_{s-1,i}+(2i+1)^{2}b_{s-1,i}-\left(m+\textstyle{{1}\over{2}}\displaystyle\right)b_{s-1,i-1}+\textstyle{{1}\over{4}}\displaystyle
b_{s-1,i-2}\right)$ (8.11)
$\displaystyle\qquad{}+k^{2}\left(i\left(i-\textstyle{{1}\over{2}}\displaystyle\right)b_{s-2,i-1}+\left(i-\textstyle{{1}\over{2}}\displaystyle\right)a_{s-2,i-1}-\textstyle{{1}\over{8}}\displaystyle(2m+1)b_{s-2,i-2}+\textstyle{{1}\over{16}}\displaystyle
b_{s-2,i-3}\right)=0.$

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.

Table 3: The orders of the polynomials. $s$ |  $A_{s}(t)$ |  $B_{s}(t)$
---|---|---
even |  $4s$ |  $4s-3$
odd |  $4s-2$ |  $4s-1$

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

$\displaystyle\mu_{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 $\mu_{s}$; in this manner the
eigenvalue terms are determined uniquely. These $\mu_{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
$\mu_{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 $\boldsymbol{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
$\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)$. We want to construct a
similar expansions which when transformed back into the $z$-plane behave like
the odd Lamé functions $\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)$. We
need our expansion to have the property

$\displaystyle
w_{\nu}^{m}\left(-K,k^{2}\right)=w_{\nu}^{m}\left(K,k^{2}\right)=0.$

The Jacobi elliptic function $\operatorname{cn}(z,k)$ is even around the
origin and odd around $z=-K$ and $z=K$. Hence again letting
$t=\sqrt{2\kappa}\operatorname{sn}(z,k)$ we consider a formal expansion for a
solution of (8.2) of the form

$\displaystyle
w_{\nu}^{m}\left(z,k^{2}\right)=\operatorname{cn}(z,k)\left(D_{m}(t)\sum_{s=0}^{\infty}\frac{P_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=1}^{\infty}\frac{Q_{s}(t)}{\kappa^{s}}\right).$

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

$\displaystyle\operatorname{cn}\left(\operatorname{arcsn}\left(\frac{t}{\sqrt{2\kappa}},k\right),k\right)=\sqrt{1-\frac{t^{2}}{2\kappa}},$

we can express the formal solution in the form

$\displaystyle
w^{m}_{\nu}\left(z,k^{2}\right)=\sqrt{1-\frac{t^{2}}{2\kappa}}\left(D_{m}(t)\sum_{s=0}^{\infty}\frac{P_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=1}^{\infty}\frac{Q_{s}(t)}{\kappa^{s}}\right).$

Expanding this square root, we can rewrite this formal expansion as

$\displaystyle
w_{\nu}^{m}\left(z,k^{2}\right)=D_{m}(t)\sum_{s=0}^{\infty}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=1}^{\infty}\frac{B_{s}(t)}{\kappa^{s}},$

where the connection between $P_{s}(t)$ and $A_{s}(t)$, and $Q_{s}(t)$ and
$B_{s}(t)$ is given by

$\displaystyle
A_{s}(t)=\sum_{j=0}^{s}\binom{\frac{1}{2}}{j}\left(-\frac{t^{2}}{2}\right)^{j}P_{s-j}(t),\qquad
B_{s}(t)=\sum_{j=0}^{s}\binom{\frac{1}{2}}{j}\left(-\frac{t^{2}}{2}\right)^{j}Q_{s-j}(t),$

where $\binom{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_{\nu}^{m}$
and $b_{\nu}^{m+1}$.

###  9 An interlude: error analysis for the eigenvalue expansions

In this section we match our eigenvalue expansions with the eigenvalues
$a_{\nu}^{m}$ and $b_{\nu}^{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

$\displaystyle\frac{d}{dx}\left[p(x)\frac{dy}{dx}\right]+\left(\lambda
w(x)-q(x)\right)y=0$

and we consider $p(x),w(x)>0$, and $p(x)$, $p^{\prime}(x)$, $q(x)$ and $w(x)$
to be continuous functions over a finite real interval $[a,b]$. When $(a,b)$
is bounded and $p(x)$ does not vanish on $[a,b]$, this is a regular SL
problem, otherwise it is singular. Consider the regular problem. We want to
find special values of $\lambda$ called eigenvalues for which there exists a
non-trivial solution satisfying the separated boundary conditions

$\displaystyle\alpha_{1}y(a)+\alpha_{2}y^{\prime}(a)=0,\qquad$
$\displaystyle\alpha_{1}^{2}+\alpha_{2}^{2}>0,$
$\displaystyle\beta_{1}y(b)+\beta_{2}y^{\prime}(b)=0,\qquad$
$\displaystyle\beta_{1}^{2}+\beta_{2}^{2}>0.$ (9.1)

The regular SL theory states that these eigenvalues are real and can be
ordered such that

$\displaystyle\lambda_{0}<\lambda_{1}<\cdots<\lambda_{\nu}<\cdots\to\infty$

and these eigenvalues $\lambda_{\nu}$ correspond to unique eigenfunctions
$y_{\nu}(x)$ with exactly $\nu$ zeros in $(a,b)$. Normalised eigenfunctions
form an orthonormal basis such that

$\displaystyle(y_{n},y_{m})=\int_{a}^{b}y_{n}(x)y_{m}(x)w(x)dx=\delta_{mn}$

in the Hilbert space $L^{2}\left[[a,b],w(x)dx\right]$.

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}\left[[a,b],w(x)dx\right]$ containing all the twice
differentiable functions which satisfy the boundary conditions (9.1) and
corresponding ${\lambda}_{\nu}$ such that

$\displaystyle({\rm SL}+{\lambda}_{\nu})y_{\nu}(x)=0.$

Let $\widetilde{y}(x)\in A$ be such that $(\widetilde{y},\widetilde{y})=1$.
Now define a remainder function $R(x)$ such that

$\displaystyle({\rm SL}+{\lambda})\widetilde{y}(x)=R(x)$

for some constant ${\lambda}$. If $R(x)\in L^{2}\left[[a,b],w(x)dx\right]$
then

$\displaystyle\min_{k}|{\lambda}-{\lambda}_{k}|<\sqrt{(R,R)}.$

###### Proof.

We have that

$\displaystyle(R,y_{k})=((SL+{\lambda})\widetilde{y},y_{k})=(\widetilde{y},SLy_{k})+{\lambda}(\widetilde{y},y_{k})=({\lambda}-{\lambda}_{k})(\widetilde{y},y_{k}),$

and using this in

$\displaystyle
1=(\widetilde{y},\widetilde{y})=\sum_{k}|(\widetilde{y},y_{k})|^{2}=\sum_{k}\frac{|(R,y_{k})|^{2}}{|{\lambda}-{\lambda}_{k}|^{2}}<\frac{(R,R)}{\min_{k}|{\lambda}-{\lambda}_{k}|^{2}},$

we have that

$\displaystyle\min_{k}|{\lambda}-{\lambda}_{k}|<\sqrt{(R,R)}.$ ∎

#### Eigenvalues of Lamé’s equation

We now apply this theory to Lamé’s equation to obtain order estimates upon
truncation of asymptotic expansions of the eigenvalues $a_{\nu}^{m}$ and
$b_{\nu}^{m+1}$.

Special eigenvalues $a_{\nu}^{m}$ correspond to non-trivial solutions
satisfying the separated boundary conditions

$\displaystyle\frac{d}{dz}\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)\bigg{|}_{z=-K}=\frac{d}{dz}\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)\bigg{|}_{z=K}=0$

and these eigenvalues are real and can be ordered such that

$\displaystyle a^{0}_{\nu}<a^{1}_{\nu}<\cdots<a^{m}_{\nu}<\cdots\to\infty.$

The functions $\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)$ have exactly
$m$ zeros in $(-K,K)$. The normalised eigenfunctions
$w_{\nu}^{m}\left(z,k^{2}\right)$ form an orthonormal basis such that

$\displaystyle\left(w_{\nu}^{i},w_{\nu}^{j}\right)=\int_{-K}^{K}w_{\nu}^{i}\left(z,k^{2}\right)w_{\nu}^{j}\left(z,k^{2}\right)dz=\delta_{ij}$

in the Hilbert space $L^{2}\left[[-K,K],dz\right]$.

Let $L_{\nu}$ be the operator

$\displaystyle
L_{\nu}:=\frac{d^{2}}{dz^{2}}-\kappa^{2}\operatorname{sn}^{2}(z,k)$

so that

$\displaystyle\left(L_{\nu}+{a}_{\nu}^{m}\right)w_{\nu}^{m}\left(z,k^{2}\right)=0.$

Letting $t=\sqrt{2\kappa}\operatorname{sn}(z,k)$, we define the truncated
expansions corresponding to $w_{\nu}^{m}\left(z,k^{2}\right)$ as

$\displaystyle
w_{\nu,n}^{m}\left(z,k^{2}\right)=c_{\nu,n}^{m}\left(D_{m}(t)\sum_{s=0}^{n}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{B_{s}(t)}{\kappa^{s}}\right),$
(9.2)

where the $A_{s}(t)$ and $B_{s}(t)$ terms were derived previously, and
$c_{\nu,n}^{m}$ is defined to be a function of $\kappa$ so that

$\displaystyle\int_{-K}^{K}\left(w_{\nu,n}^{m}\left(z,k^{2}\right)\right)^{2}dz=1.$
(9.3)

Then we can write the derivative of $w_{\nu,n}^{m}\left(z,k^{2}\right)$ with
respect to $z$ as

$\displaystyle\frac{dw_{\nu,n}^{m}\left(z,k^{2}\right)}{dz}=c_{\nu,n}^{m}\sqrt{2\kappa}\operatorname{cn}(z,k)\operatorname{dn}(z,k)$
$\displaystyle\hphantom{\frac{dw_{\nu,n}^{m}\left(z,k^{2}\right)}{dz}=}{}\times\left(D_{m}(t)\sum_{s=0}^{n}\frac{{A}^{\prime}_{s}(t)+\left(\frac{t^{2}}{4}-m-\frac{1}{2}\right)B_{s}(t)}{\kappa^{s}}+D^{\prime}_{m}(t)\sum_{s=0}^{n}\frac{A_{s}(t)+{B}^{\prime}_{s}(t)}{\kappa^{s}}\right),$

where the dash represents differentiation with respect to $t$. Thus clearly
$w_{\nu,n}^{m}\left(z,k^{2}\right)\in A$ since
$\operatorname{cn}(-K,k)=\operatorname{cn}(K,k)=0$. We also define the
truncated eigenvalue expansions

$\displaystyle
a_{\nu,n}^{m}=(2m+1)\kappa+2\sum_{s=0}^{n-1}\frac{\mu_{s+1}}{\kappa^{s}},$

where the $\mu_{s}$ were derived previously and define
$R_{\nu,n}^{m}\left(z,k^{2}\right)$ such that

$\displaystyle\left(L_{\nu}+{a}_{\nu,n}^{m}\right)w_{\nu,n}^{m}\left(z,k^{2}\right)=R_{\nu,n}^{m}\left(z,k^{2}\right).$

We consider the operator $L_{\nu}$ acting on $D_{m}(t)$ and derive

$\displaystyle
L_{\nu}\left(D_{m}(t)\right)=\left((2m+1)\left(\frac{1+k^{2}}{2}t^{2}-\kappa\right)-\left(k^{2}(2m+1)+\kappa\left(1+k^{2}\right)\right)\frac{t^{4}}{4\kappa}+\frac{k^{2}t^{6}}{8\kappa}\right)D_{m}(t)$
$\displaystyle\hphantom{L_{\nu}\left(D_{m}(t)\right)=}{}+\left(-\left(1+k^{2}\right)t+\frac{k^{2}t^{3}}{\kappa}\right)D_{m}^{\prime}(t).$
(9.4)

Using the recurrence relations

$\displaystyle tD_{m}(t)=D_{m+1}(t)+mD_{m-1}(t),\qquad
D^{\prime}_{m}(t)=mD_{m-1}(t)-\textstyle{{1}\over{2}}\displaystyle tD_{m}(t),$
(9.5)

we can rewrite $w_{\nu,n}^{m}\left(z,k^{2}\right)$ given in (9.2) such that
for $s\in\\{0,\dots,n\\}$, $A_{s}(t)D_{m}(t)$ and $B_{s}(t)D_{m}^{\prime}(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

$\displaystyle
R_{\nu,n}^{m}\left(z,k^{2}\right)=c_{\nu,n}^{m}\kappa^{-n}\left[(\cdots)+\kappa^{-1}(\cdots)+\cdots\right],$

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_{\nu,n}^{m}\left(z,k^{2}\right)\in L^{2}\left[[-K,K],dz\right]$. As
$w_{\nu,n}^{m}\left(z,k^{2}\right)\in A$ we have that

$\displaystyle\min_{k}\big{|}{a}_{\nu,n}^{m}-{a}_{\nu}^{k}\big{|}^{2}<\sqrt{\left(R_{\nu,n}^{m},R_{\nu,n}^{m}\right)}.$

In Part I, Section 5, we proved that

$\displaystyle
a_{\nu}^{m}=(2m+1)\kappa+\mathcal{O}(1)\qquad\textrm{as}\quad\kappa\to\infty,$

hence for $\kappa$ large enough, necessarily

$\displaystyle\min_{k}\big{|}{a}_{\nu,n}^{m}-{a}_{\nu}^{k}\big{|}^{2}=\big{|}{a}_{\nu,n}^{m}-{a}_{\nu}^{m}\big{|}^{2}$

and thus

$\displaystyle\big{|}{a}_{\nu,n}^{m}-{a}_{\nu}^{m}\big{|}<\sqrt{\left(R_{\nu,n}^{m},R_{\nu,n}^{m}\right)}.$

We need an order estimate for ${\left(R_{\nu,n}^{m},R_{\nu,n}^{m}\right)}$.
The integrals we must consider then are of the form

$\displaystyle
I=\int_{-K}^{K}D_{i}\big{(}\sqrt{2\kappa}\operatorname{sn}(z,k)\big{)}D_{j}\big{(}\sqrt{2\kappa}\operatorname{sn}(z,k)\big{)}dz,$

where $i,j\in\mathbb{N}_{0}$. Performing the substitution
$t=\sqrt{2\kappa}\operatorname{sn}(z,k)$ we obtain

$\displaystyle
I=\frac{1}{\sqrt{2\kappa}}\int_{-\sqrt{2\kappa}}^{\sqrt{2\kappa}}\frac{1}{\sqrt{1-\frac{t^{2}}{2\kappa}}\sqrt{1-\frac{k^{2}t^{2}}{2\kappa}}}D_{i}(t)D_{j}(t)dt.$

Since when $t$ is large $D_{m}$ is exponentially small it follows that

$\displaystyle
I=\mathcal{O}\big{(}\kappa^{-\frac{1}{2}}\big{)},\qquad\textrm{as}\quad\kappa\to\infty.$

Considering the first term of the expansion for
$w_{\nu,n}^{m}\left(z,k^{2}\right)$ in (9.2) we have

$\displaystyle\left(c_{\nu,n}^{m}\right)^{2}\int_{-K}^{K}D_{m}^{2}\big{(}\sqrt{2\kappa}\operatorname{sn}(z,k)\big{)}dz=\frac{\left(c_{\nu,n}^{m}\right)^{2}}{\sqrt{2\kappa}}\int_{-\sqrt{2\kappa}}^{\sqrt{2\kappa}}\frac{1}{\sqrt{1-\frac{t^{2}}{2\kappa}}\sqrt{1-\frac{k^{2}t^{2}}{2\kappa}}}D_{m}^{2}(t)dt$
$\displaystyle\hphantom{\left(c_{\nu,n}^{m}\right)^{2}\int_{-K}^{K}D_{m}^{2}\big{(}\sqrt{2\kappa}\operatorname{sn}(z,k\big{)})dz}{}\sim\frac{\left(c_{\nu,n}^{m}\right)^{2}}{\sqrt{2\kappa}}\int_{-\infty}^{\infty}D_{m}^{2}(t)\,dt=\left(c_{\nu,n}^{m}\right)^{2}m!\sqrt{\frac{\pi}{\kappa}},$

as $\kappa\to\infty$. Hence it follows from (9.3) that
$c_{m,n}=\mathcal{O}(\kappa^{1/4})$ as $\kappa\to\infty$. Similar observations
would give us that

$\displaystyle\int_{-K}^{K}\left(R_{m,n}(z)\right)^{2}dz=\mathcal{O}\left(\kappa^{-2n}\right),$

and so

$\displaystyle{a}_{\nu}^{m}-{a}_{\nu,n}^{m}=\mathcal{O}\left(\kappa^{-n}\right),$

as $\kappa\to\infty$. The error analysis for $b_{\nu,n}^{m}$ would in the same
manner give the order estimate

$\displaystyle{b}_{\nu}^{m+1}-{b}_{\nu,n}^{m+1}=\mathcal{O}\left(\kappa^{-n}\right)\qquad\textrm{as}\quad\kappa\to\infty.$

Hence we have from (8.1) that

$\displaystyle
a_{\nu}^{m}=(2m+1)\kappa+2\sum_{s=0}^{n-1}\frac{\mu_{s+1}}{\kappa^{s}}+\mathcal{O}\left(\kappa^{-n}\right).$

In the same manner we would deduce that

$\displaystyle
b_{\nu}^{m+1}=(2m+1)\kappa+2\sum_{s=0}^{n-1}\frac{\mu_{s+1}}{\kappa^{s}}+\mathcal{O}\left(\kappa^{-n}\right).$

This result also follows from the difference between $a_{\nu}^{m}$ and
$b_{\nu}^{m+1}$ being exponentially small as $\kappa\to\infty$ (see [24, §
28.8]).

###  10 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

$\displaystyle L_{\nu}:=\frac{d^{2}}{dt^{2}}+m+\frac{1}{2}-\frac{t^{2}}{4}$
$\displaystyle\hphantom{L_{\nu}:=}{}+\frac{1}{2\kappa}\left(\left(-t^{2}\left(1+k^{2}\right)+\frac{k^{2}t^{4}}{2\kappa}\right)\frac{d^{2}}{dt^{2}}-t\left(1+k^{2}-\frac{k^{2}t^{2}}{\kappa}\right)\frac{d}{dt}+h-\kappa(2m+1)\right),$

and consider $t\in(-t_{*},t_{*})$. We have the truncated expansion
corresponding to an even solution

$\displaystyle
w_{\nu,n}^{m}\left(t,k^{2}\right)=D_{m}(t)\sum_{s=0}^{n}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=1}^{n}\frac{B_{s}(t)}{\kappa^{s}},$

such that we have the exact solution

$\displaystyle
w_{\nu}^{m}\left(t,k^{2}\right)=w_{\nu,n}^{m}\left(t,k^{2}\right)+\epsilon_{\nu,n}^{m}\left(t,k^{2}\right).$
(10.1)

We define the remainder term $R_{\nu,n}^{m}\left(t,k^{2}\right)$ such that

$\displaystyle
L_{\nu}\left(w_{\nu,n}^{m}\left(t,k^{2}\right)\right)=R_{\nu,n}^{m}\left(t,k^{2}\right)$

and split the eigenvalue such that

$\displaystyle\frac{h}{2\kappa}=m+\textstyle{{1}\over{2}}\displaystyle+\sum_{s=1}^{n}\frac{\mu_{s}}{\kappa^{s}}+\frac{\widetilde{\mu}_{n+1}}{\kappa^{n+1}},$

and note we just proved that $\widetilde{\mu}_{n+1}=\mathcal{O}(1)$ as
$\kappa\to\infty$. Since the coefficients $A_{s}(t)$ and $B_{s}(t)$ satisfy
(8.4) and (8.5) it follows that

$\displaystyle
R_{\nu,n}^{m}\left(t,k^{2}\right)=\mathcal{O}\left(\kappa^{-n-1}\right),$

as $\kappa\to\infty$. Applying $L_{\nu}$ to (10.1) we obtain

$\displaystyle\left(\epsilon_{\nu,n}^{m}\right)^{\prime\prime}+\left(m+\frac{1}{2}-\frac{t^{2}}{4}\right)\epsilon_{\nu,n}^{m}=\frac{1}{1-\frac{t^{2}}{2\kappa}\left(1+k^{2}-\frac{k^{2}t^{2}}{2\kappa}\right)}\biggl{[}\frac{t}{2\kappa}\left(1+k^{2}-\frac{k^{2}t^{2}}{\kappa}\right)\left(\epsilon_{\nu,n}^{m}\right)^{\prime}$
$\displaystyle\qquad{}+\left(t^{2}\left(-1-k^{2}+\frac{k^{2}t^{2}}{2\kappa}\right)\left(2m+1-\frac{t^{2}}{2}\right)-h+\kappa\left(2m+1\right)\right)\frac{\epsilon_{\nu,n}^{m}}{2\kappa}-R_{\nu,n}^{m}\biggr{]},$

and denoting the right hand side of this equation
$\Omega_{\nu,n}^{m}\left(t,k^{2}\right)$, by use of variation of parameters we
have

$\displaystyle\epsilon_{\nu,n}^{m}\left(t,k^{2}\right)=\frac{\sqrt{\pi/2}}{m!}\int_{t}^{t_{*}}\left[D_{m}(t)\overline{D}_{m}(\tau)-D_{m}(\tau)\overline{D}_{m}(t)\right]\Omega_{\nu,n}^{m}\left(\tau,k^{2}\right)d\tau.$

In accordance with [22, § 6.10.2] we define $J(\tau)=1$,
$H(\tau)=1+k^{2}-\frac{k^{2}\tau^{2}}{2\kappa}$ and

$\displaystyle
K\left(t,\tau\right)=\frac{\sqrt{\pi/2}}{m!}\left(D_{m}(t)\overline{D}_{m}(\tau)-D_{m}(\tau)\overline{D}_{m}(t)\right),\qquad\phi(\tau)=\frac{-R_{\nu,n}^{m}\left(\tau,k^{2}\right)}{1-\frac{\tau^{2}}{2\kappa}H(\tau)},$
$\displaystyle\psi_{0}(\tau)=\frac{\kappa\left(2m+1\right)-h-\tau^{2}H(\tau)\left(2m+1-\frac{\tau^{2}}{2}\right)}{2\kappa-\tau^{2}H(\tau)},\qquad\psi_{1}(\tau)=\frac{\tau\left(1+k^{2}-\frac{k^{2}\tau^{2}}{\kappa}\right)}{2\kappa-\tau^{2}H(\tau)},$
$\displaystyle\Phi(t)=\int_{t}^{t_{*}}\left|\phi(\tau)\,d\tau\right|,\qquad\Psi_{0}(t)=\int_{t}^{t_{*}}\left|\psi_{0}(\tau)\,d\tau\right|,\qquad\Psi_{1}(t)=\int_{t}^{t_{*}}\left|\psi_{1}(\tau)\,d\tau\right|.$

Since we consider $t\in(-t_{*},t_{*})$ where $t_{*}=\mathcal{O}(1)$ as
$\kappa\to\infty$, the error analysis is much simpler than the analysis Olver
uses in [21] as we have

$\displaystyle|K\left(t,\tau\right)|\leq
k_{0},\qquad\textrm{and}\qquad|\partial K\left(t,\tau\right)/\partial t|\leq
k_{1},$

where $k_{0}$ and $k_{1}$ are $\mathcal{O}(1)$ as $\kappa\to\infty$. Thus
again in accordance with [22, § 6.10.2] we define

$\displaystyle P_{0}(t)=k_{0},\qquad Q(\tau)=1,\qquad P_{1}(t)=k_{1}$ (10.2)

(we do not define $P_{2}(t)$ as we do not need to bound
$|\partial^{2}K(t,\tau)/\partial t^{2}|$ to carry out our analysis), and
finally the constants

$\displaystyle\widetilde{\kappa}=1,\qquad\widetilde{\kappa}_{0}=k_{0},\qquad\widetilde{\kappa}_{1}=k_{1}.$
(10.3)

Hence it follows from Theorem 10.1 in [22, § 6.10.2] that

$\displaystyle\big{|}\epsilon_{\nu,n}^{m}\big{(}t,k^{2}\big{)}\big{|}\leq
P_{0}(t)\widetilde{\kappa}\Phi(t)\exp\big{[}\widetilde{\kappa}_{0}\Psi_{0}(t)+\widetilde{\kappa}_{1}\Psi_{1}(t)\big{]}.$
(10.4)

Since

$\displaystyle h-\kappa(2m+1)=\mathcal{O}(1)$

for both $h=a_{\nu}^{m}$ and $h=b_{\nu}^{m+1}$ we obtain

$\displaystyle\Phi(t)=\mathcal{O}\left(\kappa^{-n-1}\right),\qquad\Psi_{0}(t)=\mathcal{O}\left(\kappa^{-1}\right),\qquad\Psi_{1}(t)=\mathcal{O}\left(\kappa^{-1}\right),$
(10.5)

as $\kappa\to\infty$. Then substituting expressions from (10.5), (10.3) and
the first of (10.2) into (10.4) we have

$\displaystyle\epsilon_{\nu,n}^{m}\left(t,k^{2}\right)=\mathcal{O}\left(\kappa^{-n-1}\right)\qquad\textrm{as}\quad\kappa\to\infty,\quad\textrm{for}\quad
t\in\left(-t_{*},t_{*}\right).$

###  11 Identification of solutions

Now we identify the solutions derived in Section 10 with the Lamé functions.
We give the identification for $t\in(-t_{*},t_{*})$

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=C_{\nu}^{m}\biggl{(}D_{m}(t)\sum_{s=0}^{n}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{B_{s}(t)}{\kappa^{s}}$
$\displaystyle\hphantom{\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=}{}+\textstyle{{1}\over{2}}\displaystyle\left(\epsilon_{\nu,n}^{m}\left(t,k^{2}\right)+(-1)^{m}\epsilon_{\nu,n}^{m}\left(-t,k^{2}\right)\right)\biggr{)},$
$\displaystyle\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)=S_{\nu}^{m+1}\sqrt{1-\frac{t^{2}}{2\kappa}}\biggl{(}D_{m}(t)\sum_{s=0}^{n}\frac{P_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{Q_{s}(t)}{\kappa^{s}}$
$\displaystyle\hphantom{\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)=}{}+\textstyle{{1}\over{2}}\displaystyle\left(\epsilon_{\nu,n}^{m}\left(t,k^{2}\right)+(-1)^{m}\epsilon_{\nu,n}^{m}\left(-t,k^{2}\right)\right)\biggr{)},$

where the errors are defined according to $h=a_{\nu}^{m}$ or $h=b_{\nu}^{m+1}$
respectively. We consider now just $C_{\nu}^{m}$ since we will obtain the same
asymptotic expansion from $S_{\nu}^{m+1}$ by construction. To obtain an
asymptotic expansion for these constants we consider with respect to (3.1) the
integral

$\displaystyle\left(C_{\nu}^{m}\right)^{2}\int_{-K}^{K}\operatorname{dn}(z,k)\left(D_{m}(t)\sum_{s=0}^{\infty}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{\infty}\frac{B_{s}(t)}{\kappa^{s}}\right)^{2}dz.$

In the integral we let $t=\sqrt{2\kappa}\operatorname{sn}(z,k)$ and obtain

$\displaystyle\frac{\left(C_{\nu}^{m}\right)^{2}}{\sqrt{2\kappa}}\int_{-\sqrt{2\kappa}}^{\sqrt{2\kappa}}\frac{1}{\sqrt{1-\frac{t^{2}}{2\kappa}}}\left(D_{m}(t)\sum_{s=0}^{\infty}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{\infty}\frac{B_{s}(t)}{\kappa^{s}}\right)^{2}dt.$

Since the parabolic cylinder functions are exponentially small when the
variable is large, we consider the integral from $-\infty$ to $\infty$ and
express the integral in the form

$\displaystyle\frac{\left(C_{\nu}^{m}\right)^{2}}{\sqrt{2\kappa}}\int_{-\infty}^{\infty}\sum_{s=0}^{\infty}\kappa^{-s}\sum_{j=0}^{s}\binom{-\frac{1}{2}}{j}\left(t^{2}/2\right)^{j}\sum_{i=0}^{s-j}\biggl{(}A_{i}(t)A_{s-j-i}(t)D_{m}^{2}(t)$
$\displaystyle\qquad{}+2A_{i}(t)B_{s-j-i}(t)D_{m}(t)D_{m}^{\prime}(t)+B_{i}(t)B_{s-j-i}(t)\left(D_{m}^{\prime}(t)\right)^{2}\biggr{)}dt.$

Then from (3.1) we have the formal asymptotic expansions for the normalisation
constants

$\displaystyle\left.\begin{array}[]{@{}l@{}}C_{\nu}^{m}\\\
S_{\nu}^{m+1}\end{array}\right\\}\sim\frac{(\pi\kappa)^{1/4}}{\sqrt{2m!}}\left(1+\sum_{s=1}^{\infty}\frac{\eta_{s}}{\kappa^{s}}\right)^{-1/2}.$

To obtain analytic expressions for the $\eta_{s}$ terms we need expressions
for integrals of the form

$\displaystyle p(m,n)=\int_{-\infty}^{\infty}t^{n}D^{2}_{m}(t)dt,\qquad
q(m,n)=\int_{-\infty}^{\infty}t^{n}\left(D_{m}^{\prime}(t)\right)^{2}dt,$
$\displaystyle
r(m,n)=\int_{-\infty}^{\infty}t^{n+1}D_{m}(t)D_{m}^{\prime}(t)dt.$

We have the identities $p(m,2n+1)=q(m,2n+1)=r(m,2n+1)=0$ for all $m$, $n$, and

$\displaystyle
p(m,0)=\int_{-\infty}^{\infty}D^{2}_{m}(t)dt=m!\sqrt{2\pi}\qquad\mathrm{and}\qquad\int_{-\infty}^{\infty}D_{m}(t)D_{n}(t)dt=0\qquad\mathrm{for}\quad
m\neq n,$

and using (9.5) we deduce the expression $p(m,2)=(2m+1)m!\sqrt{2\pi}$. Using
integration by parts we obtain the recurrence relation for $p(m,n)$

$\displaystyle
p(m,n)=\frac{2(n-1)(2m+1)}{n}p(m,n-2)+\frac{(n-3)(n-2)(n-1)}{n}p(m,n-4).$

We know $p(m,0)$ and $p(m,2)$ thus consequently can determine $p(m,n)$ for any
$n$ recursively. Similarly using integration by parts we obtain expressions
for the other integrals in terms of $p(m,n)$:

$\displaystyle
q(m,n)=\frac{n+3}{4(n+1)}p(m,n+2)-\left(m+\tfrac{1}{2}\right)p(m,n),$
$\displaystyle
r(m,n)=\frac{2m+1}{n+2}p(m,n+2)-\frac{n+4}{2(n+2)(n+3)}p(m,n+4).$

###  12 Summary of results for the Lamé functions
and their respective eigenvalues

Let $\kappa=k\sqrt{\nu(\nu+1)}$ and
$t=\sqrt{2\kappa}\operatorname{sn}{(z,k)}$. Then for $m$ a non-negative
integer and $z=\mathcal{O}\big{(}\kappa^{-1/2}\big{)}$ as $\kappa\to\infty$

$\displaystyle\operatorname{Ec}_{\nu}^{m}\left(z,k^{2}\right)=C_{\nu}^{m}\left(D_{m}(t)\sum_{s=0}^{n}\frac{A_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{B_{s}(t)}{\kappa^{s}}+\mathcal{O}\left(\kappa^{-n-1}\right)\right),$
$\displaystyle\frac{\operatorname{Es}_{\nu}^{m+1}\left(z,k^{2}\right)}{\operatorname{cn}(z,k)}=S_{\nu}^{m+1}\left(D_{m}(t)\sum_{s=0}^{n}\frac{P_{s}(t)}{\kappa^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{Q_{s}(t)}{\kappa^{s}}+\mathcal{O}\left(\kappa^{-n-1}\right)\right),$

where

$\displaystyle\left.\begin{array}[]{@{}l@{}}C_{\nu}^{m}\\\
S_{\nu}^{m+1}\end{array}\right\\}\sim\frac{(\pi\kappa)^{1/4}}{\sqrt{2m!}}\left(1+\sum_{s=1}^{\infty}\frac{\eta_{s}}{\kappa^{s}}\right)^{-1/2}.$

Both $A_{s}(t)$ and $P_{s}(t)$ are even polynomials, and both $B_{s}(t)$ and
$Q_{s}(t)$ are odd polynomials. These polynomials are found recursively and we
give here the first two terms:

$\displaystyle A_{0}=1,\qquad A_{1}=\frac{k^{2}+1}{32}t^{2},\qquad
B_{0}=0,\qquad B_{1}=\frac{k^{2}+1}{16}\left(t^{3}-(2m+1)t\right),$
$\displaystyle P_{0}=1,\qquad P_{1}=\frac{k^{2}+9}{32}t^{2},\qquad
Q_{0}=0,\qquad Q_{1}=B_{1},\qquad\eta_{1}=\frac{3-k^{2}}{16}(2m+1).$

Correspondingly we have the eigenvalue expansion as $\kappa\to\infty$

$\displaystyle\left.\begin{array}[]{@{}l@{}}a_{\nu}^{m}\\\
b_{\nu}^{m+1}\end{array}\right\\}=(2m+1)\kappa+2\sum_{s=0}^{n-1}\frac{\mu_{s+1}}{\kappa^{s}}+\mathcal{O}\left(\kappa^{-n}\right),\qquad\textrm{as}\quad\kappa\to\infty.$

where the order term will be different in both cases. The $\mu_{s}$ terms are
constant coefficients which depend on $k$ and $m$, found recursively with the
eigenfunction expansions. We give here the first two terms:

$\displaystyle\mu_{1}=-\frac{k^{2}+1}{8}\left(1+2m+2m^{2}\right),\qquad\mu_{2}=-\frac{2m+1}{32}\left(\left(k^{2}-1\right)^{2}\left(1+m+m^{2}\right)-2k^{2}\right).$

These eigenvalue coefficients match the formal results given in [24, § 29.7].

###  13 A special case: Mathieu functions

#### Mathieu’s equation

Similarly to the Lamé case, if one considers the two term uniform
approximation for the Mathieu functions you would observe that the oscillatory
behaviour of the Mathieu functions happen in a shrinking neighbourhood of the
$z=\frac{\pi}{2}$ as $h\to\infty$. It can be shown that in a shrinking
neighbourhood of this point, $\zeta$ behaves approximately like $\cos(z)$.
Thus the variable in the parabolic cylinder function around this point behaves
approximately like $\sqrt{2}h\cos(z)$. This would motivate a much simpler
transformation.

Letting $t=2\sqrt{h}\cos{z}$ in (1.2) we present Mathieu’s equation in the
algebraic form

$\displaystyle\frac{d^{2}w}{dt^{2}}+\left(\frac{\lambda+2h^{2}}{4h}-\frac{t^{2}}{4}\right)w-\frac{1}{4h}\left(t^{2}\frac{d^{2}}{dt^{2}}+t\frac{d}{dt}\right)w=0,$

and since we proved that

$\displaystyle\frac{\lambda+2h^{2}}{4}-h\left(m+\textstyle{{1}\over{2}}\displaystyle\right)=\mathcal{O}(1),\qquad\textrm{as}\quad
h\to\infty$

for the special eigenvalues $\lambda$ which emit the Mathieu functions we have
Mathieu’s equation in the form

$\displaystyle\frac{d^{2}w}{dt^{2}}+\left(m+\frac{1}{2}-\frac{t^{2}}{4}\right)w+\frac{1}{h}\left(-\frac{t^{2}}{4}\frac{d^{2}}{dt^{2}}-\frac{t}{4}\frac{d}{dt}+\frac{\lambda+2h^{2}}{4}-h\left(m+\frac{1}{2}\right)\right)w=0.$

In a similar manner to the previous sections coefficients in the expansions
for the functions and eigenvalues can be computed using the same ansatz,
although they would only make sense asymptotically in a shrinking
neighbourhood of $z=\frac{\pi}{2}$. However these are also realised by
considering the Mathieu functions as a special case of the Lamé functions.

#### Summary of results as a special case of Lamé’s equation

We obtain rigorous uniform results for Mathieu’s equation, using the limiting
arguments given in (1.3) and (1.4). Thus (reader can check the details, we
will just summarise) letting $t=2\sqrt{h}\cos{z}$, for $m\geq 0$ and
$z=\frac{\pi}{2}+\mathcal{O}(h^{-1/2})$ we have as $h\to\infty$

$\displaystyle\operatorname{ce}_{m}(h,a_{m},z)=C_{m}\left(D_{m}(t)\sum_{s=0}^{n}\frac{A_{s}(t)}{h^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{B_{s}(t)}{h^{s}}+\mathcal{O}\left(h^{-n-1}\right)\right),$
$\displaystyle\frac{\operatorname{se}_{m+1}(h,b_{m+1},z)}{\sin{z}}=S_{m+1}\left(D_{m}(t)\sum_{s=0}^{n}\frac{P_{s}(t)}{h^{s}}+D_{m}^{\prime}(t)\sum_{s=0}^{n}\frac{Q_{s}(t)}{h^{s}}+\mathcal{O}\left(h^{-n-1}\right)\right),$
(13.1)

where

$\displaystyle\left.\begin{array}[]{@{}l@{}}C_{m}\\\
S_{m+1}\end{array}\right\\}\sim\left(\frac{\pi
h}{2(m!)^{2}}\right)^{1/4}\left(1+\sum_{s=1}^{\infty}\frac{\eta_{s}}{h^{s}}\right)^{-1/2}.$

Both $A_{s}(t)$ and $P_{m+1,s}(t)$ are even polynomials, and both $B_{s}(t)$
and $Q_{m+1,s}(t)$ are odd polynomials. These polynomials are found
recursively and we give here the first few terms:

$\displaystyle A_{0}=1,\qquad A_{1}=\frac{t^{2}}{2^{6}},\qquad
B_{0}=Q_{0}=0,\qquad B_{1}=Q_{1}=\frac{t^{3}}{2^{5}}-(1+2m)\frac{t}{2^{5}},$
$\displaystyle P_{0}=1,\qquad P_{1}=\frac{9t^{2}}{2^{6}},$ $\displaystyle
A_{2}=\frac{t^{8}}{2^{13}}-(1+2m)\frac{t^{6}}{2^{11}}+\left(9+10m+10m^{2}\right)\frac{t^{4}}{2^{12}}+\left(5+6m-12m^{2}-8m^{3}\right)\frac{t^{2}}{2^{12}},$
$\displaystyle
B_{2}=\frac{t^{5}}{2^{8}}-(1+2m)\frac{5t^{3}}{2^{11}}-\left(11+20m+20m^{2}\right)\frac{t}{2^{11}},$
$\displaystyle
P_{2}=\frac{t^{8}}{2^{13}}-(1+2m)\frac{t^{6}}{2^{11}}+\left(113+10m+10m^{2}\right)\frac{t^{4}}{2^{12}}+\left(5+6m-12m^{2}-8m^{3}\right)\frac{t^{2}}{2^{12}},$
$\displaystyle
Q_{2}=\frac{t^{5}}{2^{7}}-(1+2m)\frac{13t^{3}}{2^{11}}-\left(11-20m-20m^{2}\right)\frac{t}{2^{11}},\qquad\eta_{1}=\frac{6m+3}{32}.$

Correspondingly we have the eigenvalue expansion as $h\to\infty$

$\displaystyle\left.\begin{array}[]{@{}l@{}}a_{m}\\\
b_{m+1}\end{array}\right\\}=-2h^{2}+4h\sum_{s=0}^{n}\frac{\mu_{s}}{h^{s}}+\mathcal{O}\left(h^{-n}\right),$
(13.4)

where the order term will be different in both cases. The $\mu_{s}$ terms are
constant coefficients which depend on $m$, found recursively with the
eigenfunction expansions. We give here the first few terms:

$\displaystyle\mu_{0}=m+\textstyle{{1}\over{2}}\displaystyle,\qquad\mu_{1}=-\textstyle{{1}\over{16}}\displaystyle\big{(}1+2m+2m^{2}\big{)},\qquad\mu_{2}=-\textstyle{{1}\over{128}}\displaystyle\left(1+3m+3m^{2}+2m^{3}\right).$

These function and eigenvalue coefficients match the formal results given in
[24, § 28.8] and [17, § 2], and various other papers discussed at the start of
this paper. Note that similar results to (13.1) and (13.4) were given in [13],
but there the function expansions did not make sense in the interval they were
stated to hold in, and the method used there to obtain the coefficients in the
expansions for the functions and eigenvalue was very cumbersome. Our methods
were simple and enabled us to perform rigorous error analysis on these
expansions.

#### Acknowledgements

The authors thank the referees for very helpful comments and suggestions for
improving the presentation.

### References

  * [1] Al-Gwaiz M.A., Sturm–Liouville theory and its applications, Springer Undergraduate Mathematics Series, Springer-Verlag London, Ltd., London, 2008.
  * [2] Arscott F.M., Periodic differential equations. An introduction to Mathieu, Lamé, and allied functions, International Series of Monographs in Pure and Applied Mathematics, Vol. 66, The Macmillan Co., New York, 1964.
  * [3] Dingle R.B., Müller H.J.W., Asymptotic expansions of Mathieu functions and their characteristic numbers, J. Reine Angew. Math. 211 (1962), 11–32.
  * [4] Dunster T.M., Uniform asymptotic approximation of Mathieu functions, Methods Appl. Anal. 1 (1994), 143–168.
  * [5] Erdélyi A., Magnus W., Oberhettinger F., Tricomi F.G., Higher transcendental functions. Vol. III, McGraw-Hill Book Company, Inc., New York – Toronto – London, 1955.
  * [6] Goldstein S., Mathieu functions, Trans. Cambridge Philos. Soc. 23 (1927), 303–336.
  * [7] Goldstein S., On the asymptotic expansion of the characteristic numbers of the Mathieu equation, Proc. Roy. Soc. Edinburgh 49 (1930), 210–223.
  * [8] Ince E.L., Researches into the characteristic numbers of the Mathieu equation, Proc. Roy. Soc. Edinburgh 46 (1927), 20–29.
  * [9] Ince E.L., Researches into the characteristic numbers of the Mathieu equation (second paper), Proc. Roy. Soc. Edinburgh 46 (1927), 316–322.
  * [10] Ince E.L., Researches into the characteristic numbers of the Mathieu equation (third paper), Proc. Roy. Soc. Edinburgh 47 (1928), 294–301.
  * [11] Ince E.L., The periodic Lamé functions, Proc. Roy. Soc. Edinburgh 60 (1940), 47–63.
  * [12] Ince E.L., Further investigations into the periodic Lamé functions, Proc. Roy. Soc. Edinburgh 60 (1940), 83–99.
  * [13] Kurz M., Fehlerabschätzungen zu asymptotischen Entwicklungen der Eigenwerte und Eigenlösungen der Mathieuschen Differentialgleichung, Ph.D. Thesis, Universität Duisburg-Essen, Essen, 1979\.
  * [14] Lamé G., Sur les surfaces isothermes dans les corps homogènes en équilibre de température, J. Math. Pures Appl. 2 (1837), 147–188.
  * [15] Mathieu E., Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique, J. Math. Pures Appl. 13 (1868), 137–203.
  * [16] McLachlan N.W., Theory and application of Mathieu functions, Dover Publications, Inc., New York, 1964.
  * [17] Meixner J., Schäfke F.W., Mathieusche Funktionen und Sphäroidfunktionen mit Anwendungen auf physikalische und technische Probleme, Die Grundlehren der mathematischen Wissenschaften, Vol. 71, Springer-Verlag, Berlin – Göttingen – Heidelberg, 1954.
  * [18] Müller H.J.W., Asymptotic expansions of ellipsoidal wave functions and their characteristic numbers, Math. Nachr. 31 (1966), 89–101.
  * [19] Müller H.J.W., Asymptotic expansions of ellipsoidal wave functions in terms of Hermite functions, Math. Nachr. 32 (1966), 49–62.
  * [20] Müller H.J.W., On asymptotic expansions of ellipsoidal wave functions, Math. Nachr. 32 (1966), 157–172.
  * [21] Olver F.W.J., Second-order linear differential equations with two turning points, Philos. Trans. Roy. Soc. London Ser. A 278 (1975), 137–174.
  * [22] Olver F.W.J., Asymptotics and special functions, AKP Classics, A K Peters, Ltd., Wellesley, MA, 1997.
  * [23] Olver F.W.J., Lozier D.W., Boisvert R.F., Clark C.W. (Editors), NIST handbook of mathematical functions, U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC; Cambridge University Press, Cambridge, 2010, print companion to [24].
  * [24] NIST digital library of mathematical functions, Release 1.0.10 of 2015-08-07, online companion to [23], available at http://dlmf.nist.gov/.
  * [25] Sips R., Représentation asymptotique des fonctions de Mathieu et des fonctions d’onde sphéroidales, Trans. Amer. Math. Soc. 66 (1949), 93–134.
  * [26] Sips R., Représentation asymptotique des fonctions de Mathieu et des fonctions sphéroidales. II, Trans. Amer. Math. Soc. 90 (1959), 340–368.
  * [27] Sips R., Représentation asymptotique de la solution générale de l’équation de Mathieu–Hill, Bull. Cl. Sci. V Sér. 51 (1965), 1415–1446.
  * [28] Weinstein M.I., Keller J.B., Hill’s equation with a large potential, SIAM J. Appl. Math. 45 (1985), 200–214.
  * [29] Whittaker E.T., Watson G.N., A course of modern analysis, Cambridge Mathematical Library, Cambridge University Press, Cambridge, 1996.