On the eigenvalues of the spheroidal wave equation

Harald Schmid [email protected] University of Applied Sciences Amberg-Weiden, Amberg, Germany

Abstract.

This paper presents some new results on the eigenvalues of the spheroidal wave equation. We study the angular and Coulomb spheroidal wave equation as a special case of a more general linear Hamiltonian system depending on three parameters. We prove that the eigenvalues of this system satisfy a first-order quasilinear partial differential equation with respect to the parameters. This relation offers a new insight on how the eigenvalues of the spheroidal wave equation depend on the spheroidal parameter. Apart from analytical considerations, the PDE we obtain can also be used for a numerical computation of spheroidal eigenvalues.

Key words and phrases:

spheroidal eigenvalues, Coulomb spheroidal wave equation, linear Hamiltonian system, inviscid forced Burgers’ equation, deformation method

1991 Mathematics Subject Classification:

33E10, 33F05, 34L15, 35F20

1. Introduction

The angular spheroidal wave equation

\frac{\mathrm{d}}{\mathrm{d}x}\left((1-x^{2})\frac{\mathrm{d}}{\mathrm{d}x}w(x)\right)+\left(\lambda+\gamma^{2}(1-x^{2})-\frac{\mu^{2}}{1-x^{2}}\right)w(x)=0,\quad-1<x<1

(1)

appears in many fields of physics and engineering like quantum mechanics, electromagnetism, signal processing etc. Here $\mu$ is supposed to be some given real number, and since only $\mu^{2}$ occurs in (1), we may assume $\mu\geq 0$ without loss of generality. In many applications the so-called spheroidal parameter (or size parameter) $\gamma$ is either real or purely imaginary, so that $\gamma^{2}$ is a real number which may also be negative. If $\mu$ is an integer and $\gamma^{2}$ is real, then the separation of the Helmholtz equation in prolate ( $\gamma^{2}>0$ ) or oblate ( $\gamma^{2}<0$ ) spheroidal coordinates results in a second order ODE of the form (1). A slightly more general differential equation is the Coulomb spheroidal wave equation (CSWE)

\frac{\mathrm{d}}{\mathrm{d}x}\left((1-x^{2})\frac{\mathrm{d}}{\mathrm{d}x}w(x)\right)+\left(\lambda+\beta x+\gamma^{2}(1-x^{2})-\frac{\mu^{2}}{1-x^{2}}\right)w(x)=0

(2)

It differs from (1) only by the presence of a linear term $\beta x$ with some additional parameter $\beta$ , which we also assume to be real. The numbers $\lambda\in\mathbb{C}$ for which (2) has a nontrivial bounded solution $w(x)$ on $(-1,1)$ are the eigenvalues (or characteristic values) of the CSWE, and the corresponding eigenfunctions $w(x)$ are called Coulomb spheroidal functions. They appear in astrophysics and molecular physics and provide, for example, exact wave functions for a one-electron diatomic molecule with fixed nuclei (see [1, Chapter 9] or [2]). In our subsequent considerations, it does not make much difference whether we include the term $\beta x$ or not. Therefore, in this paper we will mainly deal with equation (2), while the results we obtain are obviously applicable to the angular spheroidal wave equation (1) as well.

For technical reasons we introduce the parameters $u_{1}$ , $u_{2}$ , $u_{3}$ which are related to $\beta$ , $\gamma^{2}$ , $\lambda$ in (2) by

\gamma^{2}=u_{1},\quad\beta=-u_{3}-2(\mu+1)u_{2},\quad\lambda=u_{3}+\mu(\mu+1)

(3)

Moreover, we set $\alpha:=\frac{1}{2}(\mu+1)$ . With these parameters in mind, we will see that the $2\times 2$ differential system

\begin{pmatrix}[r]0&-1\\[4.30554pt] 1&0\end{pmatrix}y^{\prime}(z)-\begin{pmatrix}2&2u_{2}-\frac{\alpha}{z}+\frac{\alpha}{1-z}\\[4.30554pt] 2u_{2}-\frac{\alpha}{z}+\frac{\alpha}{1-z}&2(u_{1}+u_{2}^{2})+\frac{u_{3}}{z}\end{pmatrix}y(z)=\Lambda\begin{pmatrix}\frac{1}{1-z}&0\\[4.30554pt] 0&\frac{1}{z}\end{pmatrix}y(z),\quad 0<z<1

(4)

is related to the Coulomb spheroidal wave equation by the following means: In section 2 we show that (4), in combination with appropriate boundary conditions, can be written as an eigenvalue problem for a self-adjoint differential operator $T=T(u_{1},u_{2},u_{3})$ . The eigenvalues $\Lambda$ of $T$ depend analytically on the parameters $(u_{1},u_{2},u_{3})$ , and in addition they have the following properties:

(a)

$\lambda=u_{3}+\mu(\mu+1)$ is an eigenvalue of the CSWE (2) for the parameters $\gamma^{2}=u_{1}$ , $\beta=-u_{3}-2(\mu+1)u_{2}$ if and only if $\Lambda=0$ is an eigenvalue of the linear Hamiltonian system (4).

(b)

The eigenvalues $\Lambda=\Lambda(u_{1},u_{2},u_{3})$ of (4) satisfy the first-order quasilinear partial differential equation

\begin{split}2u_{1}\frac{\partial\Lambda}{\partial u_{1}}&+\left((\Lambda+2)(u_{1}+u_{2}^{2})+\Lambda+u_{2}+u_{3}\right)\frac{\partial\Lambda}{\partial u_{2}}\\[-4.30554pt] &+\left((2\Lambda+u_{3}+2)(2\Lambda u_{2}-2\mu-1)-2(\mu+1)(u_{1}+u_{2}^{2})+2\mu\right)\frac{\partial\Lambda}{\partial u_{3}}\\ &=(1+2\mu-2\Lambda u_{2})(\Lambda+2)-2\mu\end{split}

(5)

According to (a), the eigenvalue problem (4) may be regarded as a generalization of the Coulomb spheroidal wave equation, and the PDE (5) describes an analytic relation between the eigenvalues and the parameters, where the eigenfunctions need not to be known. A similar observation has already been made in [3] for the Chandrasekhar-Page equation, and [4] provides a method by which one can establish such a relation between the eigenvalues and the parameters in terms of a PDE for other linear Hamiltonian systems as well. In the present paper we will use this \qqdeformation method from [4] to prove assertion (b). Finally, in section 3 we solve (5) by the method of characteristics, where we consider the prolate spheroidal wave equation ( $\beta=0$ , $\gamma^{2}>0$ ) in more detail. In this case, we can associate a linear $2\times 2$ differential system to (1) which contains only two parameters in addition to the eigenvalue parameter, so that the partial differential equation for the eigenvalues can be further simplified.

2. A linear Hamiltonian system associated to the CSWE

In this section we study the linear $2\times 2$ Hamiltonian system (4) on the interval $(0,1)$ . It depends on three parameters $u_{1},u_{2},u_{3}\in\mathbb{R}$ , where $\alpha=\frac{1}{2}(\mu+1)\geq\frac{1}{2}$ with some constant $\mu\geq 0$ , and the prime denotes the derivative with respect to $z$ . If we define

J:=\begin{pmatrix}[r]0&-1\\[4.30554pt] 1&0\end{pmatrix},\quad H(z,\mathbf{u}):=\begin{pmatrix}2&-\frac{\alpha}{z}+\frac{\alpha}{1-z}+2u_{2}\\[4.30554pt] -\frac{\alpha}{z}+\frac{\alpha}{1-z}+2u_{2}&\frac{u_{3}}{z}+2(u_{1}+u_{2}^{2})\end{pmatrix},\quad W(z):=\begin{pmatrix}\frac{1}{1-z}&0\\[4.30554pt] 0&\frac{1}{z}\end{pmatrix}

and

\tau y:=W(z)^{-1}\left(Jy^{\prime}(z)-H(z,\mathbf{u})y(z)\right)

(6)

where $\mathbf{u}:=(u_{1},u_{2},u_{3})\in\mathbb{R}^{3}$ combines the three parameters, then (4) can be written in the form $\tau y=\Lambda y$ . Moreover, as $W(z)^{\ast}=W(z)>0$ and $H(z,\mathbf{u})^{\ast}=H(z,\mathbf{u})$ for all $z\in(0,1)$ , $\tau$ defines for fixed parameters $\mathbf{u}\in\mathbb{R}^{3}$ a formally self-adjoint differential expression in the Hilbert space

