## WKB approximation and the hydrogen atom

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.14.

As another example of applying the WKB approximation to a 3-d problem with a spherically symmetric potential, we’ll look at the radial equation for the hydrogen atom.

The general wave function can be written as the product of a radial function ${R\left(r\right)}$ and a spherical harmonic ${Y\left(\theta,\phi\right)}$.

With the substitution ${u\left(r\right)\equiv rR\left(r\right)}$ the radial equation for hydrogen can be written

$\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}u}{dr^{2}}+\left(-\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{r}+\frac{\hbar^{2}}{2m}\frac{l(l+1)}{r^{2}}\right)u=Eu \ \ \ \ \ (1)$

The effective potential is the term in parentheses:

$\displaystyle V\left(r\right)=-\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{r}+\frac{\hbar^{2}}{2m}\frac{l(l+1)}{r^{2}} \ \ \ \ \ (2)$

and has the form of a well with sloping sides on both sides, at points ${r_{1}}$ and ${r_{2}}$ defined by the roots of

$\displaystyle E=-\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{r}+\frac{\hbar^{2}}{2m}\frac{l(l+1)}{r^{2}} \ \ \ \ \ (3)$

In this case, the WKB equations satisfy the condition

$\displaystyle \int_{r_{1}}^{r_{2}}p\left(r\right)\; dr=\left(n-\frac{1}{2}\right)\pi\hbar \ \ \ \ \ (4)$

Plugging in the formulas for the hydrogen atom we get

 $\displaystyle \int_{r_{1}}^{r_{2}}p\left(r\right)\; dr$ $\displaystyle =$ $\displaystyle \sqrt{2m}\int_{r_{1}}^{r_{2}}\sqrt{E+\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{r}-\frac{\hbar^{2}}{2m}\frac{l(l+1)}{r^{2}}}dr\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mE}\int_{r_{1}}^{r_{2}}\frac{1}{r}\sqrt{-r^{2}-\frac{e^{2}}{4\pi\epsilon_{0}E}r+\frac{\hbar^{2}}{2mE}l\left(l+1\right)}dr \ \ \ \ \ (6)$

where we’ve factored out ${-E}$ (since for a bound state, ${E<0}$). We can simplify the notation by defining the positive quantities

 $\displaystyle B$ $\displaystyle \equiv$ $\displaystyle -\frac{e^{2}}{4\pi\epsilon_{0}E}\ \ \ \ \ (7)$ $\displaystyle C$ $\displaystyle \equiv$ $\displaystyle -\frac{\hbar^{2}}{2mE}l\left(l+1\right) \ \ \ \ \ (8)$

The turning points ${r_{1}}$ and ${r_{2}}$ can be found from the roots of the term inside the square root in 6:

 $\displaystyle r_{1}$ $\displaystyle =$ $\displaystyle \frac{B-\sqrt{B^{2}-4C}}{2}\ \ \ \ \ (9)$ $\displaystyle r_{2}$ $\displaystyle =$ $\displaystyle \frac{B+\sqrt{B^{2}-4C}}{2} \ \ \ \ \ (10)$

so we can factor the quadratic to get

 $\displaystyle \int_{r_{1}}^{r_{2}}p\left(r\right)\; dr$ $\displaystyle =$ $\displaystyle \sqrt{-2mE}\int_{r_{1}}^{r_{2}}\frac{1}{r}\sqrt{\left(r-r_{1}\right)\left(r_{2}-r\right)}dr\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mE}\frac{\pi}{2}\left(\sqrt{r_{2}}-\sqrt{r_{1}}\right)^{2}\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mE}\frac{\pi}{2}\left(r_{1}+r_{2}-2\sqrt{r_{1}r_{2}}\right)\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mE}\frac{\pi}{2}\left(B-2\sqrt{C}\right)\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mE}\frac{\pi}{2}\left(-\frac{e^{2}}{4\pi\epsilon_{0}E}-2\sqrt{-\frac{\hbar^{2}}{2mE}l\left(l+1\right)}\right)\ \ \ \ \ (15)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\pi}{2}\left(-\frac{e^{2}\sqrt{2m}}{4\pi\epsilon_{0}\sqrt{-E}}-2\hbar\sqrt{l\left(l+1\right)}\right)\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\pi\hbar \ \ \ \ \ (17)$

where we’ve used the integral given in Griffiths’s question to get the second line, and the last line uses 4. Solving 16 for ${E}$ gives the WKB energy levels:

 $\displaystyle E$ $\displaystyle =$ $\displaystyle -\left(\frac{e^{2}}{4\pi\epsilon_{0}}\right)^{2}\frac{m}{2\hbar^{2}\left(n-\frac{1}{2}+\sqrt{l\left(l+1\right)}\right)^{2}}\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{-13.6\mbox{ eV}}{\left(n-\frac{1}{2}+\sqrt{l\left(l+1\right)}\right)^{2}} \ \ \ \ \ (19)$

since the ground state of hydrogen is given by the Bohr formula

$\displaystyle E_{0}=-\left(\frac{e^{2}}{4\pi\epsilon_{0}}\right)^{2}\frac{m}{2\hbar^{2}}=-13.6\mbox{ eV} \ \ \ \ \ (20)$

## WKB approximation and the radial equation

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.13.

So far, we’ve applied the WKB approximation to one-dimensional potential problems. It might seem that that’s all we can manage since WKB is essentially a way of approximating the solution of a one-dimensional ODE. However, we can use it on 3-d problems in those cases where the solution is separable, such as spherically symmetric potentials. For such potentials, the general wave function can be written as the product of a radial function ${R\left(r\right)}$ and a spherical harmonic ${Y\left(\theta,\phi\right)}$.

With the substitution ${u\left(r\right)\equiv rR\left(r\right)}$ we found that the radial equation can be written as

$\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}u}{dr^{2}}+\left(V\left(r\right)+\frac{\hbar^{2}}{2m}\frac{l(l+1)}{r^{2}}\right)u=Eu \ \ \ \ \ (1)$

For the simplest case, we take ${l=0}$ so the equation becomes

$\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}u}{dr^{2}}+V\left(r\right)u=Eu \ \ \ \ \ (2)$

which has exactly the same form as the one-dimensional Schrödinger equation, so we should be able to use WKB to get approximate solutions.

There is one detail we need to work out first, though. In applying WKB to a true 1-d situation, the ${x}$ coordinate extended to infinity in both directions which allowed us to discard exponential terms that blow up as we approach one extreme or the other. In this equation, the ${r}$ coordinate starts at 0. One way of handling this is to assume that there is an infinite wall at ${r=0}$ so that ${u\left(0\right)=0}$. Since ${u\left(r\right)=rR\left(r\right)}$, this is a reasonable assumption, since it requires only that ${R\left(0\right)}$ is finite.

We can therefore consider a spherically symmetric potential well with an infinite barrier at ${r=0}$ and then some potential that increases from ${r=0}$ out to ${r=\infty}$. For a given energy ${E}$, there will be one turning point ${r_{2}}$ where ${E=V\left(r_{2}\right)}$.

For a turning point where ${V}$ is increasing, we’ve seen that the WKB functions on either side of the turning point are

$\displaystyle u\left(r\right)=\begin{cases} \frac{2D}{\sqrt{p\left(r\right)}}\sin\left[\int_{r}^{r_{2}}p\left(r'\right)\; dr'/\hbar+\pi/4\right] & rr_{2} \end{cases} \ \ \ \ \ (3)$

The requirement ${u\left(0\right)=0}$ means that the sine must be zero at ${r=0}$, so

 $\displaystyle \int_{0}^{r_{2}}p\left(r\right)\; dr/\hbar+\frac{\pi}{4}$ $\displaystyle =$ $\displaystyle n\pi\ \ \ \ \ (4)$ $\displaystyle \int_{0}^{r_{2}}p\left(r\right)\; dr$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{4}\right)\pi\hbar \ \ \ \ \ (5)$

Example We can apply this formula to the potential

$\displaystyle V\left(r\right)=V_{0}\ln\frac{r}{a} \ \ \ \ \ (6)$

The turning point is defined by

$\displaystyle E=V_{0}\ln\frac{r_{2}}{a} \ \ \ \ \ (7)$

so the integral 5 is

 $\displaystyle \sqrt{2m}\int_{0}^{r_{2}}\sqrt{E-V_{0}\ln\frac{r}{a}}dr$ $\displaystyle =$ $\displaystyle \sqrt{2m}\int_{0}^{r_{2}}\sqrt{V_{0}\ln\frac{r_{2}}{a}-V_{0}\ln\frac{r}{a}}dr\ \ \ \ \ (8)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{2mV_{0}}\int_{0}^{r_{2}}\sqrt{\ln\frac{r_{2}}{r}}dr \ \ \ \ \ (9)$

We can use the substitution

 $\displaystyle v$ $\displaystyle =$ $\displaystyle \ln\frac{r_{2}}{r}\ \ \ \ \ (10)$ $\displaystyle dv$ $\displaystyle =$ $\displaystyle \frac{r}{r_{0}}\left(-\frac{r_{0}}{r^{2}}\right)dr\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{1}{r}dr\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{e^{v}}{r_{0}}dr \ \ \ \ \ (13)$

The limits on the integral in terms of ${v}$ are

 $\displaystyle r$ $\displaystyle =$ $\displaystyle 0\rightarrow u=\infty\ \ \ \ \ (14)$ $\displaystyle r$ $\displaystyle =$ $\displaystyle r_{2}\rightarrow u=0 \ \ \ \ \ (15)$

so the integral transforms as

 $\displaystyle \sqrt{2mV_{0}}\int_{0}^{r_{2}}\sqrt{\ln\frac{r_{2}}{r}}dr$ $\displaystyle =$ $\displaystyle r_{2}\sqrt{2mV_{0}}\int_{0}^{\infty}\sqrt{v}e^{-v}dv\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle r_{2}\sqrt{2mV_{0}}\Gamma\left(\frac{3}{2}\right)\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2\pi mV_{0}}r_{2}}{2}\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{4}\right)\pi\hbar \ \ \ \ \ (19)$

where we used 5 in the last line.

