Half life of a beer can

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

Whe we analyzed the WKB approximation for a particle tunneling through a barrier we came up with a formula for the transmission probability {T}:

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


where

\displaystyle  \gamma=\frac{1}{\hbar}\int_{-a}^{a}\left|p\left(x\right)\right|dx \ \ \ \ \ (2)

and the barrier extends from {x=-a} to {x=+a}. We can apply this formula to get an idea of why quantum tunneling isn’t a common phenomenon in the macroscopic world.

Suppose we have a can of beer with a mass of {m=0.5\mbox{ kg}}, a diameter of {D=0.06\mbox{ m}} and a height of {h=0.16\mbox{ m}}. If {x} is the height of the centre of mass of the can above its height when the can is standing upright (that is, at {x=h/2}), then the potential energy of the can is

\displaystyle  V\left(x\right)=mgx \ \ \ \ \ (3)

If the total energy of the can is {E=0} then

\displaystyle  \left|p\left(x\right)\right|=\sqrt{2m\left(mgx\right)}=m\sqrt{2gx} \ \ \ \ \ (4)

The can is stable in two orientations: on its end (the way it normally sits on a table) and on its side. There is a potential barrier between these two states. If the can starts in an upright position and is tilted, it will spontaneously fall over onto its side when its centre of mass is no longer within the diameter of the can. That is, with the centre of mass at the centre of the can, this happens when

\displaystyle  x_{0}=\sqrt{\left(D/2\right)^{2}+\left(h/2\right)^{2}}-h/2 \ \ \ \ \ (5)

The barrier through which we wish to tunnel thus extends from {x=0} (can standing upright) to {x=x_{0}} (can about to tip over). We get

\displaystyle   \gamma \displaystyle  = \displaystyle  \frac{m\sqrt{2g}}{\hbar}\int_{0}^{x_{0}}\sqrt{x}dx\ \ \ \ \ (6)
\displaystyle  \displaystyle  = \displaystyle  \frac{2m\sqrt{2g}}{3\hbar}x_{0}^{3/2}\ \ \ \ \ (7)
\displaystyle  \displaystyle  = \displaystyle  \frac{m\sqrt{g}}{3\hbar}\left(\sqrt{D^{2}+h^{2}}-h\right)^{3/2} \ \ \ \ \ (8)

With the values above we get

\displaystyle  \gamma=5.64\times10^{30} \ \ \ \ \ (9)

How long is the ‘half-life’ of a beer can? In our study of alpha decay we got an estimate of the half-life of a particle as (equation 34 there)

\displaystyle  \tau\approx\frac{2r_{1}}{v}e^{2\gamma} \ \ \ \ \ (10)

where {r_{1}} is the distance the particle (or beer can, in this case) must travel to reach the tipping point and {v} is the speed of the particle (can). Treating the can as a ‘particle’, we can its average speed due to thermal motion at room temperature (293 K) as

\displaystyle  v=\sqrt{\frac{k_{B}T}{m}}=8.992\times10^{-11}\mbox{ m s}^{-1} \ \ \ \ \ (11)

The critical distance is

\displaystyle  r_{1}=x_{0}=0.00544\mbox{ m} \ \ \ \ \ (12)

(Recall this is the upward distance the centre of mass must move to reach the tipping point; it’s not the horizontal distance the centre of mass moves. That’s why it’s only about 5 mm.)

Plugging in the numbers, we get

\displaystyle  \tau\approx3.306\times10^{-8}e^{10^{31}}\mbox{ s} \ \ \ \ \ (13)

That’s a very big number. We can get it in powers of 10 using base-10 logarithms:

\displaystyle  \log_{10}\tau\approx10^{31}\log_{10}e=5\times10^{30} \ \ \ \ \ (14)

So the half-life is

\displaystyle  \tau\approx10^{5\times10^{30}}\mbox{ s} \ \ \ \ \ (15)

Converting to years won’t make much difference, since it will subtract only around 7 from the exponent. The age of the universe is currently estimated to be around {10^{10}} years, or around {10^{17}} seconds. We can see from this just how unlikely quantum events are in the macroscopic world.

WKB approximation of double-well potential: wave functions

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

Continuing our application of the WKB approximation to the problem of a double potential well, we can now look at determining the allowed energies for bound states. Since the potential {V\left(x\right)} is even, the wave function is a linear combination of even and odd functions. We worked out the WKB wave functions for {x>0} and they are

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