\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2}):=\Big{\{}f:(0,1)\longrightarrow\mathbb{C}^{2}\ \big{|}\ \int_{0}^{1}f(z)^{\ast}W(z)f(z)\,\mathrm{d}z<\infty\Big{\}}

of square-integrable vector functions to the weight function $W(z)$ with scalar product

\langle f,g\rangle_{W}:=\int_{0}^{1}f(z)^{\ast}W(z)\,g(z)\,\mathrm{d}z

The maximal operator $Ty:=\tau y$ generated by $\tau$ has the domain (cf. [5, Section 3])

D(T):=\left\{y\in\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})\ :\ y\in\mathrm{AC_{loc}}((0,1),\mathbb{C}^{2})\mbox{ and }\tau y\in\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})\right\}

(7)

whereas the domain of the minimal operator $T_{0}y:=\tau y$ associated to $\tau$ is given by

D(T_{0}):=\big{\{}y\in D(T)\,:\,y\mbox{ has compact support in }(0,1)\big{\}}

In order to establish the self-adjointness of $T$ , we first prove that $\tau$ is in the limit point case at $z=0$ and $z=1$ . To this end, we consider the differential equation $\tau y=\Lambda y$ in the complex plane. It is equivalent to the $2\times 2$ system

y^{\prime}(z)=\left(\frac{1}{z}\begin{pmatrix}-\alpha&u_{3}+\Lambda\\[4.30554pt] 0&\alpha\end{pmatrix}+\frac{1}{1-z}\begin{pmatrix}\alpha&0\\[4.30554pt] -\Lambda&-\alpha\end{pmatrix}+\begin{pmatrix}2u_{2}&2(u_{1}+u_{2}^{2})\\[4.30554pt] -2&-2u_{2}\end{pmatrix}\right)y(z)

(8)

with regular singular points at $z=0$ and $z=1$ . If $\mathfrak{D}_{0}:=\{z\in\mathbb{C}:|z|<1\}$ denotes the unit disk centered at $0$ , then [3, Lemma 6] implies that for fixed $(\Lambda,\mathbf{u})\in\mathbb{C}\times\mathbb{R}^{3}$ there exists a fundamental matrix of the form $Y_{0}(z)=F_{0}(z,\Lambda,\mathbf{u})z^{A}z^{Q_{0}(\Lambda,\mathbf{u})}$ , $z\in\mathfrak{D}_{0}$ , where $F_{0}:\mathfrak{D}_{0}\times\mathbb{C}\times\mathbb{R}^{3}\longrightarrow\mathrm{M}_{2}(\mathbb{C})$ is an analytic matrix function,

F_{0}(0,\Lambda,\mathbf{u})=\begin{pmatrix}1&\frac{u_{3}+\Lambda}{2\alpha}\\[4.30554pt] 0&1\end{pmatrix},\quad A:=\begin{pmatrix}-\alpha&0\\[4.30554pt] 0&\alpha\end{pmatrix},\quad Q_{0}(\Lambda,\mathbf{u})=\begin{pmatrix}0&0\\[4.30554pt] q_{0}(\Lambda,\mathbf{u})&0\end{pmatrix}

and $q_{0}:\mathbb{C}\times\mathbb{R}^{3}\longrightarrow\mathbb{C}$ is an analytic function with $q_{0}\equiv 0$ if $2\alpha=\mu+1$ is not an integer. Furthermore, the transformation $\tilde{y}(z)=y(1-z)$ yields

\tilde{y}^{\prime}(z)=\left(\frac{1}{z}\begin{pmatrix}-\alpha&0\\[4.30554pt] \Lambda&\alpha\end{pmatrix}+\frac{1}{1-z}\begin{pmatrix}\alpha&-u_{3}-\Lambda\\[4.30554pt] 0&-\alpha\end{pmatrix}+\begin{pmatrix}-2u_{2}&-2(u_{1}+u_{2}^{2})\\[4.30554pt] 2&2u_{2}\end{pmatrix}\right)\tilde{y}(z)

and hence (8) also possesses a fundamental matrix of the form $Y_{1}(z)=F_{1}(z,\Lambda,\mathbf{u})(1-z)^{A}(1-z)^{Q_{1}(\Lambda,\mathbf{u})}$ , $z\in\mathfrak{D}_{1}:=\{z\in\mathbb{C}:|z-1|<1\}$ , where $F_{1}:\mathfrak{D}_{1}\times\mathbb{C}\times\mathbb{R}^{3}\longrightarrow\mathrm{M}_{2}(\mathbb{C})$ is an analytic matrix function,

F_{1}(1,\Lambda,\mathbf{u})=\begin{pmatrix}1&0\\[4.30554pt] -\frac{\Lambda}{2\alpha}&1\end{pmatrix},\quad Q_{1}(\Lambda,\mathbf{u})=\begin{pmatrix}0&0\\[4.30554pt] q_{1}(\Lambda,\mathbf{u})&0\end{pmatrix}

and $q_{1}:\mathbb{C}\times\mathbb{R}^{3}\longrightarrow\mathbb{C}$ being an analytic function satisfying $q_{1}\equiv 0$ if $2\alpha=\mu+1$ is not an integer. For fixed parameters $(\Lambda,\mathbf{u})\in\mathbb{C}\times\mathbb{R}^{3}$ , the solution

y_{0}(z):=Y_{0}(z)e_{1}=F_{0}(z)\begin{pmatrix}z^{-\alpha}&0\\[4.30554pt] 0&z^{\alpha}\end{pmatrix}\begin{pmatrix}1&0\\[4.30554pt] q_{0}\log z&1\end{pmatrix}\begin{pmatrix}1\\[4.30554pt] 0\end{pmatrix}

of (8) asymptotically behaves like $y_{0}(z)\sim z^{-\alpha}e_{1}$ as $z\to 0$ (here and in the following, $\{e_{1},e_{2}\}$ denotes the standard basis of $\mathbb{C}^{2}$ ). Since $2\alpha\geq 1$ and $y_{0}(z)^{\ast}W(z)y_{0}(z)\sim z^{-2\alpha}$ as $z\to 0$ , the function $y_{0}(z)^{\ast}W(z)y_{0}(z)$ is not integrable at $z=0$ , and hence $y_{0}$ does not lie left in $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ . Moreover, for the solution

y_{1}(z):=Y_{1}(z)e_{1}=F_{1}(z)\begin{pmatrix}(1-z)^{-\alpha}&0\\[4.30554pt] 0&(1-z)^{\alpha}\end{pmatrix}\begin{pmatrix}1&0\\[4.30554pt] q_{1}\log(1-z)&1\end{pmatrix}\begin{pmatrix}1\\[4.30554pt] 0\end{pmatrix}

we get $y_{1}(z)^{\ast}W(z)y_{1}(z)\sim(1-z)^{-2\alpha-1}$ as $z\to 1$ . Therefore, $y_{1}(z)^{\ast}W(z)y_{1}(z)$ is not integrable at $z=1$ , and thus $y_{1}$ does not lie right in $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ . According to Weyl’s alternative (see e. g. [5, Theorem 5.6]), the differential expression $\tau$ is in the limit point case at both endpoints. Hence, the closure of its associated minimal operator $\overline{T_{0}}$ is a self-adjoint extension of $T_{0}$ (cf. [6, Theorem 5.4]), and because of $\overline{T_{0}}=\overline{T_{0}}^{\ast}=T$ (see [5, Theorem 3.9]), this extension coincides with the maximal operator $T$ , so that $T$ is actually a self-adjoint differential operator in the Hilbert space $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ .

In the following, we will study the eigenvalues of the self-adjoint operator $T=T(\mathbf{u})$ associated to $\tau$ and their dependence on the parameters $\mathbf{u}\in\mathbb{R}^{3}$ . For this purpose, we introduce the functions $\eta_{0}(z,\Lambda,\mathbf{u}):=Y_{0}(z,\Lambda,\mathbf{u})e_{2}$ and $\eta_{1}(z,\Lambda,\mathbf{u}):=Y_{1}(z,\Lambda,\mathbf{u})e_{2}$ , which asymptotically behave like