To get the allowed energies we can substitute for ${r_{2}}$ using 7:

 $\displaystyle r_{2}$ $\displaystyle =$ $\displaystyle ae^{E/V_{0}}\ \ \ \ \ (20)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2\pi}{mV_{0}}}\left(n-\frac{1}{4}\right)\hbar\ \ \ \ \ (21)$ $\displaystyle E_{n}$ $\displaystyle =$ $\displaystyle V_{0}\ln\left(\sqrt{\frac{2\pi}{mV_{0}}}\left(n-\frac{1}{4}\right)\frac{\hbar}{a}\right) \ \ \ \ \ (22)$

The spacing between successive energy levels is

 $\displaystyle E_{n+1}-E_{n}$ $\displaystyle =$ $\displaystyle V_{0}\left[\ln\left(\sqrt{\frac{2\pi}{mV_{0}}}\left(n+\frac{3}{4}\right)\frac{\hbar}{a}\right)-\ln\left(\sqrt{\frac{2\pi}{mV_{0}}}\left(n-\frac{1}{4}\right)\frac{\hbar}{a}\right)\right]\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle V_{0}\ln\left(\frac{n+3/4}{n-1/4}\right) \ \ \ \ \ (24)$

## WKB approximation and the reflectionless potential

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.12.

We’ve seen that if we apply the WKB approximation to a potential well and require that the WKB wave functions match up in the region between the turning points, we get the condition

$\displaystyle \int_{x_{1}}^{x_{2}}p\left(x\right)\; dx=\left(n-\frac{1}{2}\right)\pi\hbar \ \ \ \ \ (1)$

where ${x_{1}}$ is the left turning point, ${x_{2}}$ is the right turning point and ${n=1,2,3,...}$ Another example of this process is an application to the reflectionless potential that we considered earlier:

$\displaystyle V(x)=-\frac{\hbar^{2}a^{2}}{m}\text{sech}^{2}(ax) \ \ \ \ \ (2)$

This gives a potential well centred at ${x=0}$ which approaches ${V\left(x\right)=0}$ as ${x\rightarrow\pm\infty}$. (It’s called ‘reflectionless’ since if ${E>0}$, an incident particle passes straight through the potential without any reflection.)

We’ll consider the bound state and compare the WKB approximation to the exact answer, which we worked out as

$\displaystyle E_{0}=-\frac{\hbar^{2}a^{2}}{2m} \ \ \ \ \ (3)$

To apply WKB, we need the turning points ${x_{1}}$ and ${x_{2}}$ where ${E=V}$. Since ${V}$ is even, ${x_{1}=-x_{2}}$ where

$\displaystyle E=-\frac{\hbar^{2}a^{2}}{m}\text{sech}^{2}(ax_{2}) \ \ \ \ \ (4)$

Since ${V}$ is an even function, we can write 1 as

$\displaystyle \int_{x_{1}}^{x_{2}}p\left(x\right)\; dx=2\int_{0}^{x_{2}}\sqrt{2m\left(E+\frac{\hbar^{2}a^{2}}{m}\text{sech}^{2}(ax)\right)}dx \ \ \ \ \ (5)$

Not surprisingly, Maple balks at this integral, so we need to help it along by trying a substitution. We can try

 $\displaystyle z$ $\displaystyle \equiv$ $\displaystyle \mbox{sech}^{2}\left(ax\right)\ \ \ \ \ (6)$ $\displaystyle dz$ $\displaystyle =$ $\displaystyle -2a\mbox{sech}^{2}\left(ax\right)\tanh\left(ax\right)dx\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -2az\sqrt{1-z}dx \ \ \ \ \ (8)$

From 4, we get the limits in terms of ${z}$. For ${x=0}$, ${z=1}$ and for ${x=x_{2}}$ we have

 $\displaystyle z_{2}$ $\displaystyle =$ $\displaystyle \text{sech}^{2}(ax_{2})\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{mE}{\hbar^{2}a^{2}}\equiv b \ \ \ \ \ (10)$

so

 $\displaystyle \int_{x_{1}}^{x_{2}}p\left(x\right)\; dx$ $\displaystyle =$ $\displaystyle -2\int_{1}^{b}\sqrt{2m\left(E+\frac{\hbar^{2}a^{2}}{m}z\right)}\frac{dz}{2az\sqrt{1-z}}\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{2}\hbar\int_{b}^{1}\frac{\sqrt{z-b}}{z\sqrt{1-z}}dz \ \ \ \ \ (12)$

Maple can do this integral provided we assume that ${b>0}$ (which is true from 10 since ${E<0}$) and ${b\ne1}$ (${b=1}$ just gives zero for the integral anyway). We get

$\displaystyle \sqrt{2}\hbar\int_{b}^{1}\frac{\sqrt{z-b}}{z\sqrt{1-z}}dz=\sqrt{2}\pi\hbar\left(1-\sqrt{b}\right) \ \ \ \ \ (13)$

From 1 we have

 $\displaystyle \sqrt{2}\pi\hbar\left(1-\sqrt{b}\right)$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\pi\hbar\ \ \ \ \ (14)$ $\displaystyle n$ $\displaystyle =$ $\displaystyle \sqrt{2}\left(1-\sqrt{b}\right)+\frac{1}{2} \ \ \ \ \ (15)$

Since the smallest ${\sqrt{b}}$ can be is zero, the largest ${n}$ can be is the greatest integer less than or equal to ${\sqrt{2}+\frac{1}{2}=1.914}$, so the only possible value of ${n}$ is ${n=1}$. In that case

 $\displaystyle \sqrt{2}\left(1-\sqrt{b}\right)$ $\displaystyle =$ $\displaystyle \frac{1}{2}\ \ \ \ \ (16)$ $\displaystyle b$ $\displaystyle =$ $\displaystyle \left(1-\frac{1}{2\sqrt{2}}\right)^{2}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(\frac{9}{8}-\frac{1}{\sqrt{2}}\right)\ \ \ \ \ (18)$ $\displaystyle E$ $\displaystyle =$ $\displaystyle -\frac{\hbar^{2}a^{2}}{m}\left(\frac{9}{8}-\frac{1}{\sqrt{2}}\right)\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle \approx$ $\displaystyle -0.418\frac{\hbar^{2}a^{2}}{m} \ \ \ \ \ (20)$

Comparing this to 3, we see that WKB gives a reasonable approximation.

## WKB approximation and the power law potential

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.11.

We’ve seen that if we apply the WKB approximation to a potential well and require that the WKB wave functions match up in the region between the turning points, we get the condition

$\displaystyle \int_{x_{1}}^{x_{2}}p\left(x\right)\; dx=\left(n-\frac{1}{2}\right)\pi\hbar \ \ \ \ \ (1)$

where ${x_{1}}$ is the left turning point, ${x_{2}}$ is the right turning point and ${n=1,2,3,...}$ We’ve applied this result to the case of the harmonic oscillator, where the potential is ${V\left(x\right)=\frac{1}{2}m\omega^{2}x^{2}}$, but it’s interesting to work out the more general case where the potential is any power law:

$\displaystyle V\left(x\right)=\alpha\left|x\right|^{\nu} \ \ \ \ \ (2)$

where ${\alpha}$ and ${\nu}$ are positive constants. The turning points for a particle with energy ${E}$ are found by solving ${E=V\left(x\right)}$ so we get

$\displaystyle x_{1,2}=\pm\left(\frac{E}{\alpha}\right)^{1/\nu} \ \ \ \ \ (3)$

Since the potential is an even function, we can write 1 as

$\displaystyle \int_{x_{1}}^{x_{2}}p\left(x\right)\; dx=2\int_{0}^{\left(E/\alpha\right)^{1/\nu}}\sqrt{2m\left(E-\alpha x^{\nu}\right)}dx \ \ \ \ \ (4)$

Maple is unable to handle the integral as it stands, so we can help it by using the substitution

 $\displaystyle u$ $\displaystyle =$ $\displaystyle x^{\nu}\ \ \ \ \ (5)$ $\displaystyle du$ $\displaystyle =$ $\displaystyle \nu x^{\nu-1}dx\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \nu u^{1-1/\nu}dx\ \ \ \ \ (7)$ $\displaystyle dx$ $\displaystyle =$ $\displaystyle \frac{1}{\nu}u^{1/\nu-1}du \ \ \ \ \ (8)$

The turning point in the ${u}$ coordinates is ${E/\alpha}$, so the integral becomes

$\displaystyle 2\int_{0}^{\left(E/\alpha\right)^{1/\nu}}\sqrt{2m\left(E-\alpha x^{\nu}\right)}dx=\frac{2\sqrt{2m}}{\nu}\int_{0}^{E/\alpha}\sqrt{E-\alpha u}u^{1/\nu-1}du \ \ \ \ \ (9)$

Maple can do this integral, with the result

 $\displaystyle \frac{2\sqrt{2m}}{\nu}\int_{0}^{E/\alpha}\sqrt{E-\alpha u}u^{1/\nu-1}du$ $\displaystyle =$ $\displaystyle \frac{\sqrt{2\pi m}}{\nu}\frac{1}{\alpha^{1/\nu}}\frac{\Gamma\left(\frac{1}{\nu}\right)}{\Gamma\left(\frac{1}{\nu}+\frac{3}{2}\right)}E^{\frac{1}{\nu}+\frac{1}{2}}\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\pi\hbar \ \ \ \ \ (11)$

Using the gamma function identity

$\displaystyle \Gamma\left(z+1\right)=z\Gamma\left(z\right) \ \ \ \ \ (12)$

with ${z=1/\nu}$, we can solve for ${E}$ and simplify the expression slightly:

 $\displaystyle E_{n}$ $\displaystyle =$ $\displaystyle \left[\left(n-\frac{1}{2}\right)\pi\hbar\frac{\Gamma\left(\frac{1}{\nu}+\frac{3}{2}\right)}{\Gamma\left(\frac{1}{\nu}\right)}\frac{\nu\alpha^{1/\nu}}{\sqrt{2\pi m}}\right]^{\frac{2\nu}{\nu+2}}\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left[\left(n-\frac{1}{2}\right)\pi\hbar\frac{\Gamma\left(\frac{1}{\nu}+\frac{3}{2}\right)}{\Gamma\left(\frac{1}{\nu}+1\right)}\frac{\alpha^{1/\nu}}{\sqrt{2\pi m}}\right]^{\frac{2\nu}{\nu+2}} \ \ \ \ \ (14)$

