Harmonic Oscillator, WKB/Airy Approximations

The WKB approximation estimates solutions of the one-dimensional Schrödinger equation \displaystyle -\frac{\hbar^2}{2m}\psi''(x)+V(x)\psi(x)=E\psi(x) when V(x) changes slowly over a wavelength. The idea is that the quantum wave is locally like a free-particle wave with slowly changing momentum p(x)=\sqrt{2m[E-V(x)]}.

In a classically allowed region, E>V(x) , WKB gives

\displaystyle \psi(x)\approx\frac{1}{\sqrt{p(x)}}\cos\left[\frac1\hbar\int^x p(t)dt+\varphi\right].

Thus the local wavelength is h/p(x) , and the amplitude grows where the particle moves slowly. In a forbidden region, E<V(x) , define \kappa(x)=\sqrt{2m[V(x)-E]} ; then the wave is approximately exponential,

\displaystyle \psi(x)\approx\frac{1}{\sqrt{\kappa(x)}}\exp\left[-\frac1\hbar\int\kappa(t)dt\right].

A turning point x_t satisfies V(x_t)=E . Classically, the particle stops and reverses there; quantum mechanically, the wave changes from oscillatory to decaying. Raw WKB fails exactly at x_t , because p(x_t)=0 makes its prefactor diverge. Near a simple turning point, the correct local approximation is an Airy function: it stays finite at the turning point, oscillates on the allowed side, and decays on the forbidden side. Thus WKB is accurate away from turning points, Airy is accurate near them, and in their overlap region the Airy asymptotics reproduce WKB—including the \pi/4 phase shift.

The harmonic oscillator is one of the best places to understand WKB approximation because we know the exact answer. Usually, in quantum mechanics, one writes down an approximation because the exact wavefunction is unknown. Here the situation is better: the exact wavefunctions are Hermite functions, and we can ask a more revealing question: When the quantum number n is large, what does the exact Hermite wavefunction actually look like?

The answer is that it naturally breaks into three regions. In the middle, where the classical particle is allowed to move, the exact wavefunction becomes an oscillating WKB wave. Near the places where the classical particle turns around, ordinary WKB fails and the exact wavefunction becomes Airy-like. Beyond those turning points, the same exact wavefunction becomes an exponentially decaying WKB tail. Thus WKB and Airy functions are different asymptotic faces of one exact solution.

The purpose of this discussion is to show that analytically. We will begin with the exact Hermite wavefunction, convert it into an exact contour integral, and study that integral when n is large. WKB, the Airy function, the classically allowed region, the forbidden region, and the phase shift \pi/4 will all emerge from the geometry of that integral.

Harmonic Oscillator

The one-dimensional harmonic oscillator has potential \displaystyle V(x)=\frac12 m\omega^2x^2.
Its time-independent Schrödinger equation is

\displaystyle -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} +\frac12m\omega^2x^2\psi = E\psi.

It is useful to measure distance in units of the oscillator length \displaystyle \ell_0=\sqrt{\frac{\hbar}{m\omega}}. Define the dimensionless coordinate \displaystyle x\longmapsto \xi=\frac{x}{\ell_0} = \sqrt{\frac{m\omega}{\hbar}}\,x. To avoid heavy notation, write \xi simply as x .

Then the exact normalized n -th eigenfunction is \displaystyle \phi_n(x) = \frac{1}{\pi^{1/4}\sqrt{2^n n!}} H_n(x)e^{-x^2/2}, where H_n(x) is the n -th Hermite polynomial. The corresponding energy is E_n=\left(n+\frac12\right)\hbar\omega. It is convenient to introduce \displaystyle N=n+\frac12, \quad A=\sqrt{2N}=\sqrt{2n+1}.
Then the exact differential equation for \phi_n can be written as

\displaystyle \phi_n''(x)+(A^2-x^2)\phi_n(x)=0.