\eta_{0}(z,\Lambda,\mathbf{u})\sim z^{\alpha}\begin{pmatrix}\frac{u_{3}+\Lambda}{2\alpha}\\[4.30554pt] 1\end{pmatrix}\quad\mbox{as $z\to 0$},\quad\eta_{1}(z,\Lambda,\mathbf{u})\sim(1-z)^{\alpha}\begin{pmatrix}0\\[4.30554pt] 1\end{pmatrix}\quad\mbox{as $z\to 1$}

(9)

Note that $y_{0}$ , $\eta_{0}$ as well as $y_{1}$ , $\eta_{1}$ form a fundamental system of (4). Since the solutions $y_{0}(z)=Y_{0}(z)e_{1}$ and $y_{1}(z)=Y_{1}(z)e_{1}$ are not square-integrable with respect to $W(z)$ , $\Lambda\in\mathbb{R}$ is an eigenvalue of $T(\mathbf{u})$ for fixed $\mathbf{u}\in\mathbb{R}^{3}$ if and only if there exists a nontrivial solution $y$ of (4) which is a constant multiple of both $\eta_{0}$ and $\eta_{1}$ . Additionally, by virtue of $\alpha\geq\frac{1}{2}$ and the asymptotic behavior of the fundamental solutions at the boundary points, a nontrivial solution $y$ of (4) is an eigenfunction if and only if $(z(1-z))^{-1/2}y(z)$ is bounded on $(0,1)$ . Furthermore, we can write $\eta_{0}$ as a linear combination

\eta_{0}(z,\Lambda,\mathbf{u})=\Theta(\Lambda,\mathbf{u})y_{1}(z,\Lambda,\mathbf{u})+\Omega(\Lambda,\mathbf{u})\eta_{1}(z,\Lambda,\mathbf{u})

(10)

with $\Theta$ , $\Omega$ called connection coefficients. Since for any fixed $z_{0}\in(0,1)$

\begin{pmatrix}\Theta(\Lambda,\mathbf{u})\\[2.15277pt] \Omega(\Lambda,\mathbf{u})\end{pmatrix}=Y_{1}(z_{0},\Lambda,\mathbf{u})^{-1}Y_{0}(z_{0},\Lambda,\mathbf{u})e_{2}

where the fundamental matrices depend analytically on the parameters, also $\Theta=\Theta(\Lambda,\mathbf{u})$ and $\Omega=\Omega(\Lambda,\mathbf{u})$ depend analytically on $(\Lambda,\mathbf{u})\in\mathbb{C}\times\mathbb{R}^{3}$ . Due to our considerations above, $\Lambda\in\mathbb{R}$ is an eigenvalue of $T(\mathbf{u})$ if and only if $\Theta(\Lambda,\mathbf{u})=0$ . As $\Theta(\Lambda,\mathbf{u})$ is an entire function for fixed $\mathbf{u}\in\mathbb{R}^{3}$ , the zeros of $\Theta(\Lambda,\mathbf{u})$ are isolated, and thus also the eigenvalues of $T(\mathbf{u})$ are isolated with multiplicity $1$ .

In a first step, we will figure out how the eigenvalues of (2) and (4) are related.

Lemma 1.

$\Lambda=0$ is an eigenvalue of the self-adjoint operator $T(u_{1},u_{2},u_{3})$ associated to the linear Hamiltonian system (4) if and only if $\lambda=u_{3}+\mu(\mu+1)$ is an eigenvalue of the Coulomb spheroidal wave equation (2) for the parameters $\gamma^{2}=u_{1}$ and $\beta=-u_{3}-2(\mu+1)u_{2}$ .

Proof.

The system (8) with $\Lambda=0$ is equivalent to

y^{\prime}(z)=\left(\frac{1}{z}\begin{pmatrix}-\alpha&u_{3}\\[4.30554pt] 0&\alpha\end{pmatrix}+\frac{1}{1-z}\begin{pmatrix}\alpha&0\\[4.30554pt] 0&-\alpha\end{pmatrix}+\begin{pmatrix}2u_{2}&2(u_{1}+u_{2}^{2})\\[4.30554pt] -2&-2u_{2}\end{pmatrix}\right)y(z)

If we introduce the vector function

\hat{y}(x)=\begin{pmatrix}v(x)\\[2.15277pt] w(x)\end{pmatrix}:=2(1-x^{2})^{-\frac{1}{2}}y\left(\tfrac{x+1}{2}\right),\quad x\in(-1,1)

(11)

then this differential system takes the form

\hat{y}^{\prime}(x)=\left(\frac{1}{1+x}\begin{pmatrix}-\alpha-\frac{1}{2}&u_{3}\\[4.30554pt] 0&\alpha-\frac{1}{2}\end{pmatrix}+\frac{1}{1-x}\begin{pmatrix}\alpha+\frac{1}{2}&0\\[4.30554pt] 0&-\alpha+\frac{1}{2}\end{pmatrix}+\begin{pmatrix}u_{2}&u_{1}+u_{2}^{2}\\[4.30554pt] -1&-u_{2}\end{pmatrix}\right)\hat{y}(x)

(12)

Since $\hat{y}(x)=(z(1-z))^{-1/2}y(z)$ with $z=\frac{1}{2}(x+1)\in(0,1)$ , we obtain that $\Lambda=0$ is an eigenvalue of $T(\mathbf{u})$ if and only if (12) has a nontrivial solution $\hat{y}(x)$ which is bounded on $(-1,1)$ . Using $\alpha=\frac{1}{2}(\mu+1)$ , we can write (12) in terms of its components $v(x)$ and $w(x)$ :

	$\displaystyle v^{\prime}(x)$	$\displaystyle=\left(u_{2}+\frac{(\mu+2)x}{1-x^{2}}\right)v(x)+\left(u_{1}+u_{2}^{2}+\frac{u_{3}}{1+x}\right)w(x)$
	$\displaystyle w^{\prime}(x)$	$\displaystyle=-v(x)-\left(u_{2}+\frac{\mu x}{1-x^{2}}\right)\,w(x)$

Solving the second equation for $v(x)$ yields

v(x)=-w^{\prime}(x)-\left(u_{2}+\frac{\mu x}{1-x^{2}}\right)\,w(x)

(13)

and inserting this expression into the first equation gives the differential equation

\frac{\mathrm{d}}{\mathrm{d}x}\left((1-x^{2})\frac{\mathrm{d}}{\mathrm{d}x}w(x)\right)+\left(u_{3}+\mu(\mu+1)-(u_{3}+2(\mu+1)u_{2})x+u_{1}(1-x^{2})-\frac{\mu^{2}}{1-x^{2}}\right)w(x)=0

for $w(x)$ which is a CSWE (2) with parameters (3). Conversely, any bounded solution $w(x)$ of this ODE is a Coulomb spheroidal function having the form $w(x)=(1-x^{2})^{\mu/2}\varphi(x)$ with some entire function $\varphi(x)$ , cf. [7, Section 3.2, Satz 1]. If we define $v(x)$ by (13), then $v(x)=-(1-x^{2})^{\mu/2}\varphi^{\prime}(x)-u_{2}(1-x^{2})^{\mu/2}\varphi(x)$ is also bounded on $(-1,1)$ , and hence the function $\hat{y}$ defined by (11) is a bounded solution of (12) which implies that $\Lambda=0$ is an eigenvalue of (4). ∎

In the next step we examine which way the eigenvalues $\Lambda$ of (4) depend on the parameters $(u_{1},u_{2},u_{3})$ .

Theorem 2.

Let $T=T(\mathbf{u})$ be the self-adjoint realization of the differential expression (6) with domain $D(T)$ given by (7), and suppose that $\Lambda_{0}$ is an eigenvalue of $T=T(\mathbf{u}_{0})$ for some given parameter vector $\mathbf{u}_{0}\in\mathbb{R}^{3}$ . Then there exists a domain $U\subset\mathbb{R}^{3}$ with $\mathbf{u}_{0}\in U$ and eigenvalues $\Lambda(\mathbf{u})$ of $T(\mathbf{u})$ with $\Lambda(\mathbf{u}_{0})=\Lambda_{0}$ , such that $\Lambda=\Lambda(\mathbf{u})$ depends analytically on $\mathbf{u}=(u_{1},u_{2},u_{3})\in U$ . In addition, $\Lambda$ satisfies the first-order partial differential equation (5) on $U$ .