We also have

$\displaystyle \left(\alpha^{1/\nu}\right)^{\frac{2\nu}{\nu+2}}=\alpha^{\frac{2}{\nu+2}}=\alpha^{\frac{\nu+2-\nu}{\nu+2}}=\alpha\cdot\alpha^{-\frac{\nu}{\nu+2}}=\alpha\cdot\left[\left(\alpha\right)^{-1/2}\right]^{\frac{2\nu}{\nu+2}} \ \ \ \ \ (15)$

so

$\displaystyle E_{n}=\alpha\left[\left(n-\frac{1}{2}\right)\hbar\frac{\Gamma\left(\frac{1}{\nu}+\frac{3}{2}\right)}{\Gamma\left(\frac{1}{\nu}+1\right)}\frac{\sqrt{\pi}}{\sqrt{2m\alpha}}\right]^{\frac{2\nu}{\nu+2}} \ \ \ \ \ (16)$

In the special case of the harmonic oscillator, ${\nu=2}$ and ${\alpha=\frac{1}{2}m\omega^{2}}$. The values of the gamma function required here are

 $\displaystyle \Gamma\left(2\right)$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (17)$ $\displaystyle \Gamma\left(\frac{3}{2}\right)$ $\displaystyle =$ $\displaystyle \frac{1}{2}\sqrt{\pi} \ \ \ \ \ (18)$

so we get

$\displaystyle E_{n}=\left(n-\frac{1}{2}\right)\hbar\omega \ \ \ \ \ (19)$

## WKB approximation – analysis of the overlap region near a turning point

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.8.

We’ve looked at the WKB approximation at the points called turning points where the particle’s total energy is close to the potential energy. A fundamental problem with WKB is that it breaks down at turning points, so we can’t connect the functions on either side by requiring that they are continuous at the turning point itself. We connected the WKB wave functions on either side of a turning point by defining a patching function ${\psi_{p}}$ which is a solution of the Schrödinger equation for a linearized potential of the form

$\displaystyle V\left(x\right)\approx E+V'\left(0\right)x \ \ \ \ \ (1)$

where the turning point occurs at ${x=0}$. In order to apply boundary conditions we assumed that there was a region on either side of the turning point where we were sufficiently far from ${x=0}$ for WKB to be a good approximation, yet close enough to ${x=0}$ that 1 is a good approximation to the potential. Here we’ll examine the harmonic oscillator to see if this assumption is valid.

We’ll look specifically at the turning point ${x_{2}}$, where the potential is increasing, so ${V'\left(x_{2}\right)>0}$. Using the exact potential:

$\displaystyle V\left(x\right)=\frac{1}{2}m\omega^{2}x^{2} \ \ \ \ \ (2)$

we can find ${x_{2}}$:

 $\displaystyle E_{n}$ $\displaystyle =$ $\displaystyle \frac{1}{2}m\omega^{2}x_{2}^{2}\ \ \ \ \ (3)$ $\displaystyle x_{2}$ $\displaystyle =$ $\displaystyle \frac{1}{\omega}\sqrt{\frac{2E_{n}}{m}}\ \ \ \ \ (4)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\left(2n-1\right)\hbar}{m\omega}} \ \ \ \ \ (5)$

where we’ve used the formula for the energy levels of the harmonic oscillator

$\displaystyle E_{n}=\left(n-\frac{1}{2}\right)\hbar\omega \ \ \ \ \ (6)$

with ${n=1,2,3,...}$

The linearized potential for a point ${x=x_{2}+d}$ is

 $\displaystyle V_{lin}$ $\displaystyle =$ $\displaystyle E_{n}+m\omega^{2}x_{2}d\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\hbar\omega+\sqrt{\left(2n-1\right)\hbar m\omega^{3}}d \ \ \ \ \ (8)$

and the exact potential at the same point is

 $\displaystyle V\left(x_{2}+d\right)$ $\displaystyle =$ $\displaystyle \frac{1}{2}m\omega^{2}\left(x_{2}+d\right)^{2}\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{2}m\omega^{2}\left(\sqrt{\frac{\left(2n-1\right)\hbar}{m\omega}}+d\right)^{2} \ \ \ \ \ (10)$

Suppose we want to find the largest value of ${d}$ such that the difference between the exact and linearized potentials is within 1%, that is

$\displaystyle \frac{V\left(x_{2}+d\right)-V_{lin}\left(x_{2}+d\right)}{V\left(x_{2}+d\right)}=0.01 \ \ \ \ \ (11)$

Solving for ${d}$ (which requires solving a quadratic equation) we get

$\displaystyle d=0.1\sqrt{\frac{\left(2n-1\right)\hbar}{m\omega}} \ \ \ \ \ (12)$

Thus the larger the energy level ${n}$, the further we can go from the turning point ${x_{2}}$ and still have a good approximation using the linearized potential.