where

\displaystyle  \theta\equiv\frac{1}{\hbar}\int_{x_{1}}^{x_{2}}p\left(x'\right)dx' \ \ \ \ \ (2)

The odd extension of 1 is thus

\displaystyle  \psi\left(-x\right)=-\psi\left(x\right) \ \ \ \ \ (3)

for {x\ge0}. For this case, we must have {\psi\left(0\right)=}0, so we have

\displaystyle   \psi\left(0\right) \displaystyle  = \displaystyle  \frac{D}{\sqrt{\left|p\left(0\right)\right|}}\left[2\cos\theta e^{\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}+\sin\theta e^{-\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}\right]\ \ \ \ \ (4)
\displaystyle  \displaystyle  = \displaystyle  \frac{D}{\sqrt{\left|p\left(0\right)\right|}}e^{-\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}\left[2\cos\theta e^{\phi}+\sin\theta\right]\ \ \ \ \ (5)
\displaystyle  \phi \displaystyle  \equiv \displaystyle  \frac{2}{\hbar}\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'=\frac{1}{\hbar}\int_{-x_{1}}^{x_{1}}\left|p\left(x'\right)\right|dx' \ \ \ \ \ (6)

where the last line follows from the fact that if {V\left(x\right)} is even, then {p\left(x\right)=\sqrt{2m\left(E-V\left(x\right)\right)}} is also even. Setting 5 to zero gives the condition

\displaystyle  \tan\theta=-2e^{\phi} \ \ \ \ \ (7)

Since both {\phi} and {\theta} depend on {p\left(x\right)}, they both contain the energy {E}, so this condition imposes constraints on {E}.

For the even extension of the WKB function, we have

\displaystyle   \psi\left(-x\right) \displaystyle  = \displaystyle  \psi\left(x\right)\ \ \ \ \ (8)
\displaystyle  \psi'\left(0\right) \displaystyle  = \displaystyle  0 \ \ \ \ \ (9)

The latter condition gives us

\displaystyle  \psi'\left(x\right)=-\frac{1}{2}\frac{D}{\left|p\left(x\right)\right|}p'\left(x\right)\psi\left(x\right)+\frac{D}{\hbar}\sqrt{\left|p\left(x\right)\right|}\left[-2\cos\theta e^{\int_{x}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}+\sin\theta e^{-\int_{x}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}\right] \ \ \ \ \ (10)