Proof.

We first prove the analytic dependence of the eigenvalues on the parameters $\mathbf{u}=(u_{1},u_{2},u_{3})$ in a neighborhood of $\mathbf{u}_{0}=:(a,b,c)$ . For this, let

B(z,\mathbf{u}):=H(z,\mathbf{u})-H(z,\mathbf{u}_{0})=\begin{pmatrix}0&2(u_{2}-b)\\[4.30554pt] 2(u_{2}-b)&2(u_{1}+u_{2}^{2}-a-b^{2})+\frac{u_{3}-c}{z}\end{pmatrix}

If $A(\mathbf{u})$ denotes the symmetric multiplication operator $A(\mathbf{u})y:=-W(z)^{-1}B(z,\mathbf{u})y$ in $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ , then $T(\mathbf{u})=T(\mathbf{u}_{0})+A(\mathbf{u})$ . Since

W(z)^{-1/2}B(z,\mathbf{u})W(z)^{-1/2}=\begin{pmatrix}0&2(u_{2}-b)\sqrt{z(1-z)}\\[4.30554pt] 2(u_{2}-b)\sqrt{z(1-z)}&2(u_{1}+u_{2}^{2}-a-b^{2})z+u_{3}-c\end{pmatrix}

we have $|W(z)^{-1/2}B(z,\mathbf{u})W(z)^{-1/2}|\leq C$ for all $z\in(0,1)$ with some constant $C=C(\mathbf{u})\geq 0$ , and

	$\displaystyle\\|Ay\\|_{W}^{2}$	$\displaystyle=\int_{0}^{1}\left(W(z)^{-1}B(z,\mathbf{u})\,y(z)\right)^{\ast}W(z)\left(W(z)^{-1}B(z,\mathbf{u})\,y(z)\right)\,\mathrm{d}z$
		$\displaystyle=\int_{0}^{1}\|W(z)^{-1/2}B(z,\mathbf{u})W(z)^{-1/2}W(z)^{1/2}y(z)\|^{2}\,\mathrm{d}z\leq\int_{0}^{1}C^{2}\|W(z)^{1/2}y(z)\|^{2}\,\mathrm{d}z=C^{2}\\|y\\|_{W}^{2}$

implies that $A=A(\mathbf{u})$ is even a bounded operator on $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ for all $\mathbf{u}\in\mathbb{C}^{3}$ . From [8, Chap. V, § 4, Theorem 4.3 and Chap. VII, § 2, Theorem 2.6] it follows that $T(\mathbf{u})$ defines a holomorphic family of self-adjoint operators, where its domain $D(T)$ is independent of $\mathbf{u}\in\mathbb{C}^{3}$ . As shown in [8, Chap. VII, § 2, Sec. 4 and § 3, Sec. 1 – 2] there exists a domain $U\subset\mathbb{R}^{3}$ with $\mathbf{u}_{0}\in U$ as well as an analytical function $\Lambda=\Lambda(u_{1},u_{2},u_{3})$ with $\Lambda(\mathbf{u}_{0})=\Lambda_{0}$ , such that $\Lambda(\mathbf{u})$ is a simple eigenvalue of $T(\mathbf{u})$ for all $\mathbf{u}\in U$ .

Now, if $\Lambda(\mathbf{u})$ is an eigenvalue of $T(\mathbf{u})$ , then each corresponding eigenfunction is a constant multiple of the fundamental solutions $\eta_{0}$ and $\eta_{1}$ , respectively. Hence, $\eta_{0}(z)=\Omega(\mathbf{u})\cdot\eta_{1}(z)$ for all $z\in(0,1)$ with some factor $\Omega(\mathbf{u})\neq 0$ which depends analytically on $\mathbf{u}\in U$ , and therefore $\eta_{0}(z)=z^{\alpha}(1-z)^{\alpha}h(z,\mathbf{u})$ , where $h:\mathfrak{D}\times U\longrightarrow\mathbb{C}^{2}$ is an analytic function on some domain $[0,1]\subset\mathfrak{D}\subset\mathbb{C}$ which satisfies

h(0,\mathbf{u})=\begin{pmatrix}\frac{u_{3}+\Lambda(\mathbf{u})}{2\alpha}\\[4.30554pt] 1\end{pmatrix}\quad\mbox{and}\quad h(1,\mathbf{u})=\begin{pmatrix}0\\[4.30554pt] \Omega(\mathbf{u})\end{pmatrix}

due to (9). Obviously $\eta_{0}\in\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ , and

\nu(\mathbf{u}):=\langle\eta_{0},\eta_{0}\rangle_{W}=\int_{0}^{1}h(z,\mathbf{u})^{\ast}\begin{pmatrix}z^{\mu+1}(1-z)^{\mu}&0\\[4.30554pt] 0&z^{\mu}(1-z)^{\mu+1}\end{pmatrix}h(z,\mathbf{u})\,\mathrm{d}z

defines a smooth function $\nu:U\longrightarrow(0,\infty)$ . Setting $\psi(z,\mathbf{u}):=\nu(\mathbf{u})^{-1/2}h(z,\mathbf{u})$ , then

y(z,\mathbf{u}):=\nu(\mathbf{u})^{-1/2}\eta_{0}(z,\Lambda(\mathbf{u}),\mathbf{u})=z^{\alpha}(1-z)^{\alpha}\psi(z,\mathbf{u})

is an eigenfunction corresponding to $\Lambda(\mathbf{u})$ for which $\langle y,y\rangle_{W}=1$ holds for all $\mathbf{u}\in U$ , and $y=y(z,\mathbf{u})$ depends infinitely differentiable on $(z,\mathbf{u})\in(0,1)\times U$ . Since the matrix functions $W:(0,1)\longrightarrow\mathrm{M}_{2}(\mathbb{C})$ and $H:(0,1)\times\mathbb{R}^{3}\longrightarrow\mathrm{M}_{2}(\mathbb{C})$ are also infinitely differentiable with respect to all variables, the assumptions (a), (b), (c) in [4] are fulfilled. In addition, if we introduce the matrix function

G(z,\Lambda,\mathbf{u}):=\begin{pmatrix}2z&2zu_{2}+\Lambda u_{2}-\mu-\frac{1}{2}\\[4.30554pt] 2zu_{2}+\Lambda u_{2}-\mu-\frac{1}{2}&2(z-1)(u_{1}+u_{2}^{2})\end{pmatrix}

then a direct calculation gives

\frac{\partial G}{\partial z}+(\Lambda W+H)JG-GJ(\Lambda W+H)=\sum_{k=1}^{3}f_{k}\frac{\partial H}{\partial u_{k}}+gW

(14)

where the functions $f_{k}=f_{k}(\Lambda,\mathbf{u})$ and $g=g(\Lambda,\mathbf{u})$ on the right hand side are defined by

	$\displaystyle f_{1}$	$\displaystyle:=2u_{1}$
	$\displaystyle f_{2}$	$\displaystyle:=(\Lambda+2)(u_{1}+u_{2}^{2})+\Lambda+u_{2}+u_{3}$
	$\displaystyle f_{3}$	$\displaystyle:=(2\Lambda+u_{3}+2)(2\Lambda u_{2}-2\mu-1)-2(\mu+1)(u_{1}+u_{2}^{2})+2\mu$
	$\displaystyle g$	$\displaystyle:=(1+2\mu-2\Lambda u_{2})(\Lambda+2)-2\mu$

Here, $G:(0,1)\times\mathbb{C}\times U\longrightarrow\mathrm{M}_{2}(\mathbb{C})$ and $f_{k},g:\mathbb{C}\times U\longrightarrow\mathbb{C}$ are differentiable functions, and therefore also assumption (d) in [4] is satisfied. In addition to the \qqdeformation equation (14), the function

\displaystyle\sum_{k=1}^{3}f_{k}(\Lambda(\mathbf{u}),\mathbf{u})\frac{\partial y}{\partial u_{k}}+JG(z,\Lambda(\mathbf{u}),\mathbf{u})y=z^{\alpha}(1-z)^{\alpha}\left(\sum_{k=1}^{3}f_{k}(\Lambda(\mathbf{u}),\mathbf{u})\frac{\partial\psi}{\partial u_{k}}+JG(z,\Lambda(\mathbf{u}),\mathbf{u})\psi\right)

