Asymptotic Behavior of Power Series Terms in the Hermite Equation

Question

When solving for the wave function under a harmonic potential $V(x)=\frac{1}{2}kx^2$, we attempt a solution in the form:

$$\psi(y)=H(y)e^{-\frac{1}{2}y^2},\quad y:=\left(\frac{m\omega_o}{\hbar}\right)^{\frac{1}{2}}x,\quad \alpha:=\frac{2E}{\hbar \omega_o},$$

where $x$ is position, $E$ is the energy. It follows by substitution into the time-independent Schrödinger equation that $H(y)$ is subject to the following:

$$H''-2yH'+(\alpha-1)H=0,$$

which is the Hermite equation. When solving this equation using Frobenius' method, we end up with the recursive relation:

$$a_{n+2}=\frac{2n+1-\alpha}{(n+1)(n+2)}a_n \quad\text{where}\quad H(y)=\sum_{n=0}^{\infty}a_ny^n.$$

When proving that this series must terminate at some finite order $n$ for the final wave function to be square-integrable, it is common to compare the limiting behavior of the ratio between successive terms in this series with that of the function $e^{2y^2}$:

$$H_\infty:\quad \frac{a_{n+2}}{a_n}\sim\frac{2}{n},$$ $$e^{2y^2}:\quad \frac{b_{n+2}}{b_n}\sim\frac{2}{n},$$

since the power series expansion for the latter is:

$$e^{2y^2}=\sum_{n=0}^\infty b_ny^n=\sum_{n=0}^\infty \frac{2^ny^{2n}}{n!}.$$

From here we assume that:

$$\frac{a_{n+2}}{a_n}\sim\frac{b_{n+2}}{b_n}\implies H_\infty(y)\sim e^{2y^2},$$

and therefore the wave function diverges as $y\to \infty$.

My question is: how do we prove that, if the ratios between successive terms in two different power series are asymptotically equivalent, then the two function defined by the power series are also asymptotically equivalent?

Answer 1

From @Gary's very useful comment:

I believe you can express the series solution in terms of the confluent hypergeometric function $_1F_1(−a,b,y^2) \ldots $

and after some research, the method I presented is not rigorous and not true in all cases. A more sensible method is as follows.

From the general physicist's Hermite differential equation:

$$\frac{\text{d}^2y}{\text{d}x^2}-2x\frac{\text{d}y}{\text{d}x}+2\lambda y=0,\quad \lambda\in\mathbb R,$$

we can use Frobenius' method to derive a series solution

$$y(x)=\sum_{k=0}^\infty a_kx^k,\quad a_{k+2}=\frac{2(k-\lambda)}{(k+1)(k+2)},$$

which implies that the solution is a linear combination of an even series $y_1(x)$ and an odd series $y_2(x)$:

$$y_1(x):=\sum_{k=0}^\infty a_{2k}\,x^{2k},\quad y_2(x):=\sum_{k=0}^\infty a_{2k+1}\,x^{2k+1},\quad \text{such that}\quad y(x)=y_1(x;a_0)+y_2(x;a_1),$$

where $a_0$ and $a_1$ are arbitrary coefficients for the above recursive relationship. From this relationship, we see that one of the series solutions ($y_1$ or $y_2$ exclusively) terminates for positive integer values of $\lambda$:

$$a_{k+2}=\frac{2(k-\lambda)}{(k+1)(k+2)}=0\implies k=\lambda.$$

---

Aside, we can define two functions:

$$H_\lambda(x):=\begin{cases} \sum_{k=0}^\infty a_{2k}\,x^{2k} & \text{if }\lambda \text{ is even} \\ \sum_{k=0}^\infty a_{2k+1}\,x^{2k+1} & \text{if }\lambda \text{ is odd} \\ 0 & \text{otherwise} \end{cases}\quad \text{where}\quad a_{k+2}=\frac{2(k-\lambda)}{(k+1)(k+2)},$$

$$h_\lambda(x):=\begin{cases} \sum_{k=0}^\infty a_{2k}\,x^{2k} & \text{if }\lambda \text{ is odd} \\ \sum_{k=0}^\infty a_{2k+1}\,x^{2k+1} & \text{if }\lambda \text{ is even} \\ \sum_{k=0}^\infty a_{k}\,x^{k} & \text{otherwise} \end{cases}\quad \text{where}\quad a_{k+2}=\frac{2(k-\lambda)}{(k+1)(k+2)};$$