The classical turning points occur where the kinetic energy vanishes. Since the dimensionless kinetic-energy factor is \displaystyle A^2-x^2, the turning points are x=-A,x=A. Inside this interval the classical momentum is real. This is the classically allowed region. Outside, \displaystyle |x|>A, we have A^2-x^2<0 , so classical motion is impossible. This is the classically forbidden region. For later use, define k(x)=\sqrt{A^2-x^2},  |x|<A, and \kappa(x)=\sqrt{x^2-A^2},  |x|>A. In physical units, k(x) is the local momentum divided by \hbar .

WKB prediction

Before deriving anything from Hermite functions, let us state what WKB would predict. Inside the allowed region, WKB predicts an oscillatory wave:

\displaystyle \psi_{\mathrm{WKB}}(x) \approx \frac{\text{constant}}{\sqrt{k(x)}} \cos\left( \int^x k(t)dt+\text{phase} \right).

Two important features appear here. First, the local wavelength is approximately \displaystyle \lambda(x)\sim \frac{2\pi}{k(x)}. Thus as x approaches a turning point and k(x)\to0 , the wavelength becomes longer and longer. Second, the amplitude has the form \displaystyle \frac{1}{\sqrt{k(x)}}. This reflects the fact that a classical particle moves more slowly near a turning point, and therefore spends more time there.

Outside the allowed region, WKB predicts an exponential tail:

\displaystyle \psi_{\mathrm{WKB}}(x) \approx \frac{\text{constant}}{\sqrt{\kappa(x)}} \exp\left( -\int\kappa(t)dt \right).

However, ordinary WKB has a problem at a turning point. Since \displaystyle k(A)=0,
the expression \displaystyle \frac{1}{\sqrt{k(x)}} diverges as x\to A . The true wavefunction does not diverge. The purpose of the Airy function is to repair this failure. The main point of what follows is that these formulas are not merely plausible guesses. For the harmonic oscillator, they can be extracted from the exact Hermite function.

Analysis of exact solutions

The Hermite polynomials have the exact generating function \displaystyle e^{2xz-z^2} = \sum_{m=0}^{\infty} \frac{H_m(x)}{m!}z^m. The coefficient of z^n can be recovered by a contour integral: \displaystyle H_n(x) = \frac{n!}{2\pi i} \oint_\Gamma \frac{e^{2xz-z^2}}{z^{n+1}}\,dz. Here \Gamma is any positively oriented closed contour surrounding z=0 . Substituting this exact formula into the normalized wavefunction gives

\displaystyle \phi_n(x) = \frac{\sqrt{n!}}{\pi^{1/4}2^{n/2}} e^{-x^2/2} \frac{1}{2\pi i} \oint_\Gamma e^{2xz-z^2}z^{-n-1}\,dz.


Nothing has been approximated yet. Now rescale position by the turning-point size: \displaystyle x=Ay, \quad A=\sqrt{2N}. Thus y=x/A measures position relative to the classical amplitude. The right turning point is now at y=1 , and the left turning point is at y=-1 . Also rescale the contour variable by \displaystyle z=\sqrt{\frac N2}\,w.
After straightforward algebra, the exact wavefunction becomes

\displaystyle \phi_n(Ay) = C_N \frac{1}{2\pi i} \oint w^{-1/2}e^{N G_y(w)}\,dw,

where \displaystyle G_y(w) = 2yw-\frac12w^2-\log w-y^2, and \displaystyle C_N = \pi^{-1/4}\sqrt{n!}\,N^{-n/2}.

The expression above is exact. The quantity N=n+\tfrac12 is large when the oscillator is highly excited. Thus the integral has the form \displaystyle \int e^{N G(w)}dw. Integrals of this kind are governed, for large N , by the stationary points of G . This is the method of steepest descent.

Saddle Point Approximation

Suppose one has an integral such as \displaystyle I_N=\int e^{Nf(z)}g(z)dz, where N is very large. Because of the factor e^{Nf(z)} , the integral is strongly concentrated near points where the phase changes least rapidly. These are the points where \displaystyle f'(z_0)=0. Such points are called saddle points or stationary points. If \displaystyle f''(z_0)\neq0, then near z_0 , \displaystyle f(z) = f(z_0) + \frac12f''(z_0)(z-z_0)^2+\cdots. Thus the integral near the saddle becomes approximately Gaussian. The Gaussian approximation produces the familiar N^{-1/2} factor and, in complex problems, a phase coming from the square root of a complex number. If two saddles merge together, then the quadratic term vanishes: \displaystyle f'(z_0)=f''(z_0)=0. The first nonzero term is often cubic: f(z)\approx f(z_0)+c(z-z_0)^3. A cubic saddle is governed by the Airy function. This is why Airy functions appear at turning points.

