TY - JOUR
AU - Daniel, Derek, J.
AB - Abstract Mathieu’s equation originally emerged while studying vibrations on an elliptical drumhead, so naturally, being a linear second-order ordinary differential equation with a Cosine periodic potential, it has many useful applications in theoretical and experimental physics. Unfortunately, there exists no closed-form analytic solution of Mathieu’s equation, so that future studies and applications of this equation, as evidenced in the literature, are inevitably fraught by numerical approximation schemes and nonlinear analysis of so-called stability charts. The present research work, therefore, avoids such analyses by making exceptional use of Laurent series expansions and four-term recurrence relations. Unexpectedly, this approach has uncovered two linearly independent solutions to Mathie’s equation, each of which is in closed form. An exact and general analytic solution to Mathieu’s equation, then, follows in the usual way of an appropriate linear combination of the two linearly independent solutions. 1. Introduction: theory and applications 1.1. Mathieu’s equation The problem of vibrations of an elliptical membrane was first seriously studied by the French mathematician Émile Léonard Mathieu in 1868. In his celebrated paper, “Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique”, Mathieu proved that the vertical displacement produced from oscillations on an elliptical drumhead is governed, in general, by the 2D Helmholtz equation [1]. Furthermore,Mathieu showed that when the 2D Helmholtz equation is transformed into elliptical co-ordinates, the differential equations, respectively, for the radial and angular coordinates of the elliptical system become separable. The result is two almost identical linear second-order differential equations, each of which is aptly named the “Mathieu equation”. So, to this end, Mathieu’s equation has the canonical form $$\begin{align}\label{eq1.1} \frac{d^2y(x)}{dx^2}+\left( {a-2q\cos \left( {2x} \right)} \right)y(x)=0, \end{align}$$(1.1) where the separation constant |$a$| and the physical parameter |$q$| are both real. As impressed in McLachlan’s detailed exposition of the theory and applications of Mathieu’s equation [2], the appearance of the cosine trigonometric function in Eq. (1.1) of Mathieu’s equation is characteristic of describing a wide range of phenomena in nonlinear vibrations theory. Various kinds of mechanical and structural dynamical systems associated with elliptical geometries, for instance, can accurately and realistically be modeled by Mathieu’s equation, such as wave propagation in pipes, electromagnetic wave guides, or oscillations of water in a lake, as well as floating vessels, to name a few [3]. Ruby, for instance, has provided several other interesting applications of Mathieu’s equation in the context of parametric oscillators [3], including a rather elaborate sketch on the wave motion and stability of floating bodies or vessels. Here, the “metacenter” of a floating body, which is the line that connects the center of buoyancy with the center of gravity of the vessel, has the ability to tilt at an angle that can suddenly undergo sinusoidal motion when a wave begins to strike the floating vessel. Furthermore, provided that damping is ignored, this situation can be shown to parallel that of a driven harmonic or parametric oscillator, whose corresponding equations of motion are consecutive with Mathieu’s equation in Eq. (1.1). In effect, Mathieu’s equation can be used meticulously to measure the “roll” stability of floating vessels or ships, i.e., the possibility of a floating vessel or ship to capsize. The Paul trap, as devised and named after Wolfgang Paul [4], or a quadrupole ion trap for charged particles, is an interesting example of a parametric oscillator. As the name suggests, charged particles can be trapped by an electric field by having their motion confined along the direction of an oscillating electric field. Once again, the resulting equations of motion lead directly to Mathieu’s equation [3], as given by Eq. (1.1). Even though Paul traps are commonly found in mass spectrometers, their real appeal is in “quantum computers”, on account that ion traps are required for storing qubits of quantum information, as originally proposed by Cirac and Zoller [5]. Bearing in mind, of course, that the present aim is to derive a general and exact solution of Mathieu’s equation, as defined in Eq. (1.1), most of the attempts in the literature [2,3], however, have yet to succeed in finding such a solution for Mathieu’s equation. For this reason, the majority of scientific efforts rests on using the most reliable numerical and approximate methods available. For example, some of the reliable methods are discussed in Ref. [6], as well as in Ref. [7] and references therein, with the latter authors also having gone to some length to provide a compendium of data tables for computing eigenvalues of Mathieu’s equation. The status of the various kinds of mathematical developments that are now laying at the “heels” of Mathieu’s equation will also be briefly mentioned next. 1.2. The mathematics The earliest historical mathematical developments of Mathieu’s equation, as defined by Eq. (1.1), following the power series expansion method of solution originally used by Mathieu in 1868 [1], has been traced out in great detail by McLachlan [2], who also references the major works of many of the original contributors over the period from 1878 to 1927. This period starts with Heine’s work in 1878 on the sine and cosine series solutions of integral order, in which infinite continued fractions were first used to determine the characteristic values or eigenvalues, |$a$|⁠. Several years later, Floquet developed a general method to treat periodic coefficients, with Floquet’s theorem emerging in 1883 as a result. The introduction of infinite determinants, as published by Hill in 1886, appeared as the next major achievement in expressing solutions to Mathieu’s equation and its generalization. But around this time, the matter of stability in the computations of the characteristic parameter, |$a,$| because of its inextricable dependence on the amplitude parameter, |$q,$| required a sharper approach. Therefore, the inception of the “stability chart” was introduced by Ince in 1927 as a useful apparatus and graphical tool to investigate the convergence behavior of the characteristic parameter |$a,$| which is known to vary as a function of the amplitude parameter |$q.$| This approach has been carefully reviewed by Kovacic et al. [8], where their use of the stability chart, also known as a “Strutt diagram” or “Strutt–Ince diagram,” has become an integral component of any computations associated with Mathieu’s equation. These authors have also provided some material from the original contributors (see also the appendices in Ref. [8]) along with a very effective timeline of the historical developments behind Mathieu’s equation, which include the tuning fork problem as one of the first applications of Mathieu’s equation by Lord Rayleigh in 1887. The boundary conditions of Mathieu’s equation in Eq. (1.1) are given by $$\begin{align}\label{eq1.2} y(x)=y(x+\pi), \end{align}$$(1.2) which reflects the periodic nature of the trigonometric term, |$\cos 2x,$| in Mathieu’s equation. Mathieu’s equation, in general, does not possess periodic solutions. But because the coefficients of Mathieu’s equation are periodic functions of |$x,$| an important result, known as Floquet’s theorem [7], $$y(x+\pi )=e^{i\nu x}y(x),$$ with the characteristic exponent, |$\nu ,$| being any real or complex constant, allows two independent solutions to be admitted [2]. Further, if |$\nu $| is a real number and is an integer, |$y(x)$| will be periodic with period |$\pi $| or |$2\pi ,$| as a consequence. Therefore, the appropriate solutions to Mathieu’s equation are those that are periodic in |$x$| with period |$\pi $| or |$2\pi ,$| and so are referred to in the literature as the Mathieu functions (see, e.g., Ref. [7]). As pointed out by Morse and Feshbach, any attempt to solve Mathieu’s equation directly will be confronted, of course, by a three-term recursion relation [6]. To examine this more closely, let the solution |$y(x)$| in Eq. (1.1) be represented by a Laurent series, or, equivalently, have the Fourier expansion $$y(x)=\sum\limits_{n=-\infty }^\infty {\,a_n^ {\kern 1pt}} {\kern 1pt}e^{2in{\kern 1pt}x},$$ which when substituted into Eq. (1.1) leads to the so-called three-term recursion, $$\begin{align} \label{eq1.3} qa_{n+1}=\left( {a-4n^2} \right)a_n +qa_{n-1}=0, \end{align}$$(1.3) for the unknown coefficients |$a_n$|⁠. Clearly, solving this three-term recursion relation is going to be far more complex than solving a two-term recursion formula. Fortunately, solving a three-term recursion relation, as shown by Morse and Feshbach, can be rendered rather effectively by the method of infinite continued fractions [6]. The formal aspects behind continued fractions begin the solving process of a three-term recursion relation by first deriving expressions directly from Eq. (1.3) for the ratios |$a_n/a_{n-1}$| and |$a_{n+1}/a_n$|⁠. After choosing some starting values for |$a_0$| and |$a_1$|⁠, an iterative process for computing the values of |$a_n/{a_{n-1}}$| and |$a_n/a_{n+1}$| can be initiated, respectively, for large positive and large negative values of |$n.$| If the values, |$a_n/a_{n-1}$| and |$a_n/a_{n+1}$|⁠, approach zero as |$n\to \pm \infty $|⁠, then there is every chance that the series as displayed above for |$y(x)$| will converge over the range |$0 < x<\infty .$| In practice, successive approximations of |$a_n/a_{n-1}$| and |$a_n/a_{n+1}$| will generate a solution for the eigenvalues |$a$| in Eq. (1.3) as an infinite continued fraction of |$a$| that varies as a function of |$q.$| This entire procedure is exemplified by Eq. (6) of Ref. [3], and has the capacity to extract an infinite sequence of roots of |$a$| as a function of |$q,$| thereby reproducing the appropriate Mathieu functions, i.e., both the odd and even periodic solutions of Mathieu’s equation. 1.3. Quantum mechanics While the classical physics of the Helmholtz equation, on the one hand, is contingent with most of the studies concerning the diffraction of light beams, the mathematics of the Helmholtz equation, on the other hand, has a solid basis with the Schrödinger equation in quantum mechanics [2]. This well-known connection with the Schrödinger equation and subsequent tie mentioned in Sect. 1.1 between the Helmholtz equation and Mathieu’s equation offers an opportunity here of exploring optical or light beams in the context of a quantum pendulum [9]. In fact, the earliest treatment of Mathieu’s equation for studying applications in quantum mechanics was carried out in 1928 by Edward Condon [10], who was one of the first pioneers in atomic physics to popularize the quantum pendulum in the literature. More accurately, he showed that Mathieu’s equation conforms absolutely to the physical pendulum in quantum mechanics. Basically, the quantum pendulum is described by a wavefunction |$\Psi$|⁠, satisfying a time-independent Schrödinger equation for the motion of a point mass |$m$| that is attached at the end of a massless rigid rod of length |$l,$| which is constrained to move in a circle under the influence of constant gravity, as specified by |$g.$| The time-independent Schrödinger equation for the angular displacement |$\theta $| of the pendulum is given by $$\begin{align}\label{eq1.4} \frac{-\hbar ^2}{2ml^2}\frac{d^2\Psi }{d\theta ^2}+V{\kern 1pt}{\kern 1pt}\Psi =E\Psi, \end{align}$$(1.4) where |$\hbar =h/2\pi $| is Planck’s constant |$h$| divided by |$2\pi , E$| is the energy eigenvalue, and |$V$| is the potential energy, as defined by $$\begin{align*} V=mgl(1-\cos \theta). \end{align*}$$ Use of the variable, |$x=\tfrac{1}{2}\theta ,$| and the abbreviations $$\varepsilon =\frac{8ml^2E}{\hbar ^2}-\frac{8m^2gl^3}{\hbar ^2},\quad \omega _{{\kern 1pt}0}^{{\kern 1pt}2} =\frac{8m^2gl^3}{\hbar ^2},$$ converts the Schrödinger equation in Eq. (1.4) into $$\begin{align}\label{eq1.5} \frac{d^2{\kern 1pt}\Psi (x)}{dx^2}+\left( {\varepsilon +\omega {\kern 1pt}_0^2 \cos \left( {2x} \right)} \right)\Psi (x)=0. \end{align}$$(1.5) Comparing this equation directly with Mathieu’s equation in Eq. (1.1) shows that |$a$| is related to the energy of the pendulum through the relation |$a=\varepsilon ,$| while |$q$| is related to the barrier height |${\omega _{0}^{2} } / 2$| of the cosine potential well, i.e., |$q=-{\omega _{0}^{2} } / 2,$| which asserts an acquaintance between Mathieu’s equation in Eq. (1.1) and the quantum pendulum in Eq. (1.5). Although Schrödinger’s equation for the pendulum, as defined by Eq. (1.5), falls into a class of trigonometric potentials [11] in quantum mechanics that cannot be solved for analytically, unlike the Pöschl–Teller [12] potential in quantum mechanics, which can be solved through the use of special functions, it will be shown, all the same, by the present author that exact and analytical solutions to Mathieu’s equation can be found.And these solutions, as expected, also necessarily solve the quantum pendulum problem in Eq. (1.5). 1.4. The method Since the preceding discussion in Sect. 1.3 makes it clear that solving the equation for the quantum pendulum in Eq. (1.5) is interchangeable with solving Mathieu’s equation in Eq. (1.1), there will be no loss of generality with the choice of adopting the time-independent Schrödinger equation in Eq. (1.5), as well as accompanying notation, as the starting point in the present analysis that follows. Consider first a change of variable in Eq. (1.5) from the independent variable |$x$| to the complex variable |$z=e^{ix}.$| In terms of the complex variable |$z,$| the cosine function in Eq. (1.5) will carry on in the form $$\cos 2x=\tfrac{1}{2}\left( {z^2+z^{-2}} \right)\!,$$ while the corresponding time-independent Schrödinger equation in Eq. (1.5) transforms as $$z^2\frac{d^2{\kern 1pt}\Psi (z)}{dz^2}+z\frac{d{\kern 1pt}{\kern 1pt}\Psi (z)}{dz}-\left[ {\frac{\omega _0^2 {\kern 1pt}}{2}\left( {z^2+\frac{1}{z^2}} \right)+\varepsilon } \right]\Psi (z)=0.$$ To see if there is more than one way to solve this second-order differential equation, it will be more useful, tactfully, to contemplate the following ansatz for |$\Psi (z),$| $$\begin{align}\label{eq1.7} \Psi (z)=z^p\exp \left( {-\frac{\alpha }{z}+\beta z} \right)f(z), \end{align}$$(1.7) where |$p, \alpha$|⁠, and |$\beta$| are arbitrary free parameters yet to be determined. Consequently, the problem has now shifted to solving the following second-order differential equation for the auxiliary function |$f(z)$| in Eq. (1.7): $$\begin{align}\label{eq1.8} f'' (z)+\left[ {\frac{2\alpha }{z^2}+\frac{\delta }{z}+2\beta } \right]{f}'(z)+\left[ {\left( {\beta ^2-\tfrac{1}{2}\omega ^2} \right)+\frac{\beta \delta }{z}+\frac{\gamma }{z^2}+\frac{\alpha (\delta -2)}{z^3}+\frac{\left( {\alpha ^2-\tfrac{1}{2}\omega ^2} \right)}{z^4}} \right] f(z)=0, \end{align}$$(1.8) where the prime |$'$| denotes differentiation with respect to the |$z$| variable, with |$\delta =2p+1$| and |$\gamma =p^2+2\alpha \beta -\varepsilon$|⁠. Curiously, the second-order differential equation defined in Eq. (1.8) is known to be a special case of Heun’s equation, which has been well studied by numerous authoritative researchers in Ref. [13] as a second-order differential equation with its various confluent forms. In comparison, it can be seen from Part C of Ref. [13] that Eq. (1.8) is related to the double-confluent Heun equation (DCHE), making it clear that the special and confluent forms of Heun’s equation can go by the names of Mathieu, Lamé, and Coulomb spheroidal equations, to indicate a few. General approaches and expansions for solving Heun’s equation can also be found in Ronveaux and chief references mentioned in Ref. [13], though fundamental solutions concerned specifically with the confluent and double-confluent Heun equations have also been described in some detail by El-Jaick and Figueiredo [14]. The latter authors explore the possibility of using both a formal power series as well as the confluent hypergeometric functions for solving the DCHE, which leads to an infinite continued fraction analysis of a three-term recursion relation for the expansion coefficients. Although the present work maintains the aim of deriving similar recursion relations from the corresponding power series expansions for solving Mathieu’s equation, the method adopted here, however, will differ strictly in its approach by avoiding three-term recursion relations and infinite continued fractions. In anticipation of offering flexibility in the analysis to be carried out below, Eq. (1.8) has been purposely armed sufficiently with a number of free parameters for the single purpose of seeking an analytical solution. The following functions will also appear from time to time in the analysis below. The hypergeometric function of general order $${}_pF_q \left( {a_1 ,\cdots ,a_p ;b_1 ,\cdots ,b_q ;z} \right)=\sum\limits_{k=0}^\infty {\frac{\left( {a_1 } \right)_k \cdots \left( {a_p } \right)_k }{\left( {b_1 } \right)_k \cdots \left( {b_q } \right)_k }} \frac{z^k}{k!},$$ where |$(a)_n$| is the Pochhammer symbol, defined as follows: $$\begin{array}{l} (a)_n =a(a+1)\cdots (a+n-1)=\frac{\Gamma (a+n)}{\Gamma (a)},\quad n\geqslant 1,\\ (a)_0 =1.\\ \end{array}$$ 2. General solution of Mathieu’s equation 2.1. First independent solution The first and last term inside the last square brackets of Eq. (1.8) can be made to vanish by choosing |$\alpha =\pm \omega / {\sqrt 2 }$| and |$\beta =\pm \omega / {\sqrt 2 },$| while the |$z^{-3}$| will vanish by assuming |$\delta =2,$| which is equivalent to assigning |$p=\tfrac{1}{2}.$| It is expedient to multiply Eq. (1.8) throughout by |$z^2,$| followed by a change of the independent variable from |$z$| to |$u$| via the substitution |$z=r-su,$| where |$r$| and |$s$| are arbitrary free parameters. Equation (1.8) now becomes a second-order differential equation with respect to the variable |$u$|⁠, $$\begin{align}\label{eq2.1} \frac{1}{s^2}\left( {r^2-2rsu+s^2u^2} \right)\frac{d^2f(u)}{du^2}-\frac{1}{s}\left( {A-Bu+Cu^2} \right)\frac{df(u)}{du}+\left( {D-Eu} \right)f(u)=0, \end{align}$$(2.1) where $$\begin{align*} &A=2\alpha +\delta r+2\beta r^2, \quad B=\delta s+4\beta rs, \quad{} C=2\beta s^2,\\ & D=\beta \delta r+\gamma, \quad E=\beta \delta s. \end{align*}$$ One seeks a solution of Eq. (2.1) by assuming that |$f(u)$| has a Laurent series expansion, $$f (u) = \sum_{n=-\infty}^\infty \Gamma (n+\nu)\,a_n u^{-n-\nu},$$ where |$a_n$| are expansion coefficients and |$\nu $| is arbitrary. Substitution of |$f(u)$| into Eq. (1.8) and equating to powers of |$u^{-n-\nu}$| leads, after some algebra, to a four-term recursion relation for |$a_n$|⁠, $$\begin{align}\label{eq2.2} & \sum_{n=-\infty}^\infty \Gamma (n+\nu) u^{-n-\nu} \,\left\{\left(\frac{C}{s}(n+\nu +1)(n+\nu)-E(n+\nu)\right) a_{n+1}\right.\nonumber\\ &\quad \quad \quad \quad \quad \quad +\left((n+\nu +1)(n + \nu)-\frac{B}{s}(n+\nu) + D\right) a_n +\nonumber\\ &\quad \quad \quad \quad \quad \quad \left. + \left(-\frac{2r}{s} (n+\nu) + \frac{A}{s}\right) a_{n-1} +\frac{r^2}{s^2} a_{n-2}\right\} = 0 \end{align}$$(2.2) It can be seen inside the curly bracket of Eq. (2.2) that a four-term recursion relation for |$a_n$| has appeared instead of the usual three-term recursion relation that was encountered in Eq. (1.3) for Mathieu’s equation. To avoid repeating the analysis of infinite continued fractions for solving a four-term recursion relation, which can be rather more tedious than that for a three-term recursion relation discussed in Sect. 1.2, one can look, instead, into the prospect of matching the recursion relation in Eq. (2.2) with a known four-term recursion formula. For instance, there are four-term recursion formulas that have already been derived by Exton [15] for many of the well-known hypergeometric functions. More specifically, Exton has derived the four-term recursion formula for the hypergeometric function |${}_3 F_1$| (Eq. (4.8) of Ref. [15]): $$\begin{align}\label{eq2.3} &\left(n^2+(c+1)n+c\right) P_{n+1} -\left(n^2+(a+b+2t)n+ct+ab\right)P_n\nonumber\\ &\quad{} +\left({2t{\kern 1pt}n+a+b+1-2t+t^2} \right)P_{n-1} -t^2P_{n-2} = 0, \end{align}$$(2.3) where |$P_n $| is defined in terms of the hypergeometric function |${ }_3F_1 $|⁠: $$\begin{align}\label{eq2.4} P_n =P_n \left(a,b;c;t\right)=t^n\frac{{ }_3F_1 \left( {a,b,-n;c;-\tfrac{1}{t}} \right)}{n!}. \end{align}$$(2.4) Note, for the sake of clarity, that the |$t$| appearing in Eq. (2.3) replaces the original variable |$x$| used in Eq. (4.8) of Ref. [15]. It is a matter of now proving that the functions |$P_n $| (up to an overall sign change) in the four-term recursion formula of Eq. (2.3) will exactly solve the four-term recursion relation for |$a_n$| in Eq. (2.2), subject to the condition that |$t$| can be treated as a parameter whose value can be fixed. If |$a_n$| is to be precisely given by |$P_n ,$| as defined by Eq. (2.4), then each of the terms in the round parentheses associated with the four-term recursion relation for |$a_n$| in Eq. (2.2) has to be matched identically with the corresponding terms in the round parentheses appearing in the recursion formula for |$P_n$| in Eq. (2.3), as demonstrated below. Consider applying the definitions for |$A,B,C,D$|⁠, and |$E,$| as given below Eq. (2.1), to the first set of terms in the round brackets before |$a_{n+1}$| in Eq. (2.2), which can be shown to simplify as $$\left(\frac{C}{s}(n+\nu +1)(n+\nu )-E(n+\nu)\right) a_{n+1} \to \;\;2\beta s\left(n^2+2\nu n+\nu ^2\right)a_{n+1}.$$ For this result to be identical (up to an overall sign change) with the first set of terms contained in the round brackets preceding |$P_{n+1} $| in Eq. (2.3), one demands that |$c+1=2\nu \;$| and |$\nu ^2=c.$| The latter two relations immediately give |$c=1,$| followed by |$v=1$| and |$s=-1/{2\beta}$|⁠. Similarly, one can simplify the second set of terms inside the round brackets of the recursion coefficient |$a_n$| in Eq. (2.2): $$\begin{align*} &\left( {(n+\nu +1)(n+\nu )-\frac{B}{s}(n+\nu )+D} \right)a_n\\ &\quad{} \to \left(n^2+(2\nu +1-\delta -4\beta r)n+\nu (\nu +1)-\delta -4\beta r+\delta \beta r+\gamma\right) a_n. \end{align*}$$ It then follows, after matching this result with the second set of terms inside the round brackets preceding |$P_n $| in Eq. (2.3), that |$t=r/s,$| since |$v=1$| and |$\delta =2$| were fixed at the outset. The sum of |$a$| and |$b$| can then be extracted next: $$a+b+2t=1+\frac{2r}{s}\quad \Rightarrow \quad a+b=1.$$ But it was shown above that |$c=1,$| so the product of |$a$| and |$b$| must satisfy $$ct+ab=\frac{r}{s}+\gamma \quad \Rightarrow \quad ab=\gamma .$$ Applying the definition below Eq. (1.1) for the |$\gamma $| term and remembering from the start of this section that |$\delta =2,$| which is equivalent to setting |$p=\tfrac{1}{2},$| permits |$\gamma $| to be expressed as $$\gamma =2\alpha \beta -\varepsilon +\frac{1}{4}.$$ Combining the previous three relations obtained above for |$a,b,$| and |$\gamma $| enables |$a$| and |$b$| to be solved explicitly: $$a^2-a+\gamma =0\quad \Rightarrow \quad a=\frac{1}{2}\pm \sqrt {\varepsilon -2\alpha \beta } \quad \;\text{and}\quad \text{ }b=1-a=-\frac{1}{2}\mp \sqrt {\varepsilon -2\alpha \beta } .$$ Here, |$a$| and |$\text{ }b$| can always be kept as real constants, provided that the signs of |$\alpha $| and |$\beta $| are taken opposite to one another. As deduced above, |$s=-1 /{2\beta}$|⁠, so the third set of terms inside the round brackets preceding |$a_{n-1} $| in Eq. (2.2) simplifies as $$\left(\frac{-2r}{s}(n+\nu )+\frac{A}{s}\right) a_{n-1} \to -\left(\frac{2r}{s} n-\frac{\left(2\alpha +\delta r+2\beta r^2-2r\nu\right)}{s}\right)a_{n-1}.$$ Comparing this expression with the third set of terms in the round brackets preceding |$P_{n-1} $| in Eq. (2.3), and remembering that |$t=r/s$|⁠, advances a relation for $$a+b+1-2t+t^2=-\frac{2\alpha }{s}-\frac{2\beta r^2}{s}\quad \Rightarrow \quad a+b+1-2t=-\frac{2\alpha }{s},$$ which will readily determine the final constant, |$r.$| Simply use the fact established above that |$a+b=1$|⁠, as well as |$s=-1 / {(2\beta )},$| to obtain $$2-2t=-\frac{2\alpha }{s}\quad \Rightarrow \quad r=\frac{2\alpha \beta -1}{2\beta }.$$ All of the unknown constants in Exton’s four-term recursion formula have now been explicitly determined and identified with all of the corresponding arbitrary parameters contained in the four-term recursion relation for |$a_n ,$| the result of which is summarized below: $$t=\frac{r}{s}=1-2\alpha \beta ,\quad a=\frac{1}{2}\pm \sqrt {\varepsilon -2\alpha \beta } ,\quad b=-\frac{1}{2}\mp \sqrt {\varepsilon -2\alpha \beta } ,\;\quad \text{and}\quad c=1.$$ This proves that the four-term recursion relation for |$a_n$| in Eq. (2.2) has a strict identity (up to an overall sign change) with the four-term recursion formula of Eq. (2.3) for |$P_n .$| Subsequently, there can be no further objection in stating that the recursion coefficients |$a_n$| are unmistakably determined by Eq. (2.4) for |$P_n$|⁠, $$a_n =t^n \frac{{ }_3F_1 \left( {a,b,-n;c;-\tfrac{1}{t}} \right)}{n!}.$$ It is worth pausing here in order to affirm that the range of summation of |$n$| from |$-\infty $| to |$\infty $| in the Laurent series expansion of |$f(u)$| is in no way affected by the |$n!$| appearing in the denominator of the above result for |$a_n .$| Simply insert the value |$\nu =1,$| as deduced from the first step above, into the gamma function |$\Gamma (n+\nu )$| appearing in the Laurent series expansion of |$f(u)$| to see that |$\Gamma (n+\nu )=n!,$| which clearly cancels the |$n!$| contained in the denominator above for |$a_n$|⁠. The analytical form above for the expansion coefficients |$a_n$| can be inserted back into the Laurent series expansion of |$f(u),$| which in combination with the ansatz for |$\Psi (z),$| defined in Eq. (1.7), together with |$(r-z)/s$| in place of |$u,$| i.e., returning back to the |$z$| variable, will lead to an explicit form of the desired solution to Mathieu’s equation. It follows from the preceding remarks that the first independent solution, denoted here by |$\Psi _I (x),$| can be written in the analytical form |$(z=e^{ix})$|⁠, $$\begin{align}\label{eq2.5} \Psi_I (x)=z^{-\tfrac{1}{2}}\exp \left( {-\frac{\alpha }{z}+\beta z} \right)\sum\limits_{n=-\infty }^\infty \, \left( {\frac{t}{s}} \right)^n{ }_3F_1 \left({a,b,-n;c;-\tfrac{1}{t}} \right)\left( {r-z} \right)^{-n}, \end{align}$$(2.5) which is an exact solution to both Mathieu’s equation, defined in Eq. (1.1), and the quantum pendulum, defined by Eq. (1.5). Here, the corresponding parameters follow as: $$\begin{align*} \alpha =\pm \frac{\omega}{\sqrt 2},\quad \beta =\pm \frac{\omega}{\sqrt 2}, \quad a=\frac{1}{2}\pm \sqrt {\varepsilon -2\alpha \beta}, \quad b=\frac{1}{2}\mp \sqrt {\varepsilon -2\alpha \beta},\\ c=1,\quad {t=\frac{r}{s},} \quad {s=-\frac{1}{2\beta },} \quad {r=\frac{2\alpha \beta -1}{2\beta},}\\ \end{align*}$$ 2.2. Second independent solution Let |$\alpha =\pm \omega / {\sqrt 2 },$| but leave |$\beta $| and |$\delta $| arbitrary in Eq. (1.8). Next, assume that the expansion $$f(z)=\sum\limits_{n=-\infty }^\infty {(-1)^n\,\lambda ^na_n } z^{1-n},$$ is a solution of the second-order differential equation (1.8), where |$a_n$| are the expansion coefficients and |$\lambda $| is a free parameter. It can be shown, if the above expansion for |$f\left( z \right)$| is substituted into the second-order differential equation of Eq. (1.8) and then multiplied throughout by |$z^3,$| that a four-term recursion relation for |$a_n$| appears: $$\begin{align}\label{eq2.6} & \sum_{n=-\infty}^\infty (-1)^n\,\lambda^{n+1} z^{1-n}\,\left\{\!\!\vphantom{\frac{2\alpha}{\lambda}} -\lambda ^2\left( {\beta ^2{-}\tfrac{1}{2}\omega ^2} \right)a_{n+3}{-}\lambda \beta \left( {2n+2{-}\delta } \right)a_{n+2} {-}\left( {n^2+(1-\delta )n+\gamma } \right)a_{n+1}\right.\nonumber\\ &\left.\quad{} - \frac{2\alpha}{\lambda}\left( {n-\tfrac{1}{2}\delta} \right) a_n\right\} = 0. \end{align}$$(2.6) The aim here is the same as in Sect. 2.1, which is to demonstrate that there is a known four-term recursion formula that can be shown to satisfy the recursion relation for the expansion coefficients |$a_n$| appearing in Eq. (2.6). Fortunately, one can appeal again to another one of Exton’s four-term recursion formulas, as given by Eq. (3.7) in Ref. [15], which this time is specified by the hypergeometric function |${ }_2F_2 $| as follows: $$\begin{align}\label{eq2.7} & t^2G_{n+3} +t\left( {2n+2-c-t} \right)G_{n+2} -\left( {n^2+\left( {c+1-2t} \right)n+c+\left( {1+a+b} \right)t} \right)G_{n+1}\nonumber\\ &\quad +\left( {\left( {1+a+b} \right)n-ab} \right)G_n = 0, \end{align}$$(2.7) where |$G_n $| is defined in terms of the hypergeometric function |${ }_2F_2 $|⁠: $$\begin{align}\label{eq2.8} G_n =G_n \left(a,b;c;t\right)=\frac{\left( a \right)_n \left( b \right)_n }{\left(c \right)_n n!}{ }_2 F_2 \left(a+n,b+n;c+n,1+n;t\right)\!. \end{align}$$(2.8) Consider now comparing the recursion identity for |$G_n $| in Eq. (2.7) side by side with the recursion relation defined in Eq. (2.6) for |$a_n$|⁠, then one needs to concentrate on matching the unknown constants in the recursion formula for |$G_n$| with those in the recursion relation for |$a_n$|⁠, as follows: $$\begin{align}\label{eq2.9} 2-c-t & = 2-\delta\nonumber\\ c+1-2t & = 1-\delta\nonumber\\ c+\left(1+a+b\right)t& = \gamma\nonumber\\ 1+a+b & = -\frac{2\alpha}{\lambda}\nonumber\\ -ab & = \frac{\delta \alpha}{\lambda}. \end{align}$$(2.9) Observe, from the relations above, that |$t$| can be solved first: $$t=\tfrac{2}{3}\delta ,$$ as seen by adding the first two relations in Eq. (2.9). Likewise, |$c$| is determined by substituting the above result for |$t,$| say, into the first relation of Eq. (2.9): $$\begin{align*} c=\tfrac{1}{3}\delta \Rightarrow c=\frac{t}{2}. \end{align*}$$ There are two more relations involving |$t,$| however. This can be seen by comparing the two leading terms in each of the four-term recursion relations for |$a_n$| and |$G_n ,$| defined respectively in Eqs. (2.6) and (2.7). That is, $$t=-\lambda \beta \text{ and }t^2=-\lambda ^2\left( {\beta ^2+\tfrac{1}{2}\omega ^2} \right)\!.$$ Obviously, for these two relations to hold requires $$\beta =\pm \frac{\omega}{2}.$$ This leaves the matter of sorting out how |$t=-\lambda \beta $| and |$t=\tfrac{2}{3}\delta $| can be satisfied simultaneously. Fortunately, |$\lambda $| is a free parameter, which means that |$t$| will be unique, so long as |$\lambda $| is chosen to satisfy $$t=-\lambda \beta \quad \text{and}\quad t=\tfrac{2}{3}\delta \quad \Rightarrow \quad \lambda =-\frac{2\delta }{3\beta }.$$ Before sorting out the three remaining relations in Eq. (2.9), consider first the above definitions given at the start for |$\delta $| and |$\gamma ,$| namely, $$\delta =2p+1\text{ and }\gamma =p^2+2\alpha \beta -\varepsilon .$$ The first relation here will serve its purpose better when recast as |$p=\tfrac{1}{2}\left( {\delta -1} \right)\!,$| so that |$\gamma ,$| in turn, can be rewritten as $$\begin{align*} \gamma =\frac{\delta ^2}{4}-\frac{\delta }{2}+2\alpha \beta -\varepsilon +\frac{1}{4}. \end{align*}$$ This form of |$\gamma $| commits itself to the intention of solving a quadratic equation involving |$\delta .$| First express the left-hand side of the third relation in Eq. (2.9) as $$c+\left( {1+a+b} \right)t=c+2\alpha \beta,$$ which results from using the fourth relation in Eq. (2.9), together with |$t=-\lambda \beta .$| But it was also established above that |$c=\tfrac{1}{3}\delta ,$| so it too can be applied here directly to give $$c+\left( {1+a+b} \right)t=\gamma \quad \Rightarrow \quad \tfrac{1}{3}\delta +2\alpha \beta =\gamma.$$ Inserting into the right-hand side of this relation the expression for |$\gamma $| derived above, involving |$\delta $| and |$\delta ^2,$| produces the sought-after quadratic equation for |$\delta ,$| with corresponding roots determined as $$\delta =\frac{10\pm 4\sqrt {9\varepsilon +4} }{3}.$$ The last two relations in Eq. (2.9) are all that remain, which will be accommodated together by applying the earlier result |$\lambda =-{2\delta} / {3\beta }$| above, like so: $$\begin{align*} 1+a+b & = \frac{3\alpha \beta}{\delta},\\ -ab & = -\frac{3\alpha \beta}{2}. \end{align*}$$ These two equations make it straightforward to solve for |$a$| and |$b$|⁠, respectively: $$\begin{align*} a & = \frac{-\left(\delta -3\alpha \beta\right) \pm \sqrt{\left(\delta -3\alpha \beta\right)^2 - 6 \alpha \beta \delta^2}}{2\delta},\\ b & =\frac{-\left(\delta -3\alpha \beta\right) \mp \sqrt{\left( {\delta -3\alpha \beta } \right)^2-6\alpha \beta \delta ^2}}{2\delta}. \end{align*}$$ As before, |$a$| and |$b$| can be kept as real constants, provided that |$\alpha $| and |$\beta $| are chosen to have opposite signs. This completes the analysis of the second independent solution. All of the unknown constants, therefore, have been fully determined in Exton’s four-term recursion formula for |$G_n ,$| to the extent that the latter recursion formula establishes a strict identity with the four-term recursion relation of |$a_n ,$| defined in Eq. (2.6). As a consequence, this means that |$a_n$| can be precisely specified by |$G_n $| in Eq. (2.7), i.e., $$a_n =\frac{\left( a \right)_n \left( b \right)_n}{\left( c \right)_n n!}{ }_2F_2 \left(a+n,b+n;c+n,1+n;t\right)\!.$$ Following the same logic adopted in Sect. 2.1, one can insert the analytical form of |$a_n$| given above into the expansion of |$f\left( z \right)$|⁠, as well as collect together the results accompanied in the previous paragraphs of this section, to obtain an explicit form for the second independent solution to Mathieu’s equation. It follows then that the second independent solution, denoted here by |$\Psi _{II} (x),$| can be written in the analytical form |$(z=e^{ix})$|⁠, $$\begin{align}\label{eq2.10} \Psi_{II} (x)=z^{\tfrac{1}{2}+\delta }\exp \left( {-\frac{\alpha }{z}+\beta z} \right)\sum\limits_{n=-\infty }^\infty {(-1)^n\,\lambda ^n\frac{\left( a \right)_n \left( b \right)_n }{\left( c \right)_n }\frac{{ }_2F_2 \left( {a+n,b+n;c+n,1+n;t} \right)}{n!}} z^{-n}, \end{align}$$(2.10) which, like the first independent solution in Sect. 2.1, corresponds to being an exact solution to Mathieu’s equation, given by Eq. (1.1), or the quantum pendulum, defined by Eq. (1.5). Here, the corresponding parameters are given by $$\begin{align*} \alpha =\pm \frac{\omega}{\sqrt 2}, \quad \beta =\pm \frac{\omega}{2}, \quad t=-\lambda \beta,\quad \lambda =-\frac{2\delta }{3\beta},\quad \delta =\frac{10\pm 4\sqrt {9\varepsilon +4} }{3},\\ a=\frac{-R\pm \sqrt {R^2-S}}{2\delta},\quad b=\frac{-R\mp \sqrt {R^2-S}}{2\delta}, \quad c=\frac{\delta}{3}, \end{align*}$$ along with |$R=\delta -3\alpha \beta $| and |$S=6\alpha \beta \delta ^2$|⁠. 2.3. The Wronskian For abbreviation, let the solutions |$\Psi _I (x)$| and |$\Psi _{II} (x),$| as defined respectively in Eqs. (2.5) and (2.10), be represented as $$\Psi _I (x)=z^{-\tfrac{1}{2}}\exp \left(\! {-\frac{\alpha }{z}+\beta z}\!\right)\!\!\sum\limits_{n=-\infty }^\infty \, d_n \left( {r-z} \right)^{-n},\quad \Psi _{II} (x)=z^{\tfrac{1}{2}+\delta }\exp \left(\! {-\frac{{\alpha }'}{z}+{\beta }'z} \!\right) \!\!\sum\limits_{m=-\infty }^\infty \, e_m z^{-m},$$ whose coefficients, |$d_n, e_m$|⁠, have been defined as $$\begin{align*} d_n &=\left( {\frac{t}{s}} \right)^n{ }_3F_1 \left( {a,b,-n;c;-\tfrac{1}{t}} \right)\!,\\ e_m &=(-1)^m\,\lambda ^m\frac{\left( {{a}'} \right)_m \left( {{b}'} \right)_m }{\left( {{c}'} \right)_m }\frac{{ }_2F_2 \left( {{a}'+m,{b}'+m;{c}'+m,1+m;{t}'} \right)}{m!}. \end{align*}$$ Here the parameters appearing above for |$\Psi _I (x)$| and |$d_n $| are those given below Eq. (2.5), while the parameters appearing above for |$\Psi _{II} (x)$|and |$e_m $| correspond to those defined directly below Eq. (2.10), with the prime |$\prime $| being used to distinguish the similar labeling of the parameters adopted in |$\Psi _{II} (x)$| from |$\Psi _I (x).$| It is customary to verify the linear independence of |$\Psi _I (x)$| and |$\Psi _{II} (x)$| by computing their Wronskian, denoted by |$W\left( {\Psi _I ,\Psi _{II} } \right)\!,$| which is defined via the determinant: $$W \left(\Psi_I, \Psi_{II} \right) = \left|\begin{array}{cc} \Psi_I (x) & \Psi _{II} (x)\\ \frac{d\Psi _I (x)}{dx} & \frac{d\Psi _{II} (x)}{dx}\\ \end{array} \right|=\Psi _I (x)\frac{d\Psi _{II} (x)}{dx}-\Psi _{II} (x)\frac{d\Psi _I (x)}{dx}.$$ Employing the above representations for |$\Psi _I (x)$| and |$\Psi _{II} (x)$|⁠, and recalling that |$z=e^{ix},$| as well as noting that |$d / {dx} =izd / {dz},$| it follows that the Wronskian can be calculated as $$\begin{array}{l} W\left( {\Psi _I ,\Psi _{II} } \right)=iz^\delta \exp \left[ {-\frac{\alpha +{\alpha }'}{z}+\left( {\beta +{\beta }'} \right)z} \right]\sum\limits_{n=-\infty }^\infty \, \sum\limits_{m=-\infty }^\infty \, d_n e_m \left( {r-z} \right)^{-n}z^{-n}\left[ {\left( {n-m} \right)+\frac{nr}{r-z}} \right. \\ \text{ }\left. {-\frac{\left( {\alpha -{\alpha }'} \right)}{z}-\left( {\beta -{\beta }'} \right)z+\delta +1} \right]\!, \\ \end{array}$$ which is further simplified into the form $$\begin{align*} &W\left( {\Psi _I ,\Psi _{II} } \right)=iz^\delta \exp \left[\!{-\frac{\alpha +{\alpha }'}{z}+\left( {\beta +{\beta }'} \right)z}\! \right]\sum\limits_{n=-\infty }^\infty {\sum\limits_{m=-\infty }^\infty {d_n e_m \left( {r-z} \right)^{-n}z^{-m}} } \left[\! {\left. {\left( {n-m} \right)+\frac{nr}{r-z}}\!\right]} \right.\\ &\quad{} + i\,\left[ {\delta +1\left. {-\frac{\left( {\alpha -{\alpha }'} \right)}{z}-\left( {\beta -{\beta }'} \right)z} \right]} \right.\Psi _I (x)\Psi _{II} (x). \end{align*}$$ It is most readily seen from the second line of the above equation that the presence of the terms inside the square brackets will persist for any |$z$| value, showing clearly that |$\Psi _I (x)$| and |$\Psi _{II} (x)$| have a nonvanishing Wronskian, i.e., |$W\left( {\Psi _I ,\Psi _{II} } \right)\ne 0.$| Accordingly, |$\Psi _I (x)$| and |$\Psi _{II} (x)$| are linearly independent. 3. Discussion and conclusion As reviewed in Sect. 1, ever since Mathieu’s original work first made an appearance in 1868, the numerical and mathematical difficulties associated with solving Mathieu’s equation, as defined by Eq. (1.1), are tantamount to the intensive mathematical analysis and research applications that both physics and engineering demand [1,2]. At best, key findings and solutions to Mathieu’s equation tend to only exist in numerical tables as drawn from infinite continued fraction methods or approximate expressions for eigenvalue calculations, which at times are represented grossly by stability charts [2,6–8]. All the same, it was seen in Sect. 1.3 that solving Mathieu’s equation is also equivalent to solving the physical pendulum in quantum mechanics, i.e., quantum pendulum, so vested interests of Mathieu’s equation also exert themselves into future quantum technology applications [5,9]. Additionally, the concern over the characteristic and amplitude constants, |$a$| and |$q$|⁠, respectively, which are well known to be inherently and inextricably coupled in Mathieu’s equation, raised complex stability issues that are attached to their numerical computation [3,8]. Therefore, Sect. 2 was primarily motivated toward finding exact analytic solutions of Mathieu’s equation that will remove, or at least alleviate, any stability issues. Specifically, the present author showed in that section how two linearly independent solutions to Mathieu’s equation could be delineated using four-term recursion relations for hypergeometric functions to represent the expansion coefficients, as specified by |$a_n .$| In this sense, the use of four-term recursion relations has permitted a way of departing from the usual three-term recursion relation for the sole purpose of solving Mathieu’s equation analytically. Actually, the ability to devise two separate expansions techniques for consolidating the first and second linearly independent solutions, respectively, in Sects. 2.1 and 2.2 owes its success to underpinning two four-term recurrence relations from the very same research article by Exton [15], which in itself turned out to be a momentous coincidence. Finally, the general analytical solution |$\Psi (x)$| to Mathieu’s equation, for arbitrary constants |$c_1 $| and |$c_2 ,$| can be written as $$\Psi (x)=c_1 \Psi _I (x)+c_2 \Psi _{II} (x),$$ where |$\Psi _I (x)$| and |$\Psi _{II} (x)$| are the first and second linearly independent solutions, given respectively by Eqs. (2.5) and (2.10). Furthermore, the nature of the series expansions for |$\Psi _I (x)$| and |$\Psi _{II} (x)$| that were developed in Sect. 2 will remain as valid independent solutions for Mathieu’s equation (1.1) even when the |$x$| variable is relaxed from its restriction of an angular coordinate to that representing a spatial coordinate over the entire real line, should the need arise. This means that there is no dilemma over the validity of the general solution |$\Psi (x)$| above when it is non-periodic in |$x.$| Of course, when |$\Psi (x)$| is periodic in |$x,$| one will need to revert to the boundary conditions given in Sect. 1.2, namely, |$\Psi (x)=\Psi (x+\pi ),$| in order to calculate the eigenvalues, |$a.$| Overall, the problem of periodic and non-periodic boundary conditions of Mathieu’s equation is intended for subsequent research, where it is expected that any of the mathematical analysis associated with the boundary conditions will definitely benefit from being able to contract infinite sums over hypergeometric functions contained in the series expansions used for representing |$\Psi _I (x)$| and |$\Psi _{II} (x).$| This makes the latter aspect of the problem a pressing objective with future intentions over Mathieu’s equation. Acknowledgements The author is grateful to the referee for pointing out the reference to Heun’s differential equation. The present work has been prepared by the author in memory of the University of The Bahamas-North, which was devastated by the storm surge of Hurricane Dorian during the writing of this article. References [1] Mathieu É. , J. Math. Pures Appl. 13 , 137 ( 1868 ). ( http://eudml.org/doc/234720 ) [2] McLachlan N. W. , Theory and Application of Mathieu Functions ( Oxford University Press , London , 1947 ). Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC [3] Ruby L. , Am. J. Phys. 64 , 39 ( 1996 ); 64 , 655 ( 1996 ) [erratum]. ( https://doi.org/10.1119/1.18290; https://doi.org/10.1119/1.18476 ) Crossref Search ADS [4] Paul W. , Rev. Mod. Phys. 62 , 531 ( 1990 ). ( http://dx.doi.org/10.1103/RevModPhys.62.531 ) Crossref Search ADS [5] Cirac J. I. and Zoller P.,Phys. Rev. Lett. 74 , 4091 ( 1995 ). ( http://dx.doi.org/10.1103/PhysRevLett.74.4091 ) Crossref Search ADS [6] Morse M. and Feshbach H., Methods of Theoretical Physics ( McGraw-Hill , New York , 1953 ), Part 1, Sect. 5.2. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC [7] Abramowitz M. and Stegun I. A., Handbook of Mathematical Functions ( US Government Printing Office , Washington, DC , 1964 ), Chap. 20. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC [8] Kovacic I. , Rand R., and Sah S. M., AMSE Appl. Mech. Rev. 70 , 020802 ( 2018 ). ( http://dx.doi.org/10.1115/1.4039144 ) Crossref Search ADS [9] Galvez E. J. , Auccapuclla F. J., Wittler K. L., and Qin Y., Proc. SPIE Complex Light Opt. Forces XIII 10935 , 1093509 ( 2019 ). ( https://doi.org/10.1117/12.2510822 ) [10] Condon E. U. ,Phys. Rev. 31 , 891 ( 1928 ). ( http://dx.doi.org/10.1103/PhysRev.31.891 ) [11] Ciftci H. , Hall R. L., and Saad N., Cent. Eur. J. Phys. 11 , 37 ( 2013 ). ( http://dx.doi.org/10.2478/s11534-012-0147-3 ) [12] Pöschl G. and Teller E., Z. Phys. 83 , 143 ( 1933 ). ( https://doi.org/10.1007/BF01331132 ) Crossref Search ADS [13] Ronveaux A. (ed.), Heun’s Differential Equations ( Oxford University Press , New York , 1995 ). [14] El-Jaick L. J. and Figueiredo B. D. B., J. Math. Phys. 49 , 083508 ( 2008 ). ( https://doi.org/10.1063/1.2970150 ) Crossref Search ADS [15] Exton H. , Collect. Math. 49 , 43 ( 1998 ). ( https://www.raco.cat/index.php/CollectaneaMathematica/article/view/56439 ) © The Author(s) 2020. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - Exact solutions of Mathieu’s equation
JO - Progress of Theoretical and Experimental Physics
DO - 10.1093/ptep/ptaa024
DA - 2020-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/exact-solutions-of-mathieu-s-equation-LZLSVmn61E
VL - 2020
IS - 4
DP - DeepDyve
ER -