The patching wave function in this region is given by the Airy function ${Ai\left(z\right)=Ai\left(\alpha\left(x-x_{2}\right)\right)=Ai\left(\alpha d\right)}$. From an analysis of the Airy functions, it is known that the asymptotic form for large ${z}$ is within 1% of the true value provided ${z=\alpha d>5}$. From our earlier analysis, we have

 $\displaystyle \alpha$ $\displaystyle =$ $\displaystyle \left(\frac{2m}{\hbar^{2}}V'\left(x_{2}\right)\right)^{1/3}\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(\frac{2m}{\hbar^{2}}\sqrt{\left(2n-1\right)\hbar m\omega^{3}}\right)^{1/3} \ \ \ \ \ (14)$

Therefore

 $\displaystyle \alpha d$ $\displaystyle =$ $\displaystyle 0.1\left(\frac{2m}{\hbar^{2}}\sqrt{\left(2n-1\right)\hbar m\omega^{3}}\right)^{1/3}\sqrt{\frac{\left(2n-1\right)\hbar}{m\omega}}\ \ \ \ \ (15)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 0.1\left(2\right)^{1/3}\left(2n-1\right)^{2/3} \ \ \ \ \ (16)$

Requiring ${\alpha d>5}$ results in ${n>125.5}$, so the smallest value of ${n}$ for which the asymptotic form of the Airy function ${Ai\left(\alpha d\right)}$ is accurate to within 1% is ${n=126}$ (assuming that the ground state is given by ${n=1}$). Thus for large ${n}$, an overlap region in which both the linearized potential and the asymptotic form of the Airy function are valid does indeed exist.

In some cases (as we saw with the harmonic oscillator) however, WKB does actually give good results for smaller ${n}$.

## WKB approximation for a barrier with sloping sides

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.9.

We’ve studied the WKB approximation in the cases where the potential is increasing and decreasing at a finite rate (that is, the potential has no step functions). We’ve applied this to the case of a potential well with sloping sides, such as the harmonic oscillator. You might think that, since we’ve worked out the behaviour of the WKB solutions at both increasing and decreasing turning points, that we could apply these solutions to the case of a potential barrier with sloping sides. However, there is one fly in the ointment here: in considering the potential well, we assumed that ${E out to infinity on both sides of the well, meaning that we could throw away the positive exponential term in the solution. With a potential barrier of finite width, we can’t do this since ${E only within the barrier. We can, however, apply the same techniques to find a solution.

The starting point is the WKB functions for the three regions we’re interested in. We’ll take the barrier to be between ${x_{1}}$ on the left and ${x_{2}}$ on the right, so we have

$\displaystyle \psi\left(x\right)\approx\begin{cases} \frac{1}{\sqrt{p\left(x\right)}}\left[Ae^{i\int_{x}^{x_{1}}p\left(x'\right)dx'/\hbar}+Be^{-i\int_{x}^{x_{1}}p\left(x'\right)dx'/\hbar}\right] & xx_{2} \end{cases} \ \ \ \ \ (1)$

There’s no negative exponential in the last term, since we’re treating the usual scattering problem where a particle comes in from the left and is reflected and transmitted by the barrier. The key difference between the scattering problem and the potential well is that we can’t set ${C=0}$ in the middle equation for the reason given above.

We’ll consider first the turning point at ${x=x_{1}}$, where the potential is increasing. As before, we’ll set ${x_{1}=0}$ temporarily to ease the notation a bit. Recall that we approximate the potential at the turning points by

$\displaystyle V\left(x\right)\approx E+V'\left(0\right)x \ \ \ \ \ (2)$

and following the same procedure as before, we get the WKB function to the right of the left turning point (that is, for ${x}$ slightly larger than 0) to be

$\displaystyle \psi\approx\frac{1}{\sqrt{\hbar}\alpha^{3/4}x^{1/4}}\left[Ce^{2\left(\alpha x\right)^{3/2}/3}+De^{-2\left(\alpha x\right)^{3/2}/3}\right] \ \ \ \ \ (3)$

where

$\displaystyle \alpha\equiv\left(\frac{2mV'\left(0\right)}{\hbar^{2}}\right)^{1/3} \ \ \ \ \ (4)$

As before, we introduce a patching wave function ${\psi_{p}}$ to span the turning point, and find that it satisfies Airy’s equation

$\displaystyle \frac{d^{2}\psi_{p}}{dz^{2}}=z\psi_{p} \ \ \ \ \ (5)$

where ${z=\alpha x}$. This has the general solution

 $\displaystyle \psi_{p}\left(z\right)$ $\displaystyle =$ $\displaystyle aAi\left(z\right)+bBi\left(z\right)\ \ \ \ \ (6)$ $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle =$ $\displaystyle aAi\left(\alpha x\right)+bBi\left(\alpha x\right) \ \ \ \ \ (7)$

The asymptotic forms of the Airy functions ${Ai\left(z\right)}$ and ${Bi\left(z\right)}$ are

$\displaystyle Ai\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\sin\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{2\sqrt{\pi}z^{1/4}}e^{-2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (8)$

$\displaystyle Bi\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\cos\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{\sqrt{\pi}z^{1/4}}e^{2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (9)$

Restricting our attention to values of ${z}$ that are large enough that we can use the asymptotic forms of the Airy functions (but still satisfy 2), we get

$\displaystyle \psi_{p}\left(x\right)\approx\frac{1}{\sqrt{\pi}\left(\alpha x\right)^{1/4}}\left[\frac{a}{2}e^{-2\left(\alpha x\right)^{3/2}/3}+be^{2\left(\alpha x\right)^{3/2}/3}\right] \ \ \ \ \ (10)$

This time we can’t set ${b=0}$, but by equating coefficients between 3 and 10 we find

 $\displaystyle a$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}D\ \ \ \ \ (11)$ $\displaystyle b$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\pi}{\alpha\hbar}}C \ \ \ \ \ (12)$

Now looking at the region just to the left of the left-hand turning point, the WKB function is

$\displaystyle \psi\approx\frac{1}{\sqrt{\hbar}\alpha^{3/4}\left(-x\right)^{1/4}}\left[Ae^{2i\left(-\alpha x\right)^{3/2}/3}+Be^{-2i\left(-\alpha x\right)^{3/2}/3}\right] \ \ \ \ \ (13)$

The patching function comes out to (this time using the asymptotic forms of the Airy functions for large negative ${x}$):

 $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle \approx$ $\displaystyle \frac{1}{\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\left[a\sin\left(\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right)+b\cos\left(\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right)\right]\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{2\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\left[e^{i\pi/4}e^{2i\left(-\alpha x\right)^{3/2}/3}\left(b-ai\right)+e^{-i\pi/4}e^{-2i\left(-\alpha x\right)^{3/2}/3}\left(b+ai\right)\right] \ \ \ \ \ (15)$

Comparing coefficients gives

 $\displaystyle b-ai$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}Ae^{-i\pi/4}\ \ \ \ \ (16)$ $\displaystyle b+ai$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}Be^{i\pi/4} \ \ \ \ \ (17)$

Solving for ${a}$ and ${b}$ and comparing with 11 and 12 we get

 $\displaystyle C$ $\displaystyle =$ $\displaystyle Ae^{-i\pi/4}+Be^{i\pi/4}\ \ \ \ \ (18)$ $\displaystyle D$ $\displaystyle =$ $\displaystyle \frac{1}{2i}\left(-Ae^{-i\pi/4}+Be^{i\pi/4}\right)\ \ \ \ \ (19)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle e^{i\pi/4}\left(\frac{C}{2}-iD\right)\ \ \ \ \ (20)$ $\displaystyle B$ $\displaystyle =$ $\displaystyle e^{-i\pi/4}\left(\frac{C}{2}+iD\right) \ \ \ \ \ (21)$

That relates the WKB functions to the left and right of the left hand turning point at ${x=x_{1}}$. Now we need to repeat the calculation for the turning point at ${x=x_{2}}$, where the potential is decreasing. First, we need to rewrite the WKB functions so that the middle equation in 1 is in terms of ${x_{2}}$ rather than ${x_{1}}$. We can do this by noting that an integral from ${x_{1}}$ to ${x}$ is the sum of integrals from ${x_{1}}$ to ${x_{2}}$ and then from ${x_{2}}$ to ${x}$. When we considered the WKB approximation in tunneling we saw that the transmission coefficient was approximately

$\displaystyle T\equiv\frac{\left|F\right|^{2}}{\left|A\right|^{2}}\approx e^{-2\gamma}=e^{-2\int_{x_{1}}^{x_{2}}\left|p\left(x\right)\right|dx/\hbar} \ \ \ \ \ (22)$

We can therefore rewrite the WKB functions at the turning point in terms of ${\gamma}$ to get, for the region just to the left of ${x=x_{2}}$ (again, we can reposition ${x_{2}}$ to zero to make things easier):

 $\displaystyle \psi\left(x\right)$ $\displaystyle \approx$ $\displaystyle \frac{1}{\sqrt{\left|p\left(x\right)\right|}}\left[Ce^{\left(\int_{x_{1}}^{x_{2}}+\int_{x_{2}}^{x}\right)\left|p\left(x'\right)\right|dx'/\hbar}+De^{-\left(\int_{x_{1}}^{x_{2}}+\int_{x_{2}}^{x}\right)\left|p\left(x'\right)\right|dx'/\hbar}\right]\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\left|p\left(x\right)\right|}}\left[Ce^{\left(\int_{x_{1}}^{0}+\int_{0}^{x}\right)\left|p\left(x'\right)\right|dx'/\hbar}+De^{-\left(\int_{x_{1}}^{0}+\int_{0}^{x}\right)\left|p\left(x'\right)\right|dx'/\hbar}\right]\ \ \ \ \ (24)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\hbar}\alpha^{3/4}x^{1/4}}\left[Ce^{\gamma}e^{-2\left(-\alpha x\right)^{3/2}/3}+De^{-\gamma}e^{2\left(-\alpha x\right)^{3/2}/3}\right] \ \ \ \ \ (25)$

For the patching function in this region we have

 $\displaystyle \alpha$ $\displaystyle \equiv$ $\displaystyle \left(\frac{2m\left|V'\left(0\right)\right|}{\hbar^{2}}\right)^{1/3}\ \ \ \ \ (26)$ $\displaystyle z$ $\displaystyle \equiv$ $\displaystyle -\alpha x \ \ \ \ \ (27)$

so the asymptotic form of the Airy functions for large positive ${z}$ (since ${x<0}$) can be used to get

$\displaystyle \psi_{p}\left(x\right)\approx\frac{1}{\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\left[\frac{c}{2}e^{-2\left(-\alpha x\right)^{3/2}/3}+de^{2\left(-\alpha x\right)^{3/2}/3}\right] \ \ \ \ \ (28)$

Comparing coefficients between ${\psi}$ and ${\psi_{p}}$ we get

 $\displaystyle c$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}Ce^{\gamma}\ \ \ \ \ (29)$ $\displaystyle d$ $\displaystyle =$ $\displaystyle \sqrt{\frac{\pi}{\alpha\hbar}}De^{-\gamma} \ \ \ \ \ (30)$

Looking now at the region just to the right of ${x_{2}}$ we get for the WKB function

$\displaystyle \psi\left(x\right)\approx\frac{1}{\sqrt{\hbar}\alpha^{3/4}x^{1/4}}\left[Fe^{2i\left(\alpha x\right)^{3/2}/3}\right] \ \ \ \ \ (31)$

Using the asymptotic forms of the Airy functions for large negative ${z}$ we get

 $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle \approx$ $\displaystyle \frac{1}{\sqrt{\pi}\left(\alpha x\right)^{1/4}}\left[c\sin\left(\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right)+d\cos\left(\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right)\right]\ \ \ \ \ (32)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{2\sqrt{\pi}\left(\alpha x\right)^{1/4}}\left[e^{i\pi/4}e^{2i\left(\alpha x\right)^{3/2}/3}\left(d-ci\right)+e^{-i\pi/4}e^{-2i\left(\alpha x\right)^{3/2}/3}\left(d+ci\right)\right] \ \ \ \ \ (33)$

Comparing coefficients:

 $\displaystyle d+ci$ $\displaystyle =$ $\displaystyle 0\ \ \ \ \ (34)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}\left(\frac{D}{2}e^{-\gamma}+iCe^{\gamma}\right)\ \ \ \ \ (35)$ $\displaystyle D$ $\displaystyle =$ $\displaystyle -2ie^{2\gamma}C\ \ \ \ \ (36)$ $\displaystyle d-ci$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}Fe^{-i\pi/4}\ \ \ \ \ (37)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{4\pi}{\alpha\hbar}}\left(\frac{D}{2}e^{-\gamma}-iCe^{\gamma}\right)\ \ \ \ \ (38)$ $\displaystyle F$ $\displaystyle =$ $\displaystyle e^{i\pi/4}\left(\frac{D}{2}e^{-\gamma}-iCe^{\gamma}\right)\ \ \ \ \ (39)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -2Cie^{\gamma}e^{i\pi/4} \ \ \ \ \ (40)$

We’ve thus managed to express ${D}$ and ${F}$ in terms of ${C}$. We can find the transmission coefficient by using 20 to get

 $\displaystyle A$ $\displaystyle =$ $\displaystyle e^{i\pi/4}\left(\frac{C}{2}-iD\right)\ \ \ \ \ (41)$ $\displaystyle$ $\displaystyle =$ $\displaystyle e^{i\pi/4}C\left(\frac{1}{2}-2e^{2\gamma}\right)\ \ \ \ \ (42)$ $\displaystyle T$ $\displaystyle =$ $\displaystyle \frac{\left|F\right|^{2}}{\left|A\right|^{2}}\ \ \ \ \ (43)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{4e^{2\gamma}}{\left(\frac{1}{2}-2e^{2\gamma}\right)^{2}}\ \ \ \ \ (44)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{e^{-2\gamma}}{\left(\frac{e^{-2\gamma}}{4}-1\right)^{2}} \ \ \ \ \ (45)$

For a wide or high barrier, ${\gamma}$ gets very large, so the denominator tends to 1, and we get

$\displaystyle T\approx e^{-2\gamma} \ \ \ \ \ (46)$

which agrees with our earlier result 22 which was derived by ignoring the effects at the turning points. (Well, actually, it wasn’t ‘derived’ as such; it was more of a plausibility argument.)

## WKB approximation of the harmonic oscillator

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.7.

We’ve seen that if we apply the WKB approximation to a region where the potential is increasing we get

$\displaystyle \psi\left(x\right)\approx\begin{cases} \frac{2D}{\sqrt{p\left(x\right)}}\sin\left[\int_{x}^{x_{2}}p\; dx'/\hbar+\pi/4\right] & xx_{2} \end{cases} \ \ \ \ \ (1)$

where ${x_{2}}$ is the turning point, where the particle’s energy ${E=V\left(x_{2}\right)}$ with ${V'\left(x_{2}\right)>0}$ and ${D}$ is a normalization constant.

At a point ${x_{1}}$ where the potential is decreasing and ${E=V\left(x_{1}\right)}$ with ${V'\left(x_{1}\right)<0}$ we get

$\displaystyle \psi\left(x\right)\approx\begin{cases} \frac{2D'}{\sqrt{p\left(x\right)}}\sin\left[\int_{x_{1}}^{x}p\; dx'/\hbar+\pi/4\right] & x>x_{1}\\ \frac{D'}{\sqrt{\left|p\left(x\right)\right|}}\exp\left[-\int_{x}^{x_{1}}\left|p\left(x'\right)\right|\; dx'/\hbar\right] & x

where ${D'}$ is another normalization constant.

If we’re applying the WKB approximation to a potential well where the energy ${E}$ is greater than ${V\left(x\right)}$ only in the region ${x_{1} then the wave function in this region can be written as the sine function from either of these two forms, that is

 $\displaystyle \psi\left(x\right)$ $\displaystyle \approx$ $\displaystyle \frac{2D}{\sqrt{p\left(x\right)}}\sin\left[\int_{x}^{x_{2}}p\; dx'/\hbar+\pi/4\right]\ \ \ \ \ (3)$ $\displaystyle$ $\displaystyle \approx$ $\displaystyle \frac{2D'}{\sqrt{p\left(x\right)}}\sin\left[\int_{x_{1}}^{x}p\; dx'/\hbar+\pi/4\right] \ \ \ \ \ (4)$

Because the sine is an odd function, we can write the second form as

$\displaystyle \psi\left(x\right)\approx-\frac{2D'}{\sqrt{p\left(x\right)}}\sin\left[-\int_{x_{1}}^{x}p\; dx'/\hbar-\pi/4\right] \ \ \ \ \ (5)$

The zeroes of the sines must match up between these two forms which means the arguments of these two sines must be equal up to a multiple of ${\pi}$:

 $\displaystyle \frac{1}{\hbar}\int_{x}^{x_{2}}p\; dx'+\frac{\pi}{4}$ $\displaystyle =$ $\displaystyle -\frac{1}{\hbar}\int_{x_{1}}^{x}p\; dx'-\frac{\pi}{4}+n\pi\ \ \ \ \ (6)$ $\displaystyle \left(\int_{x_{1}}^{x}+\int_{x}^{x_{2}}\right)p\; dx'$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\pi\hbar\ \ \ \ \ (7)$ $\displaystyle \int_{x_{1}}^{x_{2}}p\; dx$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\pi\hbar \ \ \ \ \ (8)$

where ${n=1,2,3,...}$ (we have to start at ${n=1}$ rather than ${n=0}$ to keep the integral positive).

We can make the two approximate forms of ${\psi}$ equal if we also set ${D'=-D}$.

Example As a simple illustration of this, we consider the harmonic oscillator, with a potential

 $\displaystyle V\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{1}{2}kx^{2}\ \ \ \ \ (9)$ $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{2m\left(E-\frac{1}{2}kx^{2}\right)} \ \ \ \ \ (10)$

In this case, the turning points are

 $\displaystyle x_{1}$ $\displaystyle =$ $\displaystyle -\sqrt{\frac{2E}{k}}\ \ \ \ \ (11)$ $\displaystyle x_{2}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2E}{k}} \ \ \ \ \ (12)$

and

 $\displaystyle \int_{x_{1}}^{x_{2}}p\; dx$ $\displaystyle =$ $\displaystyle \sqrt{2m}\int_{-\sqrt{2E/k}}^{\sqrt{2E/k}}\sqrt{E-\frac{kx^{2}}{2}}dx\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \pi E\sqrt{\frac{m}{k}}\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\pi\hbar\ \ \ \ \ (15)$ $\displaystyle E$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\sqrt{\frac{k}{m}}\hbar\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{2}\right)\hbar\omega \ \ \ \ \ (17)$

Since ${n}$ starts at 1, this gives the same sequence of energies as the exact analysis.

## WKB approximation at a turning point with decreasing potential

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.9.

In the last post, we worked out the WKB approximation on either side of a turning point where the potential is increasing. We can do the same analysis for a turning point where the potential is decreasing. The calculations are almost exactly the same as in the first case.

The wave function has the form

$\displaystyle \psi\left(x\right)\approx\frac{1}{\sqrt{p\left(x\right)}}\left(Be^{i\int p\; dx/\hbar}+Ce^{-i\int p\; dx/\hbar}\right) \ \ \ \ \ (1)$

if the particle’s energy ${E>V\left(x\right)}$, where

$\displaystyle p\left(x\right)=\sqrt{2m\left(E-V\left(x\right)\right)} \ \ \ \ \ (2)$

is the momentum of the particle. In the tunneling case, where ${E, the WKB approximation gives

$\displaystyle \psi\left(x\right)\approx\frac{1}{\sqrt{\left|p\right|}}\left(De^{-\frac{1}{\hbar}\int\left|p\left(x\right)\right|dx}+Fe^{\frac{1}{\hbar}\int\left|p\left(x\right)\right|dx}\right) \ \ \ \ \ (3)$

In the decreasing potential case, we have ${E for ${x and ${E>V\left(x\right)}$ for ${x>x_{1}}$. As before, we’ll move the turning point so that ${x_{1}=0}$ and shift if back after the analysis. In this case, the limits on the integrals above give the results

$\displaystyle \psi\left(x\right)=\begin{cases} \frac{1}{\sqrt{\left|p\right|}}\left(De^{-\frac{1}{\hbar}\int_{x}^{0}\left|p\left(x'\right)\right|dx'}+Fe^{\frac{1}{\hbar}\int_{x}^{0}\left|p\left(x'\right)\right|dx'}\right) & x<0\\ \frac{1}{\sqrt{p\left(x\right)}}\left(Be^{i\int_{0}^{x}p\left(x'\right)\; dx'/\hbar}+Ce^{-i\int_{0}^{x}p\left(x'\right)\; dx'/\hbar}\right) & x>0 \end{cases} \ \ \ \ \ (4)$

As the integrand for ${x<0}$ is non-negative, we must have ${F=0}$ to keep ${\psi}$ finite as ${x\rightarrow-\infty}$, so we have

$\displaystyle \psi\left(x\right)=\begin{cases} \frac{1}{\sqrt{\left|p\right|}}\left(De^{-\frac{1}{\hbar}\int_{x}^{0}\left|p\left(x'\right)\right|dx'}\right) & x<0\\ \frac{1}{\sqrt{p\left(x\right)}}\left(Be^{i\int_{0}^{x}p\left(x'\right)\; dx'/\hbar}+Ce^{-i\int_{0}^{x}p\left(x'\right)\; dx'/\hbar}\right) & x>0 \end{cases} \ \ \ \ \ (5)$

As before, we assume that around ${x=0}$

$\displaystyle V\left(x\right)\approx E+V'\left(0\right)x \ \ \ \ \ (6)$

where now ${V'\left(0\right)<0}$ since the potential is decreasing. We now define the patching function ${\psi_{p}}$ that satisfies the Schrödinger equation with this approximate potential near ${x=0}$.

 $\displaystyle \frac{d^{2}\psi_{p}}{dx^{2}}$ $\displaystyle =$ $\displaystyle \frac{2mV'\left(0\right)}{\hbar^{2}}x\psi_{p}\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{2m\left|V'\left(0\right)\right|}{\hbar^{2}}x\psi_{p} \ \ \ \ \ (8)$

By defining

 $\displaystyle \alpha$ $\displaystyle \equiv$ $\displaystyle \left(\frac{2m\left|V'\left(0\right)\right|}{\hbar^{2}}\right)^{1/3}\ \ \ \ \ (9)$ $\displaystyle z$ $\displaystyle \equiv$ $\displaystyle -\alpha x \ \ \ \ \ (10)$

we end up with Airy’s equation

$\displaystyle \frac{d^{2}\psi_{p}}{dz^{2}}=z\psi_{p} \ \ \ \ \ (11)$

As we saw earlier, this has the general solution

 $\displaystyle \psi_{p}\left(z\right)$ $\displaystyle =$ $\displaystyle aAi\left(z\right)+bBi\left(z\right)\ \ \ \ \ (12)$ $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle =$ $\displaystyle aAi\left(-\alpha x\right)+bBi\left(-\alpha x\right) \ \ \ \ \ (13)$

It is this function that we need to match up to 5.

To do this, we make a further assumption that the overlap regions (where we’re trying to match the WKB functions with ${\psi_{p}}$) are far enough from ${x=0}$ that we can use the asymptotic forms of the Airy functions. These forms are

$\displaystyle Ai\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\sin\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{2\sqrt{\pi}z^{1/4}}e^{-2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (14)$

$\displaystyle Bi\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\cos\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{\sqrt{\pi}z^{1/4}}e^{2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (15)$

The forms can’t be used if we get too close to ${x=0}$, but the hope is that there is still some range of ${x}$ where both the WKB wave function and ${\psi_{p}}$ are reasonable approximations to the true wave function. Under this assumption, we have from 13 for ${z=-\alpha x\gg0}$, or ${x\ll0}$:

 $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{a}{2\sqrt{\pi}z^{1/4}}e^{-2z^{3/2}/3}+\frac{b}{\sqrt{\pi}z^{1/4}}e^{2z^{3/2}/3}\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{a}{2\sqrt{\pi}\left(-\alpha x\right)^{1/4}}e^{-2\left(-\alpha x\right)^{3/2}/3}+\frac{b}{\sqrt{\pi}\left(\alpha x\right)^{1/4}}e^{2\left(-\alpha x\right)^{3/2}/3} \ \ \ \ \ (17)$

With the potential 6, the WKB function 5 for ${x<0}$ is

 $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{2m\left(E-V\left(x\right)\right)}\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mV'\left(0\right)x}\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \hbar\sqrt{x}\alpha^{3/2}\ \ \ \ \ (20)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\left|p\right|}}De^{-\frac{1}{\hbar}\int_{x}^{0}\left|p\left(x'\right)\right|dx'}\ \ \ \ \ (21)$ $\displaystyle \int_{0}^{x}\left|p\left(x'\right)\right|dx'$ $\displaystyle =$ $\displaystyle \sqrt{-2mV'\left(0\right)}\int_{x}^{0}\sqrt{-x'}dx'\ \ \ \ \ (22)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{2m\left|V'\left(0\right)\right|}\frac{2}{3}\left(-x\right)^{3/2}\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2}{3}\hbar\left(-\alpha x\right)^{3/2}\ \ \ \ \ (24)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{D}{\sqrt{\hbar}\alpha^{3/4}\left(-x\right)^{1/4}}e^{-2\left(-\alpha x\right)^{3/2}/3} \ \ \ \ \ (25)$

Matching this up with 17, we see that ${b=0}$ and

 $\displaystyle \frac{a}{2\sqrt{\pi}\alpha^{1/4}}$ $\displaystyle =$ $\displaystyle \frac{D}{\sqrt{\hbar}\alpha^{3/4}}\ \ \ \ \ (26)$ $\displaystyle a$ $\displaystyle =$ $\displaystyle 2\sqrt{\frac{\pi}{\hbar\alpha}}D \ \ \ \ \ (27)$

The patching function must therefore have the form

$\displaystyle \psi_{p}\left(x\right)=2\sqrt{\frac{\pi}{\hbar\alpha}}De^{-2\left(-\alpha x\right)^{3/2}/3} \ \ \ \ \ (28)$

for ${\alpha x\ll0}$.

We can do the same calculation for the overlap region for ${x>0}$, where ${E>V\left(x\right)}$. Since ${x>0}$, ${z=-\alpha x<0}$ and we can use the asymptotic form of ${Ai\left(x\right)}$ for large negative ${x}$ and we get

$\displaystyle \psi_{p}\left(x\right)=\frac{a}{\sqrt{\pi}\left(\alpha x\right)^{1/4}}\sin\left[\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right] \ \ \ \ \ (29)$

 $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{-2mV'\left(0\right)x}\ \ \ \ \ (30)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \hbar\alpha^{3/2}\sqrt{x}\ \ \ \ \ (31)$ $\displaystyle \int_{0}^{x}p\left(x'\right)dx'$ $\displaystyle =$ $\displaystyle \frac{2}{3}\hbar\left(\alpha x\right)^{3/2}\ \ \ \ \ (32)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\hbar}\alpha^{3/4}\left(x\right)^{1/4}}\left(Be^{2i\left(\alpha x\right)^{3/2}/3}+Ce^{-2i\left(\alpha x\right)^{3/2}/3}\right) \ \ \ \ \ (33)$

We can now compare this with 29 by writing the sine in terms of exponentials

$\displaystyle \psi_{p}\left(x\right)=\frac{a}{2i\sqrt{\pi}\left(\alpha x\right)^{1/4}}\left[e^{i\left[\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right]}-e^{-i\left[\frac{2}{3}\left(\alpha x\right)^{3/2}+\frac{\pi}{4}\right]}\right] \ \ \ \ \ (34)$

Comparing terms, we have

 $\displaystyle \frac{B}{\sqrt{\hbar}\alpha^{3/4}}$ $\displaystyle =$ $\displaystyle \frac{ae^{i\pi/4}}{2i\sqrt{\pi}\alpha^{1/4}}\ \ \ \ \ (35)$ $\displaystyle B$ $\displaystyle =$ $\displaystyle \frac{ae^{i\pi/4}\sqrt{\hbar\alpha}}{2i\sqrt{\pi}}\ \ \ \ \ (36)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -ie^{i\pi/4}D\ \ \ \ \ (37)$ $\displaystyle \frac{C}{\sqrt{\hbar}\alpha^{3/4}}$ $\displaystyle =$ $\displaystyle -\frac{ae^{-i\pi/4}}{2i\sqrt{\pi}\alpha^{1/4}}\ \ \ \ \ (38)$ $\displaystyle C$ $\displaystyle =$ $\displaystyle -\frac{ae^{-i\pi/4}\sqrt{\hbar\alpha}}{2i\sqrt{\pi}}\ \ \ \ \ (39)$ $\displaystyle$ $\displaystyle =$ $\displaystyle ie^{-i\pi/4}D \ \ \ \ \ (40)$

We now have the constants ${B}$ and ${C}$, belonging to the WKB wave function for ${x>0}$, in terms of the constant ${D}$ belonging to the WKB function for ${x<0}$, which is what we were after. Note that the WKB functions are still valid only for values of ${x}$ that maintain a respectable distance from the turning point; all we have done is connected the approximations on either side of the turning point. The final form of the WKB function is, for ${x>0}$:

 $\displaystyle \psi\left(x\right)$ $\displaystyle \cong$ $\displaystyle \frac{-iD}{\sqrt{p\left(x\right)}}\left(e^{i\left(\int p\; dx/\hbar+\pi/4\right)}-e^{-i\left(\int p\; dx/\hbar+\pi/4\right)}\right)\ \ \ \ \ (41)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2D}{\sqrt{p\left(x\right)}}\frac{1}{2i}\left(e^{i\left(\int p\; dx/\hbar+\pi/4\right)}-e^{-i\left(\int p\; dx/\hbar+\pi/4\right)}\right)\ \ \ \ \ (42)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2D}{\sqrt{p\left(x\right)}}\sin\left[\int_{0}^{x}p\; dx'/\hbar+\pi/4\right] \ \ \ \ \ (43)$

And for ${x<0}$

$\displaystyle \psi\left(x\right)\cong\frac{D}{\sqrt{\left|p\left(x\right)\right|}}\exp\left[-\int_{x}^{0}\left|p\left(x'\right)\right|\; dx'/\hbar\right] \ \ \ \ \ (44)$

If the turning point is at the general location of ${x_{1}}$, we can just replace the 0 in the limits of the integrals by ${x_{1}}$. The only difference between the cases of increasing and decreasing potentials is that the limits of integration are reversed.

## WKB approximation: turning points

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.6.

We’ve seen that the WKB approximation for the wave function has the form

$\displaystyle \psi\left(x\right)\approx\frac{1}{\sqrt{p\left(x\right)}}\left(Be^{i\int p\; dx/\hbar}+Ce^{-i\int p\; dx/\hbar}\right) \ \ \ \ \ (1)$

if the particle’s energy ${E>V\left(x\right)}$, where

$\displaystyle p\left(x\right)=\sqrt{2m\left(E-V\left(x\right)\right)} \ \ \ \ \ (2)$

is the momentum of the particle. In the tunneling case, where ${E, the WKB approximation gives

$\displaystyle \psi\left(x\right)\approx\frac{1}{\sqrt{\left|p\right|}}\left(De^{-\frac{1}{\hbar}\int\left|p\left(x\right)\right|dx}+Fe^{\frac{1}{\hbar}\int\left|p\left(x\right)\right|dx}\right) \ \ \ \ \ (3)$

If the region in which ${E extends out to ${x\rightarrow\infty}$, the positive exponential term must vanish, so we have ${F=0}$ and

$\displaystyle \psi\left(x\right)\approx\frac{1}{\sqrt{\left|p\right|}}De^{-\frac{1}{\hbar}\int\left|p\left(x\right)\right|dx} \ \ \ \ \ (4)$

The WKB approximation was derived by approximating the differential equation

$\displaystyle \frac{A^{\prime\prime}}{A}=\left(\phi'\right)^{2}-\frac{p^{2}}{\hbar^{2}} \ \ \ \ \ (5)$

by assuming that the LHS is very small compared to both terms on the RHS, in particular, that

$\displaystyle \frac{A''}{A}\ll\frac{p^{2}\left(x\right)}{\hbar^{2}} \ \ \ \ \ (6)$

Now suppose we have a potential that increases at some finite rate (that is, there are no vertical walls as in the infinite square well). Then if the particle’s energy is ${E}$, there is a point ${x_{2}}$, say, where ${E=V\left(x_{2}\right)}$ and at this point ${p\left(x_{2}\right)=0}$ so the condition 6 is not valid. Griffiths calls ${x_{2}}$ the turning point, and it’s clear that the WKB approximation breaks down at turning points. We can see this from 1 and 4 as well, since if ${p\left(x_{2}\right)=0}$, ${\psi\left(x\right)\rightarrow\infty}$ as ${x\rightarrow x_{2}}$, which can’t be allowed since an infinite wave function cannot be normalized. Because of this problem, we can’t satisfy the condition that the wave function must be continuous at the turning point.

The way around this is a bit of a fudge, but the way it works is as follows. If the potential varies continuously at the turning point (as all physical potential do), then for a small region around the turning point we can approximate the potential by its tangent line at ${x=x_{2}}$. (For convenience, we’ll take ${x_{2}=0}$ in the argument that follows.) That is, around ${x=0}$ we assume

$\displaystyle V\left(x\right)\approx E+V'\left(0\right)x \ \ \ \ \ (7)$

With a linear potential, we have a situation similar to the problem we solved earlier, with a particle moving in the gravitational field at the Earth’s surface. The Schrödinger equation is

$\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}\psi_{p}}{dx^{2}}+\left(E+V'\left(0\right)x\right)\psi_{p}=E\psi_{p} \ \ \ \ \ (8)$

where ${\psi_{p}\left(x\right)}$ is a ‘patching’ wave function that spans a small distance either side of the turning point. The idea is that this patching wave function is continuous across the turning point, and that we can make it continuous with one WKB wave function on one side and with the other WKB function on the other side of the turning point.

To make things definite, suppose the potential is increasing at ${x=0}$, so that ${V'\left(0\right)>0}$, ${E>V\left(x\right)}$ for ${x<0}$ and ${E for ${x>0}$. Then the WKB function 1 applies for ${x<0}$ and the function 4 applies for ${x>0}$. The idea is that we should be able to make ${\psi_{p}\left(x\right)}$ equal to 1 for some range of ${x<0}$ and also make ${\psi_{p}\left(x\right)}$ equal to 4 for some range of ${x>0}$.

It’s important to understand that we’re not going to make the wave functions match up right at ${x=0}$. Rather, we will take a range of ${x}$ on the ${x<0}$ side of the turning point that is close enough to the turning point that we can still use 7 as an approximation for the potential (and thus that the patching wave function is still valid), but far enough from ${x=0}$ that the WKB approximation is also valid. It’s a bit of a leap of faith that such a region can actually be found, although in practice it turns out that in most cases it can be.

We must first find ${\psi_{p}\left(x\right)}$. We can rewrite 8 as

$\displaystyle \frac{d^{2}\psi_{p}}{dx^{2}}=\frac{2mV'\left(0\right)}{\hbar^{2}}x\psi_{p} \ \ \ \ \ (9)$

By defining

 $\displaystyle \alpha$ $\displaystyle \equiv$ $\displaystyle \left(\frac{2mV'\left(0\right)}{\hbar^{2}}\right)^{1/3}\ \ \ \ \ (10)$ $\displaystyle z$ $\displaystyle \equiv$ $\displaystyle \alpha x \ \ \ \ \ (11)$

we end up with Airy’s equation

$\displaystyle \frac{d^{2}\psi_{p}}{dz^{2}}=z\psi_{p} \ \ \ \ \ (12)$

As we saw earlier, this has the general solution

 $\displaystyle \psi_{p}\left(z\right)$ $\displaystyle =$ $\displaystyle aAi\left(z\right)+bBi\left(z\right)\ \ \ \ \ (13)$ $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle =$ $\displaystyle aAi\left(\alpha x\right)+bBi\left(\alpha x\right) \ \ \ \ \ (14)$

It is this function that we need to match up to 1 for ${x<0}$ and to 4 for ${x>0}$.

To do this, we make a further assumption that the overlap regions (where we’re trying to match the WKB functions with ${\psi_{p}}$) are far enough from ${x=0}$ that we can use the asymptotic forms of the Airy functions. These forms are

$\displaystyle Ai\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\sin\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{2\sqrt{\pi}z^{1/4}}e^{-2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (15)$

$\displaystyle Bi\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\cos\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{\sqrt{\pi}z^{1/4}}e^{2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (16)$

The forms can’t be used if we get too close to ${x=0}$, but the hope is that there is still some range of ${x}$ where both the WKB wave function and ${\psi_{p}}$ are reasonable approximations to the true wave function. Under this assumption, we have from 14 for ${z=\alpha x\gg0}$:

 $\displaystyle \psi_{p}\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{a}{2\sqrt{\pi}z^{1/4}}e^{-2z^{3/2}/3}+\frac{b}{\sqrt{\pi}z^{1/4}}e^{2z^{3/2}/3}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{a}{2\sqrt{\pi}\left(\alpha x\right)^{1/4}}e^{-2\left(\alpha x\right)^{3/2}/3}+\frac{b}{\sqrt{\pi}\left(\alpha x\right)^{1/4}}e^{2\left(\alpha x\right)^{3/2}/3} \ \ \ \ \ (18)$

With the potential 7, the WKB function 4 is

 $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{2m\left(E-V\left(x\right)\right)}\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{-2mV'\left(0\right)x}\ \ \ \ \ (20)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\left|p\right|}}De^{-\frac{1}{\hbar}\int_{0}^{x}\left|p\left(x'\right)\right|dx'}\ \ \ \ \ (21)$ $\displaystyle \int_{0}^{x}\left|p\left(x'\right)\right|dx'$ $\displaystyle =$ $\displaystyle \int_{0}^{x}\left|p\left(x'\right)\right|dx'\ \ \ \ \ (22)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{2mV'\left(0\right)}\int_{0}^{x}\sqrt{x'}dx'\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \sqrt{2mV'\left(0\right)}\frac{2}{3}x^{3/2}\ \ \ \ \ (24)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2}{3}\hbar\left(\alpha x\right)^{3/2}\ \ \ \ \ (25)$ $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{-2mV'\left(0\right)x}\ \ \ \ \ (26)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \hbar\alpha^{3/2}\sqrt{-x}\ \ \ \ \ (27)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{D}{\sqrt{\hbar}\alpha^{3/4}x^{1/4}}e^{-2\left(\alpha x\right)^{3/2}/3} \ \ \ \ \ (28)$

Matching this up with 18, we see that ${b=0}$ and

 $\displaystyle \frac{a}{2\sqrt{\pi}\alpha^{1/4}}$ $\displaystyle =$ $\displaystyle \frac{D}{\sqrt{\hbar}\alpha^{3/4}}\ \ \ \ \ (29)$ $\displaystyle a$ $\displaystyle =$ $\displaystyle 2\sqrt{\frac{\pi}{\hbar\alpha}}D \ \ \ \ \ (30)$

The patching function must therefore have the form

$\displaystyle \psi_{p}\left(x\right)=2\sqrt{\frac{\pi}{\hbar\alpha}}De^{-2\left(\alpha x\right)^{3/2}/3} \ \ \ \ \ (31)$

for ${\alpha x\gg0}$.

We can do the same calculation for the overlap region for ${x<0}$, where ${E>V\left(x\right)}$. Since ${x<0}$, we can use the asymptotic form of ${Ai\left(x\right)}$ for large negative ${x}$ and we get

$\displaystyle \psi_{p}\left(x\right)=\frac{a}{\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\sin\left[\frac{2}{3}\left(-\alpha x\right)^{3/2}+\frac{\pi}{4}\right] \ \ \ \ \ (32)$

 $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{-2mV'\left(0\right)x}\ \ \ \ \ (33)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \hbar\alpha^{3/2}\sqrt{-x}\ \ \ \ \ (34)$ $\displaystyle \int_{x}^{0}p\left(x\right)dx$ $\displaystyle =$ $\displaystyle \frac{2}{3}\hbar\alpha^{3/2}\left(-x\right)^{3/2}\ \ \ \ \ (35)$ $\displaystyle \psi\left(x\right)$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\hbar}\alpha^{3/4}\left(-x\right)^{1/4}}\left(Be^{2i\left(-\alpha x\right)^{3/2}/3}+Ce^{-2i\left(-\alpha x\right)^{3/2}/3}\right) \ \ \ \ \ (36)$

We can now compare this with 32 by writing the sine in terms of exponentials

$\displaystyle \psi_{p}\left(x\right)=\frac{a}{2i\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\left[e^{i\left[\frac{2}{3}\left(-\alpha x\right)^{3/2}+\frac{\pi}{4}\right]}-e^{-i\left[\frac{2}{3}\left(-\alpha x\right)^{3/2}+\frac{\pi}{4}\right]}\right] \ \ \ \ \ (37)$

Comparing terms, we have

 $\displaystyle \frac{B}{\sqrt{\hbar}\alpha^{3/4}}$ $\displaystyle =$ $\displaystyle \frac{ae^{i\pi/4}}{2i\sqrt{\pi}\alpha^{1/4}}\ \ \ \ \ (38)$ $\displaystyle B$ $\displaystyle =$ $\displaystyle \frac{ae^{i\pi/4}\sqrt{\hbar\alpha}}{2i\sqrt{\pi}}\ \ \ \ \ (39)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -ie^{i\pi/4}D\ \ \ \ \ (40)$ $\displaystyle \frac{C}{\sqrt{\hbar}\alpha^{3/4}}$ $\displaystyle =$ $\displaystyle -\frac{ae^{-i\pi/4}}{2i\sqrt{\pi}\alpha^{1/4}}\ \ \ \ \ (41)$ $\displaystyle C$ $\displaystyle =$ $\displaystyle -\frac{ae^{-i\pi/4}\sqrt{\hbar\alpha}}{2i\sqrt{\pi}}\ \ \ \ \ (42)$ $\displaystyle$ $\displaystyle =$ $\displaystyle ie^{-i\pi/4}D \ \ \ \ \ (43)$

We now have the constants ${B}$ and ${C}$, belonging to the WKB wave function for ${x<0}$, in terms of the constant ${D}$ belonging to the WKB function for ${x>0}$, which is what we were after. Note that the WKB functions are still valid only for values of ${x}$ that maintain a respectable distance from the turning point; all we have done is connected the approximations on either side of the turning point. The final form of the WKB function is, for ${x<0}$:

 $\displaystyle \psi\left(x\right)$ $\displaystyle \cong$ $\displaystyle \frac{-iD}{\sqrt{p\left(x\right)}}\left(e^{i\left(\int p\; dx/\hbar+\pi/4\right)}-e^{-i\left(\int p\; dx/\hbar+\pi/4\right)}\right)\ \ \ \ \ (44)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2D}{\sqrt{p\left(x\right)}}\frac{1}{2i}\left(e^{i\left(\int p\; dx/\hbar+\pi/4\right)}-e^{-i\left(\int p\; dx/\hbar+\pi/4\right)}\right)\ \ \ \ \ (45)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{2D}{\sqrt{p\left(x\right)}}\sin\left[\int_{x}^{0}p\; dx'/\hbar+\pi/4\right] \ \ \ \ \ (46)$

And for ${x>0}$

$\displaystyle \psi\left(x\right)\cong\frac{D}{\sqrt{\left|p\left(x\right)\right|}}\exp\left[-\int_{0}^{x}\left|p\left(x'\right)\right|\; dx'/\hbar\right] \ \ \ \ \ (47)$

If the turning point is at the general location of ${x_{2}}$, we can just replace the 0 in the limits of the integrals by ${x_{2}}$.

As an example, we can revisit the case of the particle moving under gravity at the Earth’s surface. The potential has a vertical wall at ${x=0}$ and for a particle of energy ${E}$, there is a turning point at

$\displaystyle x_{2}=\frac{E}{mg} \ \ \ \ \ (48)$

The potential can then be written as

 $\displaystyle p\left(x\right)$ $\displaystyle =$ $\displaystyle \sqrt{2m\left(E-mgx\right)}\ \ \ \ \ (49)$ $\displaystyle$ $\displaystyle =$ $\displaystyle m\sqrt{2g\left(x_{2}-x\right)} \ \ \ \ \ (50)$

We can get the allowed energy levels from the boundary condition that ${\psi\left(0\right)=0}$. From 46 (replacing the 0 in the limit by ${x_{2}}$) we have

 $\displaystyle \psi\left(0\right)$ $\displaystyle =$ $\displaystyle \frac{2D}{\sqrt{p\left(x\right)}}\sin\left[\int_{0}^{x_{2}}p\; dx'/\hbar+\pi/4\right]=0\ \ \ \ \ (51)$ $\displaystyle \int_{0}^{x_{2}}p\; dx'/\hbar+\pi/4$ $\displaystyle =$ $\displaystyle n\pi\ \ \ \ \ (52)$ $\displaystyle \int_{0}^{x_{2}}p\; dx'$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{4}\right)\pi\hbar \ \ \ \ \ (53)$

Doing the integral, we get

 $\displaystyle m\sqrt{2g}\int_{0}^{x_{2}}\sqrt{x_{2}-x'}dx'$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{4}\right)\pi\hbar\ \ \ \ \ (54)$ $\displaystyle \frac{m\sqrt{g}}{3}\left(2x_{2}\right)^{3/2}$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{4}\right)\pi\hbar\ \ \ \ \ (55)$ $\displaystyle \frac{m\sqrt{g}}{3}\left(2\frac{E_{n}}{mg}\right)^{3/2}$ $\displaystyle =$ $\displaystyle \left(n-\frac{1}{4}\right)\pi\hbar\ \ \ \ \ (56)$ $\displaystyle E_{n}$ $\displaystyle =$ $\displaystyle \frac{1}{8}\left(6\pi\hbar g\left(4n-1\right)\right)^{2/3}m^{1/3} \ \ \ \ \ (57)$

Using values from the previous problem (${g=9.8\mbox{ m s}^{-2}}$ and a mass of ${0.1\mbox{ kg}}$), we get the WKB estimate of the first four energy levels:

 $\displaystyle E_{1}$ $\displaystyle =$ $\displaystyle 8.738\times10^{-23}\mbox{ J}\ \ \ \ \ (58)$ $\displaystyle E_{2}$ $\displaystyle =$ $\displaystyle 1.537\times10^{-22}\mbox{ J}\ \ \ \ \ (59)$ $\displaystyle E_{3}$ $\displaystyle =$ $\displaystyle 2.078\times10^{-22}\mbox{ J}\ \ \ \ \ (60)$ $\displaystyle E_{4}$ $\displaystyle =$ $\displaystyle 2.555\times10^{-22}\mbox{ J} \ \ \ \ \ (61)$

Comparing with the exact results shows the WKB values are quite good in this case:

 $\displaystyle E_{1}$ $\displaystyle =$ $\displaystyle 8.805\times10^{-23}\mbox{ J}\ \ \ \ \ (62)$ $\displaystyle E_{2}$ $\displaystyle =$ $\displaystyle 1.539\times10^{-22}\mbox{ J}\ \ \ \ \ (63)$ $\displaystyle E_{3}$ $\displaystyle =$ $\displaystyle 2.079\times10^{-22}\mbox{ J}\ \ \ \ \ (64)$ $\displaystyle E_{4}$ $\displaystyle =$ $\displaystyle 2.556\times10^{-22}\mbox{ J} \ \ \ \ \ (65)$

From the previous problem, we can find the average height of the particle for a given energy level:

$\displaystyle \left\langle x_{n}\right\rangle =\frac{2E_{n}}{3mg} \ \ \ \ \ (66)$

If we want ${\left\langle x_{n}\right\rangle =1\mbox{ m}}$, then

 $\displaystyle \frac{1}{12g}\left(6\pi\hbar g\left(4n-1\right)\right)^{2/3}m^{-2/3}$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (67)$ $\displaystyle n$ $\displaystyle =$ $\displaystyle 1.6366\times10^{33} \ \ \ \ \ (68)$

## Airy functions and the bouncing electron

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 8.5.

In order to progress further with our study of the WKB approximation we’ll need to study a differential equation known as Airy’s equation, named after George Biddell Airy (1801 – 1892), an English mathematician and astronomer (he served as Astronomer Royal for many years). We’ll leave the application to the WKB approximation for the next post; here we’ll look at how Airy’s equation arises in quantum mechanics and do an example of its solution.

Suppose we have a particle of mass ${m}$ that moves vertically under the influence of gravity at the Earth’s surface, and that it bounces elastically off some hard surface which occupies the plane ${x=0}$. The potential is

$\displaystyle V\left(x\right)=\begin{cases} \infty & x<0\\ mgx & x>0 \end{cases} \ \ \ \ \ (1)$

The Schrödinger equation is

 $\displaystyle -\frac{\hbar^{2}}{2m}\frac{\partial^{2}\psi}{\partial x^{2}}+mgx\psi$ $\displaystyle =$ $\displaystyle E\psi\ \ \ \ \ (2)$ $\displaystyle \frac{\partial^{2}\psi}{\partial x^{2}}$ $\displaystyle =$ $\displaystyle \frac{2m^{2}g}{\hbar^{2}}\left(x-\frac{E}{mg}\right)\psi\ \ \ \ \ (3)$ $\displaystyle$ $\displaystyle \equiv$ $\displaystyle \alpha^{3}\left(x-\frac{E}{mg}\right)\psi\ \ \ \ \ (4)$ $\displaystyle \alpha$ $\displaystyle \equiv$ $\displaystyle \left(\frac{2m^{2}g}{\hbar^{2}}\right)^{1/3} \ \ \ \ \ (5)$

Now we use the substitution

 $\displaystyle z$ $\displaystyle \equiv$ $\displaystyle \alpha\left(x-\frac{E}{mg}\right)\ \ \ \ \ (6)$ $\displaystyle \frac{\partial^{2}\psi}{\partial x^{2}}$ $\displaystyle =$ $\displaystyle \alpha^{2}\frac{\partial^{2}\psi}{\partial z^{2}}\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \alpha^{3}\frac{z}{\alpha}\psi\ \ \ \ \ (8)$ $\displaystyle \frac{\partial^{2}\psi}{\partial z^{2}}$ $\displaystyle =$ $\displaystyle z\psi \ \ \ \ \ (9)$

Equation 9 is Airy’s equation. Its solutions are called Airy functions. As this is a second order differential equation, there are two linearly independent solutions, known as ${Ai\left(z\right)}$ and ${Bi\left(z\right)}$. They are not easily expressed in terms of anything simpler, but they do have some useful properties that aren’t too hard to state.

Both ${Ai\left(z\right)}$ and ${Bi\left(z\right)}$ oscillate about the ${z}$ axis for ${z<0}$, so there are a number of zeroes for these functions in that region. As ${z\rightarrow+\infty}$, ${Ai\left(z\right)\rightarrow0}$ and ${Bi\left(z\right)\rightarrow\infty}$. Their asymptotic forms are, for reference:

$\displaystyle Ai\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\sin\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{2\sqrt{\pi}z^{1/4}}e^{-2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (10)$

$\displaystyle Bi\left(z\right)\sim\begin{cases} \frac{1}{\sqrt{\pi}\left(-z\right)^{1/4}}\cos\left[\frac{2}{3}\left(-z\right)^{3/2}+\frac{\pi}{4}\right] & z\ll0\\ \frac{1}{\sqrt{\pi}z^{1/4}}e^{2z^{3/2}/3} & z\gg0 \end{cases} \ \ \ \ \ (11)$

For our example, we can exclude ${Bi\left(z\right)}$ since it blows up as ${z}$ gets large, so the general solution is

$\displaystyle \psi\left(z\right)=\begin{cases} 0 & z<-\frac{E}{mg}\\ aAi\left(z\right) & z>\frac{E}{mg} \end{cases} \ \ \ \ \ (12)$

or, in terms of ${x}$:

$\displaystyle \psi\left(\alpha\left(x-\frac{E}{mg}\right)\right)=\begin{cases} 0 & x<0\\ aAi\left(\alpha\left(x-\frac{E}{mg}\right)\right) & x>0 \end{cases} \ \ \ \ \ (13)$

To find the allowed energies, we impose the boundary condition that ${\psi\left(x\right)}$ is continuous at ${x=0}$, which gives us

$\displaystyle Ai\left(-\frac{\alpha E}{mg}\right)=0 \ \ \ \ \ (14)$

To solve this, we need the zeroes of ${Ai\left(x\right)}$ which can be found from tables, or by using Maple’s fsolve method on the AiryAi function. For ${g=9.8\mbox{ m s}^{-2}}$ and a mass of ${0.1\mbox{ kg}}$we have

 $\displaystyle \alpha$ $\displaystyle =$ $\displaystyle 2.602\times10^{22}\mbox{ m}^{-1}\ \ \ \ \ (15)$ $\displaystyle E_{1}$ $\displaystyle =$ $\displaystyle 8.805\times10^{-23}\mbox{ J}\ \ \ \ \ (16)$ $\displaystyle E_{2}$ $\displaystyle =$ $\displaystyle 1.539\times10^{-22}\mbox{ J}\ \ \ \ \ (17)$ $\displaystyle E_{3}$ $\displaystyle =$ $\displaystyle 2.079\times10^{-22}\mbox{ J}\ \ \ \ \ (18)$ $\displaystyle E_{4}$ $\displaystyle =$ $\displaystyle 2.556\times10^{-22}\mbox{ J} \ \ \ \ \ (19)$

If the particle is an electron, ${m=9.11\times10^{-31}\mbox{ kg}}$ and we get

 $\displaystyle \alpha$ $\displaystyle =$ $\displaystyle 1135\mbox{ m}^{-1}\ \ \ \ \ (20)$ $\displaystyle E_{1}$ $\displaystyle =$ $\displaystyle 1.839\times10^{-32}\mbox{ J}\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 1.148\times10^{-13}\mbox{ eV} \ \ \ \ \ (22)$

Using the virial theorem we can find the electron’s average height above the floor in its ground state:

 $\displaystyle 2\left\langle T\right\rangle$ $\displaystyle =$ $\displaystyle \left\langle x\frac{dV}{dx}\right\rangle \ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle mg\left\langle x\right\rangle \ \ \ \ \ (24)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left\langle V\right\rangle \ \ \ \ \ (25)$ $\displaystyle E_{1}$ $\displaystyle =$ $\displaystyle \left\langle T\right\rangle +\left\langle V\right\rangle \ \ \ \ \ (26)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3}{2}\left\langle V\right\rangle \ \ \ \ \ (27)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3}{2}mg\left\langle x\right\rangle \ \ \ \ \ (28)$ $\displaystyle \left\langle x\right\rangle$ $\displaystyle =$ $\displaystyle \frac{2E_{1}}{3mg}\ \ \ \ \ (29)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 0.00137\mbox{ m}\ \ \ \ \ (30)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 1.37\mbox{ mm} \ \ \ \ \ (31)$