lies in $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ . Finally, from [4, Theorem 2.2] it follows that $\Lambda=\Lambda(u_{1},u_{2},u_{3})$ solves

\sum_{k=1}^{3}f_{k}(\Lambda,u_{1},u_{2},u_{3})\frac{\partial\Lambda}{\partial u_{k}}=g(\Lambda,u_{1},u_{2},u_{3})

on $U$ , and this is exactly the first-order quasilinear PDE (5). ∎

The eigenvalues of (4) coincide with the zeros of the function $\Theta(\Lambda,\mathbf{u})$ , and this is one of the connection coefficients appearing in (10). The question arises how this function can be calculated numerically for given parameter values $\mathbf{u}\in\mathbb{R}^{3}$ and $\Lambda\in\mathbb{C}$ . In [9], the Coulomb spheroidal wave equation has already been associated to a $2\times 2$ differential system, which is, however, different from the linear Hamiltonian system studied in the present paper. Nevertheless, the method suggested in [9] for calculating the connection coefficient $\Theta$ can also be applied to the system (4). To this end, we apply the transformation $\hat{y}(z)=(1-z)^{\alpha-1}y(z)$ , which turns (4) resp. (8) into the differential system

\hat{y}^{\prime}(z)=\left(\frac{1}{z}\,A_{0}+\frac{1}{z-1}\,A_{1}+C\right)\hat{y}(z)

(15)

where

A_{0}=\begin{pmatrix}-\alpha&u_{3}+\Lambda\\[4.30554pt] 0&\alpha\end{pmatrix},\quad A_{1}=\begin{pmatrix}-1&0\\[4.30554pt] \Lambda&2\alpha-1\end{pmatrix},\quad C=\begin{pmatrix}2u_{2}&2(u_{1}+u_{2}^{2})\\[4.30554pt] -2&-2u_{2}\end{pmatrix}

For fixed $(\lambda,\mathbf{u})\in\mathbb{C}\times\mathbb{R}^{3}$ it has the Floquet solution

\hat{\eta}_{0}(z):=(1-z)^{\alpha-1}\eta_{0}(z)=z^{\alpha}h(z),\quad h(z)=\sum_{k=0}^{\infty}z^{k}d_{k},\quad d_{0}:=\begin{pmatrix}\frac{u_{3}+\Lambda}{2\alpha}\\[4.30554pt] 1\end{pmatrix}

where $h$ is holomorphic in $\mathfrak{D}_{0}$ , and $d_{0}$ is an eigenvector of $A_{0}$ corresponding to the eigenvalue $\alpha$ . Finally, once we have determined the series coefficients $d_{k}$ using the recurrence relation

b_{k}:=\begin{pmatrix}\frac{2(ku_{2}-u_{3})}{k(k+\mu+1)}&\frac{2k(u_{1}+u_{2}^{2})-2u_{2}u_{3}}{k(k+\mu+1)}\\[4.30554pt] -\frac{2}{k}&-\frac{2u_{2}}{k}\end{pmatrix}b_{k-1}-\begin{pmatrix}0&\frac{(\mu+1)u_{3}}{k(k+\mu+1)}\\[4.30554pt] 0&\frac{\mu+1}{k}\end{pmatrix}d_{k-1},\quad d_{k}:=d_{k-1}+b_{k}

(16)

for $k=1,2,3,\ldots$ starting with $b_{0}:=d_{0}$ , then we can compute $\Theta$ with arbitrary precision by means of

Lemma 3.

Let $\mu\geq 0$ and $\mathbf{u}=(u_{1},u_{2},u_{3})\in\mathbb{R}^{3}$ be fixed. Then for each $\Lambda\in\mathbb{C}$ we have

\Theta(\Lambda,\mathbf{u}):=\lim_{k\to\infty}e_{1}^{\operatorname{T}}d_{k}(\Lambda,\mathbf{u})

More precisely, we obtain $e_{1}^{\operatorname{T}}d_{k}=\Theta+\text{O}(k^{\delta-\mu-2})$ as $k\to\infty$ for arbitrary small $\delta>0$ .

This assertion, as well as formula (16), is obtained by a similar reasoning as given in [9, Section 3], and therefore the proof is omitted here. If we set $\Lambda=0$ and calculate the zeros $\mathbf{u}=(u_{1},u_{2},u_{3})$ of the function $\Theta(0,\mathbf{u})$ , then we get the eigenvalues $\lambda=u_{3}+\mu(\mu+1)$ of the Coulomb spheroidal wave equation (2) for the parameters $\gamma^{2}=u_{1}$ and $\beta=-u_{3}-2(\mu+1)u_{2}$ according to Lemma 1. This approach, like the method presented in [9], provides a simple and efficient way to obtain the eigenvalues of a CSWE. However, with increasing parameter values, numerical effects such as rounding errors or digit cancellation may occur, which require some computational effort to overcome. In this context the partial differential equation (5) may be helpful. It not only provides us with new insights into the analytical dependence of the eigenvalues $\Lambda=\Lambda(\mathbf{u})$ on the parameters, but can also be used to determine the Coulomb spheroidal eigenvalues numerically. For this purpose, we need to examine the PDE and its characteristics in more detail.

3. Analysis of the PDE for the eigenvalues

The PDE (5) looks very complicated at first sight, but it has some remarkable properties that simplify the calculation of the solution by the method of characteristics. The characteristic curves related to the quasilinear PDE (5) satisfy the (nonlinear) first-order differential system

\begin{split}\frac{\mathrm{d}u_{1}}{\mathrm{d}t}&=2u_{1}(t)\\ \frac{\mathrm{d}u_{2}}{\mathrm{d}t}&=(\Lambda(t)+2)(u_{1}(t)+u_{2}(t)^{2})+\Lambda(t)+u_{2}(t)+u_{3}(t)\\ \frac{\mathrm{d}u_{3}}{\mathrm{d}t}&=(2\Lambda(t)+u_{3}(t)+2)(2\Lambda(t)u_{2}(t)-2\mu-1)-2(\mu+1)(u_{1}(t)+u_{2}(t)^{2})+2\mu\\ \frac{\mathrm{d}\Lambda}{\mathrm{d}t}&=(1+2\mu-2\Lambda(t)u_{2}(t))(\Lambda(t)+2)-2\mu\end{split}

(17)

We can choose the parameter for the characteristic curve in such a way that the initial conditions are given at $t=0$ . If we define the scalar function

\Psi(t):=\Lambda(t)(u_{1}(t)+u_{2}(t)^{2}-1)-2(\mu+1)u_{2}(t)-u_{3}(t)

then (17) and a straightforward calculation yields $\frac{\mathrm{d}}{\mathrm{d}t}\Psi=\Psi$ . In particular, if $\Psi(0)=0$ , then $\Psi\equiv 0$ , and therefore $u_{3}(t)=\Lambda(t)(u_{1}(t)+u_{2}(t)^{2}-1)-2(\mu+1)u_{2}(t)$ holds for all $t$ along such a characteristic curve. If, in addition, $\Lambda(t_{0})=0$ for some $t_{0}\in\mathbb{R}$ , then we get $u_{3}(t_{0})=-2(\mu+1)u_{2}(t_{0})$ . According to Lemma 1, $\Lambda(t_{0})=0$ also implies that $\lambda=u_{3}(t_{0})+\mu(\mu+1)$ is an eigenvalue of the Coulomb spheroidal wave equation (2) for the parameters $\gamma^{2}=u_{1}(t_{0})$ and $\beta=-u_{3}(t_{0})-2(\mu+1)u_{2}(t_{0})=0$ , where $\beta=0$ in turn reduces the CSWE to the ordinary (angular) spheroidal wave equation. All in all, the above considerations provide an alternative method for the numerical computation of spheroidal eigenvalues:

Theorem 4.

Suppose that $u_{1}(t)$ , $u_{2}(t)$ and $\Lambda(t)$ are solutions of the differential system

\begin{split}\frac{\mathrm{d}u_{1}}{\mathrm{d}t}&=2u_{1}(t)\\ \frac{\mathrm{d}u_{2}}{\mathrm{d}t}&=2(\Lambda(t)+1)(u_{1}(t)+u_{2}(t)^{2})-(2\mu+1)u_{2}(t)\\ \frac{\mathrm{d}\Lambda}{\mathrm{d}t}&=(1+2\mu-2\Lambda(t)u_{2}(t))(\Lambda(t)+2)-2\mu\end{split}