where we call $H_\lambda$ the physicist's Hermite polynomial of first kind, and $h_\lambda$ the physicist's Hermite function of second kind. We can also see that the latter can be expressed as a confluent hypergeometric function of first kind:

$$h_\lambda(x)={}_1F_1(-\frac{\lambda}{2};\frac{1}{2};x^2).$$

Lastly, for $H_\lambda$ by convention we specify the initial coefficients $a_0=1$ and $a_1=2$.

---

Now, instead of partitioning the solution based of parity, we split the solution in two functions: a finite-order polynomial, and an infinite power series, such that the solution is a linear combination of the two,

$$y(x)=C_1H_\lambda(x)+C_2h_\lambda(x),\quad C_1,C_2\in\mathbb C.$$

In conclusion: applying this to the problem I originally asked in the question, we require $y\to 0$ as $|x|\to\infty$, as to render the final wavefunction square-integrable (or at least have $y$ grow slower than the rate $e^{-\frac{1}{2}x^2}$ decreases). This can only hold if $C_2=0$, since $h_\lambda(x)$ diverges as $|x|\to\infty$. And if this is the case, we are left with:

$$y(x)=C_1 H_n(x),\quad n\in\mathbb N_0,$$

where we changed $\lambda$ to $n$ to further signify that this parameter is strictly a positive integer. Finally, in the time-independent Schrödinger equation in the original question, we specified $n=\lambda=\frac{1}{2}(\alpha-1)$; applying this change of variables, and substituting into the original ansatz for the differential equation, we arrive at the final solution.

Please edit this answer if something is not quite right.

Answer 2

This material can be found in Bender and Orszag, Chapter 3. (They use the probabilist's convention, though, and derive $He(x)$ rather than $H(x)$.)

---

There is an alternative way to go about analyzing this problem in which we compute an asymptotic series about the irregular point at infinity. This first involves extracting the leading order behavior of the solutions to the original (scaled) differential equation, given by $$ -\frac{1}{2}y''+\frac{1}{2}x^2y-\epsilon y=0\,, $$ by first making the transformation $y(x)=e^{S(x)}$ and then using the method of dominant balance. This method is systematic if not entirely rigorous, and the result is $$ y(x)\sim e^{-x^2/2}x^{\nu}\,, $$ where $\nu=\epsilon-1/2$. At this point, we guess that the solution can be written as $$ y(x) = e^{-x^2/2}x^{\nu}\sum_{n=0}^{\infty}a_nx^{-n}\,. $$ (Note that this is a power series in negative powers of $x$!) This results in a recurrence relation for the $a_n$'s, given by $$ a_1=0,~~~~~ a_{n+2} = -\frac{(n-\nu)(n-\nu+1)}{2(n+2)}a_n\,, $$ so that the solutions looks like $$ y(x) = e^{-x^2/2}x^{\nu} \left( 1-\frac{\nu(\nu-1)}{4x^2} + \frac{\nu(\nu-1)(\nu-2)(\nu-3)}{32x^4} - \frac{\nu(\nu-1)(\nu-2)(\nu-3)(\nu-4)(\nu-5)}{384x^6} + \cdots \right) $$ Crucially, we can compute the radius of convergence of this series, and it is zero. It is still a good asymptotic solution for any $\nu$, but it is not good for the situation we are concerned with, which is for the function to go to zero at infinity. This computation of the radius of convergence of the asymptotic series replaces the hand-wavy "the power series acts like $e^{x^2}$" with a slightly less hand-wavy "the radius of convergence of the series expansion about $x\to\infty$ is zero".

In any case, this means that we only get a "good" solution when the power series truncates to a finite number of terms (because in that case the radius of convergence is of course infinite!), and this occurs when $\nu$ is a non-negative integer, which can be clearly seen from both the recursion relation and the first few terms of the expansion. After multiplying through by the factor of $x^{\nu}$, we get exactly the physicist Hermite polynomials $H_n(x)$, as we should.

---

The problem is that it's not clear to what extent the radius of convergence being zero is really related to the function blowing up at infinity, so this doesn't entirely answer the question, unfortunately.