Thus the entire WKB story will come from following the saddle points of G_y(w) as y moves from inside to outside the interval [-1,1] . Differentiate the exact exponent: \displaystyle G_y'(w) = 2y-w-\frac1w. The saddle points satisfy \displaystyle 2y-w-\frac1w=0. Multiplying by w , \displaystyle w^2-2yw+1=0. Therefore \displaystyle w_\pm = y\pm\sqrt{y^2-1}.
This formula contains the whole geometry.

When |y|<1 , the quantity y^2-1 is negative, so the saddles are complex conjugates. Their contributions interfere, producing oscillations. When y>1 , the saddles are real. One produces a decaying exponential and the other corresponds to a growing exponential. When y=1 , the two saddles collide. This is the turning point, and this collision produces the Airy function.

Classical Region

Suppose |y|<1. Write \displaystyle y=\cos\theta,  0<\theta<\pi. Then \sqrt{y^2-1}=i\sin\theta, and the two saddle points become w_+=e^{i\theta}, w_-=e^{-i\theta}. The saddles are complex conjugates. This matters because the two contributions have equal magnitude but opposite phase. Adding them produces a cosine. At the first saddle, G_y(e^{i\theta})=\frac12-i\left( \theta-\frac12\sin2\theta \right). At the conjugate saddle, G_y(e^{-i\theta}) = \frac12 + i\left( \theta-\frac12\sin2\theta \right). Define S(\theta) = N\left( \theta-\frac12\sin2\theta \right). The two exponential factors are then e^{-iS(\theta)} and e^{+iS(\theta)}. These are the two semiclassical traveling waves.

To estimate the integral near each saddle, one expands G_y(w) to second order. For example, \displaystyle G_y(w) = G_y(w_+) + \frac12G_y''(w_+)(w-w_+)^2+\cdots. The relevant second derivatives are \displaystyle G_y''(e^{i\theta}) = -2ie^{-i\theta}\sin\theta, and \displaystyle G_y''(e^{-i\theta}) = 2ie^{i\theta}\sin\theta. The local integrals are Gaussian, but they are complex Gaussians. The direction of steepest descent and the square root of the complex curvature contribute a quarter-turn phase. After carrying out the standard saddle-point calculation, the two contributions combine into \displaystyle \cos\left( S(\theta)-\frac{\pi}{4} \right). The \pi/4 is therefore phase produced by the local Gaussian geometry of the exact contour integral.

After including the amplitude factors and Stirling’s approximation for n! , one obtains the asymptotic formula

\displaystyle \phi_n(A\cos\theta) = \sqrt{\frac{2}{\pi A\sin\theta}} \left[ \cos\left( N\left(\theta-\frac12\sin2\theta\right) -\frac{\pi}{4} \right) + O_\delta(N^{-1}) \right].


This is valid uniformly when \displaystyle \delta\leq\theta\leq\pi-\delta. It is valid everywhere in the allowed region except in small neighborhoods of the turning points.

The exact saddle-point formula already looks like WKB, but let us show that it is precisely WKB.

The local wave number is k(x)=\sqrt{A^2-x^2}. Set x=A\cos\theta. Then k(x)=A\sin\theta. Thus the amplitude factor becomes \sqrt{\frac{2}{\pi A\sin\theta}} = \sqrt{\frac{2}{\pi}} \frac{1}{\sqrt{k(x)}}. This is exactly the WKB amplitude. Now examine the phase. The WKB action from x to the right turning point A is \mathcal S(x) = \int_x^A \sqrt{A^2-t^2}dt. Put t=A\cos u . Then dt=-A\sin u\,du, and \sqrt{A^2-t^2}=A\sin u. Therefore \mathcal S(x) = A^2 \int_0^\theta \sin^2u du = \frac{A^2}{2} \left( \theta-\frac12\sin2\theta \right).
But A^2=2N. Hence \mathcal S(x) = N\left( \theta-\frac12\sin2\theta \right). This is exactly the phase that emerged from the saddle points. Thus the exact Hermite asymptotic becomes