(18)

on an interval $0\in I\subset\mathbb{R}$ satisfying $u_{1}(0)=a$ , $u_{2}(0)=b$ , $\Lambda(0)=\Lambda_{0}$ with some given initial values $a,b\in\mathbb{R}$ such that $\Theta(\Lambda_{0},a,b,(a+b^{2}-1)\Lambda_{0}-2(\mu+1)b)=0$ . If $\Lambda(t_{0})=0$ holds for some $t_{0}\in I$ , then $\lambda=(\mu+1)(\mu-2u_{2}(t_{0}))$ is an eigenvalue of the angular spheroidal wave equation (1) for the parameter $\gamma^{2}=u_{1}(t_{0})$ .

Proof.

If we set $c:=(a+b^{2}-1)\Lambda_{0}-2(\mu+1)b$ , then $\Theta(\Lambda_{0},a,b,c)=0$ , and hence $\Lambda_{0}$ is an isolated eigenvalue of the operator $T(\mathbf{u}_{0})$ for $\mathbf{u}_{0}:=(a,b,c)$ . According to Theorem 2, we can find an analytic continuation $\Lambda_{0}:U\longrightarrow\mathbb{R}$ on some domain $U\subset\mathbb{R}^{3}$ with $\mathbf{u}_{0}\in U$ , such that $\Lambda_{0}(\mathbf{u})$ is an eigenvalue of $T(\mathbf{u})$ satisfying, in addition, the partial differential equation (5). Defining $u_{3}(t):=\Lambda(t)(u_{1}(t)+u_{2}(t)^{2}-1)-2(\mu+1)u_{2}(t)$ , then we can write the second equation in (18) as

\frac{\mathrm{d}u_{2}}{\mathrm{d}t}=(\Lambda(t)+2)(u_{1}(t)+u_{2}(t)^{2})+u_{2}(t)+u_{3}(t)+\Lambda(t)

Moreover, the differential equations (18) imply

	$\displaystyle\frac{\mathrm{d}u_{3}}{\mathrm{d}t}$	$\displaystyle=\frac{\mathrm{d}\Lambda}{\mathrm{d}t}(u_{1}(t)+u_{2}(t)^{2}-1)+\Lambda(t)\left(\frac{\mathrm{d}u_{1}}{\mathrm{d}t}+2u_{2}(t)\frac{\mathrm{d}u_{2}}{\mathrm{d}t}\right)-2(\mu+1)\frac{\mathrm{d}u_{2}}{\mathrm{d}t}$
		$\displaystyle=(2\Lambda(t)+u_{3}(t)+2)(2\Lambda(t)u_{2}(t)-2\mu-1)-2(\mu+1)(u_{1}(t)+u_{2}(t)^{2})+2\mu$

Hence, the functions $u_{1}(t)$ , $u_{2}(t)$ , $u_{3}(t)$ , $\Lambda(t)$ satisfy the characteristic equations (17), and from the initial values $u_{1}(0)=a$ , $u_{2}(0)=b$ , $u_{3}(0)=c$ , $\Lambda(0)=\Lambda_{0}=\Lambda_{0}(a,b,c)$ it follows that this characteristic curve lies completely on the integral surface $\Lambda_{0}(\mathbf{u})$ . In particular, $\Lambda(t)=\Lambda_{0}(u_{1}(t),u_{2}(t),u_{3}(t))$ , and therefore $\Lambda(t)$ is an eigenvalue of (4) for the parameter values $u_{1}(t)$ , $u_{2}(t)$ , $u_{3}(t)$ . Finally, if $\Lambda(t_{0})=0$ , then $u_{3}(t_{0})=-2(\mu+1)u_{2}(t_{0})$ implies $\beta=-u_{3}(t_{0})-2(\mu+1)u_{2}(t_{0})=0$ , and $\lambda=(\mu+1)(\mu-2u_{2}(t_{0}))$ is an eigenvalue of (1) for the spheroidal parameter $\gamma^{2}=u_{1}(t_{0})$ . ∎

We may, for example, choose $b=0$ in Theorem 4, and then we need to find a zero $\Lambda_{0}$ of the function $\tilde{\Theta}(\Lambda):=\Theta(\Lambda,a,0,(a-1)\Lambda)$ for some given initial value $a\in\mathbb{R}$ . Then we have to solve the differential system (18) with $u_{1}(0)=a$ , $u_{2}(0)=0$ , $\Lambda(0)=\Lambda_{0}$ on some interval $I$ with $0\in I$ by an appropriate numerical method, and we follow the path of this characteristic curve until we get $\Lambda(t_{0})=0$ for some $t_{0}\in I$ . Finally, we arrive at an eigenvalue $\lambda=(\mu+1)(\mu-2u_{2}(t_{0}))$ of the angular spheroidal wave equation (1) for $\gamma^{2}=u_{1}(t_{0})$ . In the appendix, this method is illustrated by a concrete numerical example, namely for the values $\mu=1$ and $a=5$ .

The differential system (17) in Theorem 4 may also be regarded as the characteristic equation of a PDE for a function $\Lambda(u_{1},u_{2})$ with only two variables $u_{1}$ and $u_{2}$ . The question arises whether one can associate a more simple linear Hamiltonian system with the spheroidal wave equation ( $\beta=0$ ), which then yields a PDE for the eigenvalues that is less complicated than (5). In the remainder of this paper we address this problem, where we restrict our consideration to the case of the prolate spheroidal wave equation (2) with $\beta=0$ and $\gamma^{2}>0$ . For this purpose, we fix $u_{3}=\Lambda(u_{1}+u_{2}^{2}-1)-2(\mu+1)u_{2}$ in (4) resp. (8), and we assume $u_{1}=\gamma^{2}>0$ , in which case also $u_{1}+u_{2}^{2}>0$ holds.

First, let us introduce new parameters that make the coefficient matrix of the differential system (8) appear more symmetric, and these are

v_{1}:=\frac{u_{2}}{\sqrt{u_{1}+u_{2}^{2}}},\quad v_{2}:=\sqrt{u_{1}+u_{2}^{2}},\quad\zeta:=\frac{\Lambda(u_{1}+u_{2}^{2})-2\alpha u_{2}}{\sqrt{u_{1}+u_{2}^{2}}}

(19)

where $v_{1}\in(-1,1)$ and $v_{2}\in(0,\infty)$ . Note that conversely

u_{1}+u_{2}^{2}=v_{2}^{2},\quad u_{2}=v_{1}v_{2},\quad\Lambda=\frac{\zeta+2\alpha v_{1}}{v_{2}},\quad u_{3}+\Lambda=(\zeta-2\alpha v_{1})v_{2}

By means of the transformation

\hat{y}(z)=\begin{pmatrix}1&0\\[4.30554pt] 0&v_{2}\end{pmatrix}y(z)

(20)

the differential system (4) resp. (8) is equivalent to the $2\times 2$ system

\hat{y}^{\prime}(z)=\left(\frac{1}{z}\begin{pmatrix}-\alpha&\zeta-2\alpha v_{1}\\[4.30554pt] 0&\alpha\end{pmatrix}+\frac{1}{1-z}\begin{pmatrix}\alpha&0\\[4.30554pt] -\zeta-2\alpha v_{1}&-\alpha\end{pmatrix}+\begin{pmatrix}2v_{1}v_{2}&2v_{2}\\[4.30554pt] -2v_{2}&-2v_{1}v_{2}\end{pmatrix}\right)\hat{y}(z)

which can be written in the form $\hat{\tau}\hat{y}=\zeta\hat{y}$ , where

\hat{\tau}\hat{y}:=W(z)^{-1}\left(J\hat{y}^{\prime}(z)-\hat{H}(z)\hat{y}(z)\right),\quad\hat{H}(z):=\begin{pmatrix}2v_{2}+\frac{2\alpha v_{1}}{1-z}&-\frac{\alpha}{z}+\frac{\alpha}{1-z}+2v_{1}v_{2}\\[4.30554pt] -\frac{\alpha}{z}+\frac{\alpha}{1-z}+2v_{1}v_{2}&2v_{2}-\frac{2\alpha v_{1}}{z}\end{pmatrix}