Because {p\left(x\right)} is even, {p'\left(0\right)=0} so

\displaystyle   \psi'\left(0\right) \displaystyle  = \displaystyle  \frac{D}{\hbar}\sqrt{\left|p\left(0\right)\right|}\left[-2\cos\theta e^{\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}+\sin\theta e^{-\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}\right]\ \ \ \ \ (11)
\displaystyle  \displaystyle  = \displaystyle  \frac{D}{\hbar}\sqrt{\left|p\left(0\right)\right|}e^{-\int_{0}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}\left[-2\cos\theta e^{\phi}+\sin\theta\right] \ \ \ \ \ (12)

Setting this to zero gives the other energy quantization condition

\displaystyle  \tan\theta=2e^{\phi} \ \ \ \ \ (13)

so the combined conditions are

\displaystyle  \boxed{\tan\theta=\pm2e^{\phi}} \ \ \ \ \ (14)

If the central part of the potential (between {x=-x_{1}} and {x=+x_{1}}) is high and/or broad then from 6, {\phi} will become large, since {p\left(x\right)} depends on {V\left(x\right)}. In that case we see from 14 that {\tan\theta\rightarrow\pm\infty} so {\theta\rightarrow\left(n+\frac{1}{2}\right)\frac{\pi}{2}} for some integer {n}. Rewriting 14 we get

\displaystyle   \theta \displaystyle  = \displaystyle  \mbox{arccot}\left(\pm\frac{1}{2}e^{-\phi}\right)\ \ \ \ \ (15)
\displaystyle  \displaystyle  = \displaystyle  \left(n+\frac{1}{2}\right)\pi\mp\frac{1}{2}e^{-\phi}+\mathcal{O}\left(e^{-3\phi}\right) \ \ \ \ \ (16)

Now suppose we give the potential a specific formula:

\displaystyle  V\left(x\right)=\begin{cases} \frac{1}{2}m\omega^{2}\left(x+a\right)^{2} & x<0\\ \frac{1}{2}m\omega^{2}\left(x-a\right)^{2} & x>0 \end{cases} \ \ \ \ \ (17)

That is, we’re dealing with two linked harmonic oscillator potentials. The turning points are found from the condition {p\left(x\right)=0} and are, for {x>0}:

\displaystyle   x_{1} \displaystyle  = \displaystyle  -\sqrt{\frac{2E}{m\omega^{2}}}+a\ \ \ \ \ (18)
\displaystyle  x_{2} \displaystyle  = \displaystyle  \sqrt{\frac{2E}{m\omega^{2}}}+a \ \ \ \ \ (19)

We can now work out {\theta} from 2 by first working out the integral using Maple

\displaystyle  \frac{1}{\hbar}\int p\left(x\right)dx=-\frac{m^{3/2}\omega^{2}}{2\hbar}\left(a-x\right)\left(-x^{2}+2ax+\frac{2E}{m\omega^{2}}-a^{2}\right)^{1/2}-\frac{E}{\hbar\omega}\sin^{-1}\left(\sqrt{\frac{m\omega^{2}}{2E}}\left(a-x\right)\right) \ \ \ \ \ (20)

The first term is zero at both {x_{1}} and {x_{2}}, and the argument of the arcsine is {-1} at {x_{2}} and {+1} at {x_{1}} so

\displaystyle   \theta \displaystyle  = \displaystyle  \frac{1}{\hbar}\int_{x_{1}}^{x_{2}}p\left(x\right)dx\ \ \ \ \ (21)
\displaystyle  \displaystyle  = \displaystyle  \frac{\pi E}{\hbar\omega} \ \ \ \ \ (22)

Using the approximation 16 we get

\displaystyle  E_{n}^{\pm}\approx\left(n+\frac{1}{2}\right)\hbar\omega\mp\frac{\hbar\omega}{2\pi}e^{-\phi} \ \ \ \ \ (23)

The + energy is for the even wave function, and corresponds to the minus sign on the RHS, so the energies of particles in the even state {\psi_{n}^{+}} are slightly lower than those in the odd state {\psi_{n}^{-}}.

We can get a full time-dependent (approximate) wave function for a particle that starts out in some linear combination of {\psi_{n}^{+}} and {\psi_{n}^{-}} by using the usual technique. Suppose we start the particle out in the following state:

\displaystyle  \Psi\left(x,0\right)=\frac{1}{\sqrt{2}}\left(\psi_{n}^{+}+\psi_{n}^{-}\right) \ \ \ \ \ (24)

This particle is entirely within the well in the region {x>0}, since for {x<0}, {\psi_{n}^{+}\left(x\right)=-\psi_{n}^{-}\left(x\right)} and for {x>0}, {\psi_{n}^{+}\left(x\right)=\psi_{n}^{-}\left(x\right)} so the total wave function is zero for {x<0}. Then the full time-dependent wave function is, using 23

\displaystyle   \Psi\left(x,t\right) \displaystyle  = \displaystyle  \frac{1}{\sqrt{2}}\left(\psi_{n}^{+}e^{-iE_{n}^{+}t/\hbar}+\psi_{n}^{-}e^{-iE_{n}^{-}t/\hbar}\right)\ \ \ \ \ (25)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{\sqrt{2}}e^{-i\left(n+\frac{1}{2}\right)\omega t}\left(\psi_{n}^{+}e^{i\omega te^{-\phi}/2\pi}+\psi_{n}^{-}e^{-i\omega te^{-\phi}/2\pi}\right) \ \ \ \ \ (26)

The probability density is {\left|\Psi\left(x,t\right)\right|^{2}}, which can be calculated fairly easily since {\psi_{n}^{+}} and {\psi_{n}^{-}} are both real functions, as we can see from 1. We get

\displaystyle   \left|\Psi\left(x,t\right)\right|^{2} \displaystyle  = \displaystyle  \frac{1}{2}\left[\left|\psi_{n}^{+}\right|^{2}+\left|\psi_{n}^{-}\right|^{2}+\psi_{n}^{+}\psi_{n}^{-}\left(e^{i\omega te^{-\phi}/\pi}+e^{-i\omega te^{-\phi}/\pi}\right)\right]\ \ \ \ \ (27)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{2}\left(\left|\psi_{n}^{+}\right|^{2}+\left|\psi_{n}^{-}\right|^{2}\right)+\psi_{n}^{+}\psi_{n}^{-}\cos\frac{\omega te^{-\phi}}{\pi} \ \ \ \ \ (28)

If {x>0}, {\psi_{n}^{+}\psi_{n}^{-}=\left|\psi_{n}^{+}\right|^{2}} while if {x<0}, {\psi_{n}^{+}\psi_{n}^{-}=-\left|\psi_{n}^{+}\right|^{2}} so as time progresses, the probability that the particle is one well or the other oscillates between the two wells, with a period given by

\displaystyle   \frac{\omega e^{-\phi}}{\pi} \displaystyle  = \displaystyle  \frac{2\pi}{\tau}\ \ \ \ \ (29)
\displaystyle  \tau \displaystyle  = \displaystyle  \frac{2\pi^{2}e^{\phi}}{\omega} \ \ \ \ \ (30)

Finally, we can calculate {\phi} for the double harmonic oscillator potential, using 6. We get

\displaystyle   \phi \displaystyle  = \displaystyle  \frac{2}{\hbar}\int_{0}^{x_{1}}\sqrt{2m\left(\frac{1}{2}m\omega^{2}\left(x-a\right)^{2}-E\right)}dx\ \ \ \ \ (31)
\displaystyle  \displaystyle  = \displaystyle  \frac{2\sqrt{2mE}}{\hbar}\int_{0}^{x_{1}}\sqrt{\left(\frac{m\omega^{2}}{2E}\left(x-a\right)^{2}-1\right)}dx \ \ \ \ \ (32)

To transform this to a form that Maple can handle, we use the substitution

\displaystyle   u \displaystyle  = \displaystyle  \sqrt{\frac{m\omega^{2}}{2E}}\left(x-a\right)\ \ \ \ \ (33)
\displaystyle  du \displaystyle  = \displaystyle  \sqrt{\frac{m\omega^{2}}{2E}}dx \ \ \ \ \ (34)

The limits transform as

\displaystyle   x \displaystyle  = \displaystyle  0\rightarrow u=-a\sqrt{\frac{m\omega^{2}}{2E}}\ \ \ \ \ (35)
\displaystyle  x \displaystyle  = \displaystyle  x_{1}=-\sqrt{\frac{2E}{m\omega^{2}}}+a\rightarrow u=-1 \ \ \ \ \ (36)

So

\displaystyle   \phi \displaystyle  = \displaystyle  \frac{2\sqrt{2mE}}{\hbar}\sqrt{\frac{2E}{m\omega^{2}}}\int_{-a\sqrt{\frac{m\omega^{2}}{2E}}}^{-1}\sqrt{u^{2}-1}du\ \ \ \ \ (37)
\displaystyle  \displaystyle  = \displaystyle  \frac{4E}{\hbar\omega}\int_{1}^{a\sqrt{\frac{m\omega^{2}}{2E}}}\sqrt{u^{2}-1}du \ \ \ \ \ (38)

where the last step uses the fact that the integrand is even. Since {V\left(0\right)=\frac{1}{2}m\omega^{2}a^{2}\equiv V_{0}} we can write this as

\displaystyle   \phi \displaystyle  = \displaystyle  \frac{4E}{\hbar\omega}\int_{1}^{a\sqrt{V_{0}/E}}\sqrt{u^{2}-1}du\ \ \ \ \ (39)
\displaystyle  \displaystyle  = \displaystyle  \frac{2E}{\hbar\omega}\left.\left(u\sqrt{u^{2}-1}-\ln\left(u+\sqrt{u^{2}-1}\right)\right)\right|_{1}^{\sqrt{V_{0}/E}}\ \ \ \ \ (40)
\displaystyle  \displaystyle  = \displaystyle  \frac{2E}{\hbar\omega}\left[\sqrt{\frac{V_{0}}{E}}\sqrt{\frac{V_{0}}{E}-1}-\ln\left(\sqrt{\frac{V_{0}}{E}}+\sqrt{\frac{V_{0}}{E}-1}\right)\right] \ \ \ \ \ (41)

For a high central barrier, {V_{0}\gg E} and we can approximate this by

\displaystyle   \phi \displaystyle  \approx \displaystyle  \frac{2E}{\hbar\omega}\left[\frac{V_{0}}{E}-\ln\left(2\sqrt{\frac{V_{0}}{E}}\right)\right]\ \ \ \ \ (42)
\displaystyle  \displaystyle  \approx \displaystyle  \frac{2V_{0}}{\hbar\omega}=\frac{m\omega a^{2}}{\hbar} \ \ \ \ \ (43)

where we’ve dropped the logarithm term as it’s much smaller than {V_{0}/E}.

WKB approximation of a double potential well: turning points

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

In this post, we’ll start applying the WKB approximation to the problem of a double potential well, in which the potential is an even function with two minima, at {\pm x_{0}} and a finite maximum {V_{0}} at {x=0}. The potential will tend to {+\infty} as {x\rightarrow\pm\infty}, and we’ll be looking for bound states with a total energy {E<V_{0}}, so classically, the particle would be confined to one of the two wells. To make things definite, we’ll look at the wave function in the right-hand well, where the turning points are at {x_{1}} (where {V'\left(x_{1}\right)<0}) and {x_{2}} (where {V'\left(x_{2}\right)>0}).

The situation at turning point {x_{2}} is the same as with the single-well potential, so we get

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


where {D} is a normalization constant.

To get the wave function for {0<x<x_{1}} we need to use the patching function to connect the WKB functions on either side of {x=x_{1}}. We’ve done a similar analysis for the downward sloping turning point in a single well potential, but in that case we assumed that {V\left(x\right)\rightarrow\infty} as {x\rightarrow-\infty} so we could throw away the positive exponential term in the wave function. In this case, we can’t do that since the potential is finite to the left of {x=0}, at least up to the point where it decreases and goes below {E} in the left-hand well.

The plan is essentially the same as that used in analyzing the finite barrier with sloping sides. First, we linearize the potential near the point {x=x_{1}}, then we work out the patching function in the region {x>x_{1}} and relate it to the WKB function in that region. Since we already know the WKB function for {x_{1}<x<x_{2}}, we make the WKB function from the left turning point match that given in 1. Finally, we use the patching function to find the WKB function for {0<x<x_{1}}.

As usual, we’ll do the analysis by shifting {x_{1}} to the origin. Then the WKB functions near {x=0} are

\displaystyle  \psi=\begin{cases} \frac{1}{\sqrt{\left|p\right|}}\left[Ge^{\int_{x}^{0}\left|p\left(x'\right)\right|dx'/\hbar}+Fe^{-\int_{x}^{0}\left|p\left(x'\right)\right|dx'/\hbar}\right] & x<0\\ \frac{1}{\sqrt{p}}\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} \ \ \ \ \ (2)

Linearizing the potential, we have

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

As before, solving the Schrödinger equation for this linearized potential gives us the patching wave function {\psi_{p}}:

\displaystyle   \psi_{p}\left(x\right) \displaystyle  = \displaystyle  aAi\left(\alpha x\right)+bBi\left(\alpha x\right)\ \ \ \ \ (4)
\displaystyle  \alpha \displaystyle  \equiv \displaystyle  \left(\frac{2mV'\left(0\right)}{\hbar^{2}}\right)^{1/3} \ \ \ \ \ (5)

In the overlap region to the right of {x_{1}}, {x>0} so

\displaystyle  p\left(x\right)=\sqrt{2m\left(E-V\left(x\right)\right)}=\hbar\sqrt{-\alpha^{3}x} \ \ \ \ \ (6)

Note that {p} is real, since {\alpha<0} (due to {V'\left(0\right)<0}) and {x>0}. Therefore

\displaystyle  \int_{0}^{x}p\left(x'\right)dx'=\frac{2}{3}\hbar\left(-\alpha x\right)^{3/2} \ \ \ \ \ (7)

and the WKB function from 2 for {x>0} is

\displaystyle   \psi\left(x\right) \displaystyle  = \displaystyle  \frac{1}{\sqrt{\hbar}\left(-\alpha\right)^{3/4}x^{1/4}}\left(Be^{iq}+Ce^{-iq}\right)\ \ \ \ \ (8)
\displaystyle  q \displaystyle  \equiv \displaystyle  \frac{1}{\hbar}\int_{0}^{x}p\left(x'\right)dx'=\frac{2}{3}\left(-\alpha x\right)^{3/2} \ \ \ \ \ (9)

Now for the patching function {\psi_{p}} for {x>0} from 4, we see that the argument is negative since {\alpha<0}, so we can apply the large negative asymptotic forms for the Airy functions:

\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)

This gives us

\displaystyle   \psi_{p}\left(x\right) \displaystyle  \approx \displaystyle  \frac{1}{\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\left[a\sin\left(q+\frac{\pi}{4}\right)+b\cos\left(q+\frac{\pi}{4}\right)\right]\ \ \ \ \ (12)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{\sqrt{\pi}\left(-\alpha x\right)^{1/4}}\left[\frac{a}{2i}\left(e^{iq+i\pi/4}-e^{-iq-i\pi/4}\right)+\frac{b}{2}\left(e^{iq+i\pi/4}+e^{-iq-i\pi/4}\right)\right] \ \ \ \ \ (13)

Equating coefficients of {e^{\pm iq}} in this equation and the WKB function 8 we get

\displaystyle   a \displaystyle  = \displaystyle  i\sqrt{\frac{\pi}{-\hbar\alpha}}\left(Be^{-i\pi/4}-Ce^{i\pi/4}\right)\ \ \ \ \ (14)
\displaystyle  b \displaystyle  = \displaystyle  \sqrt{\frac{\pi}{-\hbar\alpha}}\left(Be^{-i\pi/4}+Ce^{i\pi/4}\right) \ \ \ \ \ (15)

The WKB function 8 must be the same as 1 in the region {x_{1}<x<x_{2}}. To use this fact, we note that

\displaystyle  \int_{x}^{x_{2}}p\left(x'\right)dx'=\int_{x_{1}}^{x_{2}}p\left(x'\right)dx'-\int_{x_{1}}^{x}p\left(x'\right)dx' \ \ \ \ \ (16)

Applying this to 1 we get

\displaystyle   \psi\left(x\right) \displaystyle  \approx \displaystyle  \frac{2D}{\sqrt{p\left(x\right)}}\sin\left[\frac{1}{\hbar}\int_{x_{1}}^{x_{2}}p\left(x'\right)dx'-\frac{1}{\hbar}\int_{x_{1}}^{x}p\left(x'\right)dx'+\frac{\pi}{4}\right]\ \ \ \ \ (17)
\displaystyle  \displaystyle  = \displaystyle  \frac{2D}{\sqrt{p\left(x\right)}}\sin\left(\theta-q+\frac{\pi}{4}\right)\ \ \ \ \ (18)
\displaystyle  \displaystyle  = \displaystyle  \frac{2D}{\sqrt{p\left(x\right)}}\frac{1}{2i}\left(e^{i\theta-iq+i\pi/4}-e^{-i\theta+iq-i\pi/4}\right) \ \ \ \ \ (19)

where

\displaystyle  \theta\equiv\frac{1}{\hbar}\int_{x_{1}}^{x_{2}}p\left(x'\right)dx' \ \ \ \ \ (20)

This function must be the same as 8 so we can get the constants {B} and {C} in terms of {D}:

\displaystyle   B \displaystyle  = \displaystyle  -\frac{D}{i}e^{-i\theta-i\pi/4}\ \ \ \ \ (21)
\displaystyle  C \displaystyle  = \displaystyle  \frac{D}{i}e^{i\theta+i\pi/4} \ \ \ \ \ (22)

We can insert these into 14 and 15 to get

\displaystyle   a \displaystyle  = \displaystyle  -\sqrt{\frac{\pi}{-\hbar\alpha}}D\left(e^{-i\theta-i\pi/2}+e^{i\theta+i\pi/2}\right)\ \ \ \ \ (23)
\displaystyle  \displaystyle  = \displaystyle  -2\sqrt{\frac{\pi}{-\hbar\alpha}}D\cos\left(\theta+\frac{\pi}{2}\right)\ \ \ \ \ (24)
\displaystyle  \displaystyle  = \displaystyle  2\sqrt{\frac{\pi}{-\hbar\alpha}}D\sin\theta\ \ \ \ \ (25)
\displaystyle  b \displaystyle  = \displaystyle  \sqrt{\frac{\pi}{-\hbar\alpha}}\frac{D}{i}\left(-e^{-i\theta-i\pi/2}+e^{i\theta+i\pi/2}\right)\ \ \ \ \ (26)
\displaystyle  \displaystyle  = \displaystyle  2\sqrt{\frac{\pi}{-\hbar\alpha}}D\sin\left(\theta+\frac{\pi}{2}\right)\ \ \ \ \ (27)
\displaystyle  \displaystyle  = \displaystyle  2\sqrt{\frac{\pi}{-\hbar\alpha}}D\cos\theta \ \ \ \ \ (28)

Now we can look at the other side of the turning point, where {x<x_{1}} (or {x<0} in our calculations). First, we note that, from 6

\displaystyle  \left|p\left(x\right)\right|=\left|\hbar\sqrt{-\alpha^{3}x}\right|=\hbar\sqrt{\alpha^{3}x} \ \ \ \ \ (29)

since {\alpha<0} and {x<0}. Also, since {\left|p\left(x\right)\right|} is an even function and {x<0}

\displaystyle  \int_{x}^{0}\left|p\left(x'\right)\right|dx'=\int_{0}^{-x}p\left(x'\right)dx'=\hbar q \ \ \ \ \ (30)

Therefore, from 2 with {x<0} we get

\displaystyle  \psi\left(x\right)=\frac{1}{\sqrt{\hbar}\left(-\alpha\right)^{3/4}\left(-x\right)^{1/4}}\left(Ge^{q}+Fe^{-q}\right) \ \ \ \ \ (31)

Using the same patching function as before but now taking its large positive asymptotic form (since {\alpha<0} and {x<0}, so {\alpha x>0}), we get

\displaystyle  \psi_{p}\left(x\right)\approx\frac{1}{\sqrt{\pi}\left(\alpha x\right)^{1/4}}\left(\frac{a}{2}e^{-q}+be^{q}\right) \ \ \ \ \ (32)

Comparing coefficients with the WKB form:

\displaystyle   \frac{a}{2} \displaystyle  = \displaystyle  \sqrt{\frac{\pi}{-\hbar\alpha}}F\ \ \ \ \ (33)
\displaystyle  b \displaystyle  = \displaystyle  \sqrt{\frac{\pi}{-\hbar\alpha}}G \ \ \ \ \ (34)

Using 25 and 28 we can write {F} and {G} in terms of {D}:

\displaystyle   F \displaystyle  = \displaystyle  D\sin\theta\ \ \ \ \ (35)
\displaystyle  G \displaystyle  = \displaystyle  2D\cos\theta \ \ \ \ \ (36)

Substituting back into the first of 2 and restoring the turning point from {x=0} to {x=x_{1}} we get for {0<x<x_{1}}:

\displaystyle  \psi\left(x\right)\approx\frac{D}{\sqrt{\left|p\left(x\right)\right|}}\left[2\cos\theta e^{\int_{x}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}+\sin\theta e^{-\int_{x}^{x_{1}}\left|p\left(x'\right)\right|dx'/\hbar}\right] \ \ \ \ \ (37)

This, together with 1, gives us the complete WKB wave function for {x>0}.

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] & r<r_{2}\\ \frac{D}{\sqrt{\left|p\left(r\right)\right|}}\exp\left[-\int_{r_{2}}^{r}\left|p\left(r'\right)\right|\; dr'/\hbar\right] & r>r_{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<V\left(x\right)} 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<V\left(x\right)} 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] & x<x_{1}\\ \frac{1}{\sqrt{\left|p\left(x\right)\right|}}\left[Ce^{\int_{x_{1}}^{x}\left|p\left(x'\right)\right|dx'/\hbar}+De^{-\int_{x_{1}}^{x}\left|p\left(x'\right)\right|dx'/\hbar}\right] & x_{1}<x<x_{2}\\ \frac{1}{\sqrt{p\left(x\right)}}\left[Fe^{i\int_{x_{2}}^{x}p\left(x'\right)dx'/\hbar}\right] & x>x_{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] & x<x_{2}\\ \frac{D}{\sqrt{\left|p\left(x\right)\right|}}\exp\left[-\int_{x_{2}}^{x}\left|p\left(x'\right)\right|\; dx'/\hbar\right] & x>x_{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<x_{1} \end{cases} \ \ \ \ \ (2)

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}<x<x_{2}} 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.

Follow

Get every new post delivered to your Inbox.

Join 351 other followers