\displaystyle \phi_n(x) = \sqrt{\frac{2}{\pi}} \frac{1}{\sqrt{k(x)}} \left[ \cos\left( \int_x^A k(t)dt-\frac{\pi}{4} \right) + O(N^{-1}) \right].


This is exactly the oscillatory WKB wavefunction.

Forbidden Region

Now suppose y>1. Write y=\cosh\eta,  \eta>0. Then the two saddle points are real: w_+=e^\eta, w_-=e^{-\eta}. The exact normalized harmonic oscillator wavefunction must decay as x\to+\infty . Therefore its asymptotics must select the recessive, or decaying, saddle. A careful steepest-descent deformation of the original coefficient contour passes through the saddle w=e^{-\eta} . The other saddle corresponds to the exponentially growing solution, which is excluded by normalizability.

At the relevant saddle, G_y(e^{-\eta})=\frac12-\left( \sinh\eta\cosh\eta-\eta \right). Thus the exponential part of the wavefunction is \exp\left[ -N\left( \sinh\eta\cosh\eta-\eta \right) \right]. The full saddle-point calculation gives

\displaystyle \phi_n(A\cosh\eta) = \frac{1}{\sqrt{2\pi A\sinh\eta}} \exp\left[ -N\left( \sinh\eta\cosh\eta-\eta \right) \right] \left[ 1+O_\delta(N^{-1}) \right].


Again, this formula is valid away from the turning point.

Outside the turning point, \kappa(x)=\sqrt{x^2-A^2}. If x=A\cosh\eta, then \kappa(x)=A\sinh\eta. Thus the amplitude is \frac{1}{\sqrt{2\pi A\sinh\eta}} = \frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{\kappa(x)}}. Now compute the forbidden WKB action: \mathcal K(x) = \int_A^x\sqrt{t^2-A^2}dt. Put t=A\cosh u. Then dt=A\sinh u du, and \sqrt{t^2-A^2}=A\sinh u. Therefore \mathcal K(x) = A^2\int_0^\eta\sinh^2u du. Using \int_0^\eta\sinh^2u du = \frac12\sinh\eta\cosh\eta-\frac12\eta, we get \mathcal K(x) = N\left( \sinh\eta\cosh\eta-\eta \right). Hence the exact Hermite asymptotic becomes

\displaystyle \phi_n(x) = \frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{\kappa(x)}} \exp\left[ -\int_A^x\kappa(t)dt \right] \left[ 1+O(N^{-1}) \right].


This is the decaying WKB solution. Thus, outside the classical region, the exact Hermite function becomes the WKB exponential tail.

Turning Point

The WKB formulas we have derived are very good away from x=A . But both contain singular amplitudes: \frac{1}{\sqrt{k(x)}} \quad\text{or}\quad \frac{1}{\sqrt{\kappa(x)}}. At the turning point, k(A)=0, \kappa(A)=0. Thus raw WKB diverges.

The exact Hermite function does not diverge. It remains smooth and finite. The contour-integral explanation is equally clear. At y=1 , the two saddle points become w_+=w_-=1. The two saddles merge. Ordinary saddle-point approximation assumes that the second derivative is nonzero. But at the merged saddle, G_1'(1)=0, G_1''(1)=0. The Gaussian approximation fails. This is exactly analogous to WKB failing because k=0 . We need a new local approximation.

Near the right turning point, write \displaystyle x=A+\frac{s}{\alpha}, where \displaystyle \alpha=(2A)^{1/3}. Since \displaystyle A\sim n^{1/2}, we have \displaystyle \alpha\sim n^{1/6}. Therefore the turning-point layer has width \displaystyle |x-A|\sim \alpha^{-1}\sim n^{-1/6}. This is much smaller than the total oscillator size \displaystyle A\sim n^{1/2}. For highly excited states, the oscillator therefore has a large bulk oscillatory region, a narrow Airy transition layer, and then a forbidden exponential tail. The special n^{-1/6} scale is not arbitrary. It comes from balancing the cubic saddle behavior against the small displacement from the turning point.