By a similar reasoning as for (6), it follows that the differential expression $\hat{\tau}$ is in the limit point case at $z=0$ and $z=1$ . Hence, the closure of its associated minimal operator is the only self-adjoint realization of $\hat{\tau}$ in the Hilbert space $\mathrm{L}_{W}^{2}((0,1),\mathbb{C}^{2})$ . We denote this self-adjoint operator by $\hat{T}=\hat{T}(v_{1},v_{2})$ . By virtue of (20), we have $\langle y,y\rangle_{W}<\infty$ if and only if $\langle\hat{y},\hat{y}\rangle_{W}<\infty$ . Hence, $\Lambda$ is an eigenvalue of $T(\mathbf{u})$ if and only if $\zeta$ is an eigenvalue of $\hat{T}(v_{1},v_{2})$ .

Now, like in Theorem 2, suppose that $\zeta=\zeta(v_{1},v_{2})$ is an eigenvalue of $\hat{T}(v_{1},v_{2})$ , which depends analytically on the parameters $(v_{1},v_{2})\in V$ in some domain $V\subset(-1,1)\times(0,\infty)$ . If we introduce the matrix function

\hat{G}(z,\zeta,v_{1},v_{2}):=\begin{pmatrix}v_{2}z&v_{1}v_{2}z+\alpha(v_{1}^{2}-1)-\frac{1}{2}v_{1}v_{2}\\[4.30554pt] v_{1}v_{2}z+\alpha(v_{1}^{2}-1)-\frac{1}{2}v_{1}v_{2}&v_{2}(z-1)\end{pmatrix}

then a direct calculation shows that the deformation equation

\frac{\partial\hat{G}}{\partial z}+(\zeta W+\hat{H})J\hat{G}-\hat{G}J(\zeta W+\hat{H})=\hat{f}_{1}\,\frac{\partial H}{\partial v_{1}}+\hat{f}_{2}\,\frac{\partial H}{\partial v_{2}}+\hat{g}W

is fulfilled with the functions

\displaystyle\hat{f}_{1}:=(\zeta+v_{2})(1-v_{1}^{2}),\quad\hat{f}_{2}:=(\zeta v_{1}+v_{1}v_{2}+\tfrac{1}{2})v_{2},\quad\hat{g}:=4\alpha^{2}v_{1}(1-v_{1}^{2})-\zeta v_{1}v_{2}

From [4, Theorem 2.2] it follows that the eigenvalues $\zeta(v_{1},v_{2})$ satisfy the first-order quasilinear PDE

(\zeta+v_{2})(1-v_{1}^{2})\frac{\partial\zeta}{\partial v_{1}}+(\zeta v_{1}+v_{1}v_{2}+\tfrac{1}{2})v_{2}\frac{\partial\zeta}{\partial v_{2}}=4\alpha^{2}v_{1}(1-v_{1}^{2})-\zeta v_{1}v_{2}

(21)

Using Lemma 1 with $u_{3}=\Lambda(u_{1}+u_{2}^{2}-1)-2(\mu+1)u_{2}$ , we obtain that $\lambda=\Lambda(u_{1}+u_{2}^{2}-1)+(\mu+1)(\mu-2u_{2})$ is an eigenvalue of (2) for the parameters $\gamma^{2}=u_{1}$ , $\beta=-\Lambda(u_{1}+u_{2}^{2}-1)$ if and only if $\Lambda=0$ is an eigenvalue of $T(\mathbf{u})$ . In this case $\lambda=(\mu+1)(\mu-2u_{2})$ and $\beta=0$ , i.e., the CSWE reduces to the angular spheroidal wave equation (1). Note that $\Lambda=0$ is equivalent to $\zeta=-2\alpha v_{1}=-(\mu+1)v_{1}$ . Since $u_{1}=(1-v_{1}^{2})v_{2}^{2}$ and $u_{2}=v_{1}v_{2}$ , it follows that $\lambda=(\mu+1)(\mu-2v_{1}v_{2})$ is an eigenvalue of (1) for $\gamma^{2}=(1-v_{1}^{2})v_{2}^{2}$ if and only if $\zeta=-(\mu+1)v_{1}$ is an eigenvalue of $\hat{T}(v_{1},v_{2})$ .

Finally, let us simplify the PDE (21) even further. To this end, we need to modify the parameters introduced in (19) once more. We set

v_{1}=\tanh u,\quad v_{2}=\mathrm{e}^{t}\cosh u,\quad\zeta=\omega-\mathrm{e}^{t}\cosh u

With these parameters $(t,u)\in\mathbb{R}^{2}$ we convert the eigenvalue problem $\hat{\tau}\hat{y}=\zeta\hat{y}$ to

W(z)^{-1}\left(J\hat{y}^{\prime}(z)-\Phi(z,t,u)\hat{y}(z)\right)=\omega\hat{y}(z),\quad 0<z<1

(22)

where now $\omega$ is the eigenvalue parameter, and

\Phi(z,t,u):=\begin{pmatrix}2\mathrm{e}^{t}\cosh u+\frac{2\alpha\tanh u-\mathrm{e}^{t}\cosh u}{1-z}&-\frac{\alpha}{z}+\frac{\alpha}{1-z}+2\mathrm{e}^{t}\sinh u\\[4.30554pt] -\frac{\alpha}{z}+\frac{\alpha}{1-z}+2\mathrm{e}^{t}\sinh u&2\mathrm{e}^{t}\cosh u-\frac{2\alpha\tanh u+\mathrm{e}^{t}\cosh u}{z}\end{pmatrix}

Theorem 5.

Let $\omega=\omega(t,u)$ be an eigenvalue of the linear Hamiltonian system (22), and suppose that $\omega(t,u)$ depends analytically on the parameters $(t,u)\in S$ in some domain $S\subset\mathbb{R}^{2}$ . Then $\omega(t,u)$ solves the quasilinear partial differential equation

\frac{\partial\omega}{\partial t}+2\omega\cdot\frac{\partial\omega}{\partial u}=2(\mu+1)^{2}\frac{\tanh u}{\cosh^{2}u}+\mathrm{e}^{2t}\sinh(2u)+\mathrm{e}^{t}\cosh u

(23)

on $S$ . Moreover, $\lambda=(\mu+1)(\mu-2\mathrm{e}^{t}\sinh u)$ is an eigenvalue of the spheroidal wave equation (1) for $\gamma=\mathrm{e}^{t}$ if and only if the equation $\omega(t,u)=\mathrm{e}^{t}\cosh u-(\mu+1)\tanh u$ is satisfied for some $(t,u)\in S$ .

Proof.

Note that $\omega=\zeta+v_{2}$ , and that $v_{1}$ does not depend on $t$ . Hence,

	$\displaystyle\frac{\partial\omega}{\partial t}$	$\displaystyle=\frac{\partial\zeta}{\partial v_{2}}\frac{\partial v_{2}}{\partial t}+\frac{\partial v_{2}}{\partial t}=\mathrm{e}^{t}\cosh u\,\frac{\partial\zeta}{\partial v_{2}}+\mathrm{e}^{t}\cosh u$
	$\displaystyle\frac{\partial\omega}{\partial u}$	$\displaystyle=\frac{\partial\zeta}{\partial v_{1}}\frac{\partial v_{1}}{\partial u}+\frac{\partial\zeta}{\partial v_{2}}\frac{\partial v_{2}}{\partial u}+\frac{\partial v_{2}}{\partial u}=\frac{1}{\cosh^{2}u}\,\frac{\partial\zeta}{\partial v_{1}}+\mathrm{e}^{t}\sinh u\,\frac{\partial\zeta}{\partial v_{2}}+\mathrm{e}^{t}\sinh u$

and therefore

\displaystyle\frac{\partial\omega}{\partial t}+2\omega\cdot\frac{\partial\omega}{\partial u}

\displaystyle=\frac{2\omega}{\cosh^{2}u}\,\frac{\partial\zeta}{\partial v_{1}}+(\mathrm{e}^{t}\cosh u+2\omega\,\mathrm{e}^{t}\sinh u)\frac{\partial\zeta}{\partial v_{2}}+\mathrm{e}^{t}\cosh u+2\omega\,\mathrm{e}^{t}\sinh u

From

\frac{2\omega}{\cosh^{2}u}=2(\zeta+v_{2})(1-v_{1}^{2}),\quad\mathrm{e}^{t}\cosh u+2\omega\,\mathrm{e}^{t}\sinh u=2(\zeta v_{1}+v_{1}v_{2}+\tfrac{1}{2})v_{2}

and the partial differential equation (21) it follows that

	$\displaystyle\frac{\partial\omega}{\partial t}+2\omega\cdot\frac{\partial\omega}{\partial u}$	$\displaystyle=2(\zeta+v_{2})(1-v_{1}^{2})\frac{\partial\zeta}{\partial v_{1}}+2(\zeta v_{1}+v_{1}v_{2}+\tfrac{1}{2})v_{2}\frac{\partial\zeta}{\partial v_{2}}+2(\zeta v_{1}+v_{1}v_{2}+\tfrac{1}{2})v_{2}$
		$\displaystyle=8\alpha^{2}v_{1}(1-v_{1}^{2})-2\zeta v_{1}v_{2}+2(\zeta v_{1}+v_{1}v_{2}+\tfrac{1}{2})v_{2}$
		$\displaystyle=8\alpha^{2}v_{1}(1-v_{1}^{2})+(2v_{1}v_{2}+1)v_{2}=2(\mu+1)^{2}\frac{\tanh u}{\cosh^{2}u}+2\mathrm{e}^{2t}\sinh u\cosh u+\mathrm{e}^{t}\cosh u$

and this is exactly the PDE (23). As we have already noted, $\lambda=(\mu+1)(\mu-2v_{1}v_{2})=(\mu+1)(\mu-2\mathrm{e}^{t}\sinh u)$ is an eigenvalue of (1) for $\gamma^{2}=(1-v_{1}^{2})v_{2}^{2}=\mathrm{e}^{2t}$ if and only if $\omega=\zeta+v_{2}=-(\mu+1)\tanh u+\mathrm{e}^{t}\cosh u$ is an eigenvalue of $\hat{T}(v_{1},v_{2})$ and hence of (22). ∎

It is worth mentioning that the partial differential equation (23) is basically a forced inviscid Burgers’ equation. In fact, if we set $t=\frac{s}{2}$ , then (23) becomes

\frac{\partial\omega}{\partial s}+\omega\cdot\frac{\partial\omega}{\partial u}=(\mu+1)^{2}\frac{\tanh u}{\cosh^{2}u}+\tfrac{1}{2}\mathrm{e}^{s}\sinh(2u)+\tfrac{1}{2}\mathrm{e}^{\frac{s}{2}}\cosh u

and this quasilinear PDE can be written in the form $\frac{\partial\omega}{\partial s}+\omega\cdot\frac{\partial\omega}{\partial u}=f(s,u)$ . The forced inviscid Burgers’ equation is commonly encountered in fluid mechanics or gas dynamics, where it is used, for example, to study nonlinear interactions of dispersive waves. It is quite surprising that this PDE is related to the linear Hamiltonian system (22) and thus as well to the eigenvalues of the spheroidal wave equation.

4. Conclusion

Calculating the eigenvalues of a spheroidal wave equation is known to be a tricky task. There is an ongoing effort in finding new techniques to compute these eigenvalues in a convenient way and with high precision (see, for example, [1], [9], [10], [11], [12], [13] and the references therein). In the present work, we used an approach via a $2\times 2$ differential system. Apart from some numerical issues, we have mainly studied the analytical dependence of the eigenvalues on the parameters $\gamma$ and $\beta$ (or $u_{1}$ , $u_{2}$ , $u_{3}$ , respectively). Results such as those in Theorem 5 and relations like the PDE (5) may help to obtain new insights into the analytic properties of the angular or Coulomb spheroidal eigenvalues. The CSWE, in turn, is a special case of the confluent Heun differential equation (CHE), which is likewise a subject of current research. The methods presented here may also apply to the CHE and to even more general eigenvalue problems.

Appendix: A numerical example

In section 3, as an application of Theorem 4, it has been suggested to compute an eigenvalue of the spheroidal wave equation by starting from a zero of the function $\tilde{\Theta}$ and then following a characteristic curve of the PDE (5) until a zero of the function $\Lambda(t)$ is located. As a concrete numerical example, let us consider the function $\tilde{\Theta}(\Lambda):=\Theta(\Lambda,a,0,(a-1)\Lambda)$ in the case $\mu=1$ for $a=5$ . The limit formula in Lemma 3 provides the function shown in Fig. 1. It has a zero at $\Lambda_{0}=-0.8417200168449013$ , which can be determined, for instance, with the secant method. Here, all results are displayed rounded to $16$ decimal places, while the calculations were internally done with higher precision. Solving the differential system (18) for the initial values $u_{1}(0)=5$ , $u_{2}(0)=0$ , $\Lambda(0)=\Lambda_{0}$ by means of the Runge-Kutta method, we obtain $u_{1}(t)$ , $u_{2}(t)$ and $\Lambda(t)$ . The function $\Lambda(t)$ is illustrated in Fig. 2; it has a zero at $t_{0}=0.2793371978706399$ . Evaluating $u_{1}$ , $u_{2}$ at $t_{0}$ finally gives the eigenvalue $\lambda=2(1-2u_{2}(t_{0}))=-5.2736106330552739$ for the spheroidal wave equation (1) with parameter $\gamma^{2}=u_{1}(t_{0})=8.7417666942941543$ .

Refer to caption — Figure 1. The function $\tilde{\Theta}(\Lambda)$ for $\mu=1$ and $a=5$ with a zero at $\Lambda_{0}=-0.841720\ldots$

References

[1] P. E. Falloon, Theory and computation of spheroidal harmonics with general arguments, Master’s thesis, The University of Western Australia, Department of Physics (2001).
[2] T. Kereselidze, J. F. Ogilvie, The hydrogen-atom problem and Coulomb Sturmian functions in spheroidal coordinates, Vol. 77 of Advances in Quantum Chemistry, Academic Press, 2018, pp. 391–421. doi:10.1016/bs.aiq.2018.02.002.
[3] D. Batic, H. Schmid, M. Winklmeier, On the eigenvalues of the Chandrasekhar-Page angular equation, Journal of Mathematical Physics 46 (1) (2005) 012504. doi:10.1063/1.1818720.
[4] H. Schmid, On the deformation of linear Hamiltonian systems, Journal of Mathematical Analysis and Applications 499 (2) (2021) 125051. doi:10.1016/j.jmaa.2021.125051.
[5] J. Weidmann, Spectral Theory of Ordinary Differential Operators, vol. 1258 of Lecture Notes in Mathematics, Springer, Berlin – Heidelberg, 1987. doi:10.1007/BFb0077960.
[6] H. Sun, Y. Shi, Self-adjoint extensions for linear Hamiltonian systems with two singular endpoints, Journal of Functional Analysis 259 (2010) 2003–2027. doi:10.1016/j.jfa.2010.06.008.
[7] J. Meixner, F. W. Schäfke, Mathieusche Funktionen und Sphäroidfunktionen, Springer, New York, 1954, in German. doi:10.1007/978-3-662-00941-3.
[8] T. Kato, Perturbation theory for linear operators, 2nd Edition, vol. 1258 of Lecture Notes in Mathematics, Springer, Berlin – Heidelberg – New York, 1995. doi:10.1007/978-3-642-66282-9.
[9] H. Schmid, Computation of the eigenvalues for the angular and Coulomb spheroidal wave equation, Applied Numerical Mathematics 185 (2023) 101–119. doi:10.1016/j.apnum.2022.11.018.
[10] C. Flammer, Spheroidal Wave Functions, Stanford University Press, Stanford, CA, 1957.
[11] D. B. Hodge, Eigenvalues and eigenfunctions of the spheroidal wave equation, Journal of Mathematical Physics 11 (8) (1970) 2308. doi:10.1063/1.1665398.
[12] P. Kirby, Calculation of spheroidal wave functions, Computer Physics Communications 175 (7) (2006) 465–472. doi:10.1016/j.cpc.2006.06.006.
[13] S. L. Skorokhodov, Calculation of eigenvalues and eigenfunctions of the Coulomb spheroidal wave equation, Matematicheskoe modelirovanie 27 (7) (2015) 111–116, in Russian.