Near the turning point, write \displaystyle y=1+\frac{s}{2N^{2/3}}, and zoom into the merged saddle using \displaystyle w=1-N^{-1/3}t. The exponent is \displaystyle G_y(w)=2yw-\frac12w^2-\log w-y^2. At y=1 , its Taylor expansion around w=1 is \displaystyle G_1(1+u) = \frac12 -\frac{u^3}{3} +\frac{u^4}{4} +O(u^5).
Notice what has happened: the linear and quadratic terms vanish. The first nontrivial term is cubic. Substituting \displaystyle u=-N^{-1/3}t, and including the small displacement of y away from 1 , one finds N G_y(w) = \frac N2 + \frac{t^3}{3} -st + O\left( N^{-1/3}(1+|t|^4+|s|^2) \right). Therefore the leading local integral has the form

\displaystyle \frac{1}{2\pi i} \int_{\mathcal C} \exp\left( \frac{t^3}{3}-st \right)dt.

But this is exactly the contour integral defining the Airy function:

\displaystyle {\text{Ai}}(s) = \frac{1}{2\pi i} \int_{\mathcal C} \exp\left( \frac{t^3}{3}-st \right)dt.

After restoring the prefactor, we obtain

\displaystyle \phi_n\left(A+\frac{s}{\alpha}\right) = \sqrt{\frac{2}{\alpha}} \left[ {\text{Ai}}(s) + O_S(N^{-1/3}) \right].

This is valid when s remains in a fixed bounded interval. Thus the Airy function is the correct local approximation to the exact Hermite wavefunction near the turning point.

The Airy variable is \displaystyle s=\alpha(x-A). Therefore: \displaystyle s<0 \Longleftrightarrow\quad x<A, which is the allowed side. \displaystyle s=0 \Longleftrightarrow\quad x=A, which is the turning point. \displaystyle s>0 \Longleftrightarrow\quad x>A,
which is the forbidden side. The single Airy function {\text{Ai}}(s) does all three things: {\text{Ai}}(s<0 oscillates, {\text{Ai}}(0), is finite and {Ai}(s>0 decays exponentially. This is why it is the right local function.

Raw WKB gives the correct behavior away from the turning point, but it fails precisely where the character of the wave changes. The Airy function is the smooth bridge through that change.

For large positive r , {\text{Ai}}(-r) = \frac{1}{\sqrt{\pi}} r^{-1/4} \left[ \sin\left( \frac23r^{3/2}+\frac{\pi}{4} \right) + O(r^{-3/2}) \right]. Take a point slightly inside the turning point: x=A-d, \quad d>0. Then s=-\alpha d=-r. Near the turning point, k(x)=\sqrt{A^2-(A-d)^2}=\sqrt{2Ad-d^2} \approx \sqrt{2Ad}. Since r=\alpha d, \qquad \alpha^3=2A, we obtain \frac23r^{3/2} = \frac23\sqrt{2A}\,d^{3/2}. But this is the leading approximation to the WKB action: \int_{A-d}^{A}k(t)dt = \frac23\sqrt{2A}\,d^{3/2} \left[ 1+O\left(\frac dA\right) \right]. The amplitude also matches: \sqrt{\frac{2}{\alpha}} \frac{1}{\sqrt{\pi}}r^{-1/4} = \sqrt{\frac{2}{\pi}} \frac{1}{\sqrt{k(x)}} \left[ 1+O\left(\frac dA\right) \right]. Therefore the Airy formula becomes

\displaystyle \phi_n(x) \sim \sqrt{\frac{2}{\pi}} \frac{1}{\sqrt{k(x)}} \sin\left[ \int_x^A k(t)dt+\frac{\pi}{4} \right].

Since \sin\left(S+\frac{\pi}{4}\right) = \cos\left(S-\frac{\pi}{4}\right), this is exactly the WKB formula obtained earlier from the bulk saddle analysis.

For large positive r , {\text{Ai}}(r) = \frac{1}{2\sqrt{\pi}} r^{-1/4} \exp\left( -\frac23r^{3/2} \right) \left[ 1+O(r^{-3/2}) \right]. Now take a point slightly outside the turning point: x=A+d, \quad d>0. Then \displaystyle s=\alpha d=r. Near the turning point, \kappa(x)=\sqrt{(A+d)^2-A^2}=\sqrt{2Ad+d^2} \approx \sqrt{2Ad}. The WKB exponent is \int_A^{A+d}\kappa(t)dt = \frac23\sqrt{2A}\,d^{3/2} \left[ 1+O\left(\frac dA\right) \right]. But \frac23r^{3/2} = \frac23\sqrt{2A}\,d^{3/2}. The Airy prefactor also becomes the WKB prefactor: \sqrt{\frac{2}{\alpha}} \frac{1}{2\sqrt{\pi}}r^{-1/4} = \frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{\kappa(x)}} \left[ 1+O\left(\frac dA\right) \right]. Therefore

\displaystyle \phi_n(x) \sim \frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{\kappa(x)}} \exp\left[ -\int_A^x\kappa(t)dt \right].


This is precisely the decaying WKB tail.

The Airy approximation is valid near the turning point, while raw WKB is valid away from it. There is a region where both are valid. The Airy asymptotic formulas require \displaystyle |s|=\alpha|x-A|\gg1. Thus \displaystyle |x-A|\gg \alpha^{-1}. But the local turning-point expansion also requires \displaystyle |x-A|\ll A. So the overlap region is \displaystyle \alpha^{-1}\ll |x-A|\ll A.
In that region, Airy asymptotics and WKB give the same formula. This is the precise meaning of “matching.” They are not the same formula at the turning point. Raw WKB diverges there, while Airy stays finite. But there is a region sufficiently far from the turning point for Airy to simplify, and sufficiently near it for the local geometry to remain relevant. In that overlap region, the two approximations agree.

Quantization

The harmonic oscillator has two turning points, at -A and A . Each smooth turning point contributes a phase shift of \pi/4 . The total action across the allowed region is \int_{-A}^{A}\sqrt{A^2-x^2}\,dx. This is the area of a semicircle of radius A : \int_{-A}^{A}\sqrt{A^2-x^2}\,dx = \frac{\pi A^2}{2}. Since A^2=2N , \int_{-A}^{A}\sqrt{A^2-x^2}\,dx=\pi N=\pi\left(n+\frac12\right). Thus the WKB quantization rule is

\displaystyle \int_{x_-}^{x_+}k(x)\,dx = \left(n+\frac12\right)\pi.

For the harmonic oscillator, this semiclassical quantization formula gives the exact energies. The extra 1/2 is not an arbitrary correction. It comes from the two turning points: \displaystyle \frac{\pi}{4} + \frac{\pi}{4} = \frac{\pi}{2}. That total phase shift changes the naive integer rule into the correct half-integer rule.

Conclusion:

The three formulas describe one exact Hermite wavefunction in different regions. Inside |x|<A , two complex saddle points contribute and interfere, giving the oscillatory WKB wave:

\displaystyle \phi_n(x)\approx\sqrt{\frac{2}{\pi}}\frac{1}{\sqrt{k(x)}}\cos\left[\int_x^A k(t)\,dt-\frac{\pi}{4}\right].

Near x=A , the saddles merge, so the usual Gaussian saddle approximation fails. The resulting cubic approximation is the Airy function:

\displaystyle \phi_n\left(A+\frac{s}{(2A)^{1/3}}\right)\approx\sqrt{\frac{2}{(2A)^{1/3}}}{\text{Ai}}(s).

Outside x>A , one real saddle gives the decaying WKB tail:

\displaystyle \phi_n(x)\approx\frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{\kappa(x)}}\exp\left[-\int_A^x\kappa(t)\,dt\right].

So for the exact Hermite wavefunction, WKB is the correct large-n approximation away from the turning points, while the Airy function is the correct approximation near them.

Leave a comment