Author Archives: growescience

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.

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<V\left(x\right)}, 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<V\left(x\right)} for {x<x_{1}} 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)

For the WKB function, we start with 5

\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<V\left(x\right)}, 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<V\left(x\right)} 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<V\left(x\right)} 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)

For the WKB function, we start with 1

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

Alpha decay using the WKB approximation

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

One of the early applications of the WKB approximation was George Gamov’s 1928 theory of alpha decay. In a large nucleus, the nucleons (protons and neutrons) are held together by the strong nuclear force which at short range is stronger than the electric repulsion between the protons. The strong nuclear force is very short range, however, so if some nucleons can tunnel through the potential barrier, the electric force rapidly takes over resulting in the nucleons being ejected from the nucleus.

A common mode of decay is the emission of an alpha particle, consisting of 2 neutrons and 2 protons (a helium-4 nucleus). Gamow’s model proposed that the nuclear force be modelled as a square well at potential {-V_{0}} and width {r_{1}} followed by the Coulomb repulsive potential for {r>r_{1}}. That is

\displaystyle  V\left(x\right)=\begin{cases} -V_{0} & 0<r<r_{1}\\ \frac{2Ze^{2}}{4\pi\epsilon_{0}r} & r>r_{1} \end{cases} \ \ \ \ \ (1)

where {Z} is the number of protons in the nucleus of the atom remaining after the emission of the alpha particle. The 2 in the numerator is the number of protons in the alpha particle.

If the alpha particle has an energy {E} such that {0<E<\frac{2Ze^{2}}{4\pi\epsilon_{0}r_{1}^{2}}} then its energy is below the maximum of {V} which occurs at {r=r_{1}}, so in order for it to escape, it must tunnel through the potential in the region {r_{1}<r<r_{2}} where the outer radius {r_{2}} is determined by

\displaystyle   E \displaystyle  = \displaystyle  \frac{2Ze^{2}}{4\pi\epsilon_{0}r_{2}}\ \ \ \ \ (2)
\displaystyle  r_{2} \displaystyle  = \displaystyle  \frac{2Ze^{2}}{4\pi\epsilon_{0}E} \ \ \ \ \ (3)

We saw that applying the WKB approximation to tunneling gave us a transmission probability of

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


where

\displaystyle   \gamma \displaystyle  = \displaystyle  \frac{1}{\hbar}\int_{r_{1}}^{r_{2}}\sqrt{2m\left(V\left(r\right)-E\right)}dr\ \ \ \ \ (5)
\displaystyle  \displaystyle  = \displaystyle  \frac{\sqrt{2m}}{\hbar}\int_{r_{1}}^{r_{2}}\sqrt{\frac{2Ze^{2}}{4\pi\epsilon_{0}}}\sqrt{\frac{1}{r}-\frac{1}{r_{2}}}dr\ \ \ \ \ (6)
\displaystyle  \displaystyle  = \displaystyle  \frac{\sqrt{2m}}{\hbar}\sqrt{\frac{2Ze^{2}}{4\pi\epsilon_{0}r_{2}}}\int_{r_{1}}^{r_{2}}\sqrt{\frac{r_{2}}{r}-1}dr\ \ \ \ \ (7)
\displaystyle  \displaystyle  = \displaystyle  \frac{\sqrt{2mE}}{\hbar}\int_{r_{1}}^{r_{2}}\sqrt{\frac{r_{2}}{r}-1}dr \ \ \ \ \ (8)

Doing the integral with Maple gives

\displaystyle  \gamma=\frac{\sqrt{2mE}}{\hbar}\left[\frac{\pi}{4}r_{2}-\sqrt{r_{1}\left(r_{2}-r_{1}\right)}-\frac{r_{2}}{2}\arcsin\left(\frac{2r_{1}}{r_{2}}-1\right)\right] \ \ \ \ \ (9)

We can write this in terms of {\rho\equiv r_{1}/r_{2}}:

\displaystyle  \gamma=\frac{\sqrt{2mE}}{\hbar}\frac{r_{2}}{4}\left(\pi-4\sqrt{\rho\left(1-\rho\right)}-2\arcsin\left(2\rho-1\right)\right) \ \ \ \ \ (10)

If we assume that {r_{1}\ll r_{2}} (that is, the range of the nuclear force is much less than the point where the alpha particle breaks free of the nucleus, which is usually the case), then a series expansion about {\rho=0} gives

\displaystyle   \gamma \displaystyle  = \displaystyle  \frac{\sqrt{2mE}}{\hbar}\left[\frac{\pi}{2}r_{2}-2r_{2}\sqrt{\rho}+\frac{r_{2}}{3}\rho^{3/2}+...\right]\ \ \ \ \ (11)
\displaystyle  \displaystyle  = \displaystyle  \frac{\sqrt{2mE}}{\hbar}\left[\frac{\pi}{2}r_{2}-2\sqrt{r_{1}r_{2}}+\frac{r_{1}}{3}\sqrt{\frac{r_{1}}{r_{2}}}+...\right] \ \ \ \ \ (12)

Taking just the first two terms, we have, using 3

\displaystyle   \gamma \displaystyle  \cong \displaystyle  \frac{\sqrt{2mE}}{\hbar}\frac{Ze^{2}}{4\epsilon_{0}E}-\frac{2\sqrt{2mE}}{\hbar}\sqrt{\frac{2Ze^{2}}{4\pi\epsilon_{0}E}}\sqrt{r_{1}}\ \ \ \ \ (13)
\displaystyle  \displaystyle  = \displaystyle  \frac{e^{2}\pi\sqrt{2m}}{4\pi\epsilon_{0}\hbar}\frac{Z}{\sqrt{E}}-\frac{4\sqrt{m}}{\hbar}\sqrt{\frac{e^{2}}{4\pi\epsilon_{0}}}\sqrt{Zr_{1}}\ \ \ \ \ (14)
\displaystyle  \displaystyle  = \displaystyle  K_{1}\frac{Z}{\sqrt{E}}-K_{2}\sqrt{Zr_{1}} \ \ \ \ \ (15)

where the constants are

\displaystyle   K_{1} \displaystyle  = \displaystyle  \frac{e^{2}\pi\sqrt{2m}}{4\pi\epsilon_{0}\hbar}\ \ \ \ \ (16)
\displaystyle  K_{2} \displaystyle  = \displaystyle  \frac{4\sqrt{m}}{\hbar}\sqrt{\frac{e^{2}}{4\pi\epsilon_{0}}} \ \ \ \ \ (17)

Plugging in the numbers gives (in SI units; remember that {m} is the mass of the emitted alpha particle):

\displaystyle   \frac{e^{2}}{4\pi\epsilon_{0}} \displaystyle  = \displaystyle  2.307\times10^{-28}\mbox{m}^{3}\mbox{kg s}^{-2}\ \ \ \ \ (18)
\displaystyle  \hbar \displaystyle  = \displaystyle  1.0546\times10^{-34}\mbox{ kg m}^{2}\mbox{ s}^{-1}\ \ \ \ \ (19)
\displaystyle  m \displaystyle  = \displaystyle  6.645\times10^{-27}\mbox{ kg}\ \ \ \ \ (20)
\displaystyle  K_{1} \displaystyle  = \displaystyle  7.923\times10^{-7}\mbox{ kg}^{1/2}\mbox{m s}^{-1}\ \ \ \ \ (21)
\displaystyle  \displaystyle  = \displaystyle  7.923\times10^{-7}\mbox{ J}^{1/2}\ \ \ \ \ (22)
\displaystyle  \displaystyle  = \displaystyle  1.979\mbox{ MeV}^{1/2}\ \ \ \ \ (23)
\displaystyle  K_{2} \displaystyle  = \displaystyle  4.696\times10^{7}\mbox{ m}^{-1/2}\ \ \ \ \ (24)
\displaystyle  \displaystyle  = \displaystyle  1.485\mbox{ fm}^{-1/2} \ \ \ \ \ (25)

where the fm is a fermi or {10^{-15}\mbox{ m}}, so {10^{15}\mbox{ m}^{-1}=1\mbox{ fm}^{-1}} and {\sqrt{10^{15}}\mbox{ m}^{-1/2}=3.16\times10^{7}\mbox{ m}^{-1/2}=1\mbox{ fm}^{-1/2}}.

Example 1 As a specific example, we’ll look at {^{238}\mbox{U}}. When this isotope of uranium emits an alpha particle, the remaining nucleus is {^{234}\mbox{Th}} (thorium). The masses can be found at this site, and we get

\displaystyle   ^{238}\mbox{U}:\;238.050782
\displaystyle  ^{234}\mbox{Th}:\;234.043595
\displaystyle  ^{4}\mbox{He}:\;4.002603 \ \ \ \ \ (26)

The masses are in unified atomic mass units (u), where one u is 1/12 the mass of {^{12}\mbox{C}} or {1.66\times10^{-27}\mbox{ kg}} or {931.5\mbox{ MeV}}.

Returning to {^{238}\mbox{U}} we need values for {r_{1}} and {E}. From experiment, we have

\displaystyle  r_{1}\cong1.07A^{1/3}\mbox{ fm} \ \ \ \ \ (27)

where {A} is the total number of protons and neutrons in the nucleus. From special relativity, we can get the energy of the emitted alpha particle as the difference between the rest masses of the parent nucleus {m_{p}} and those of the residual nucleus {m_{r}} plus alpha particle {m_{\alpha}}:

\displaystyle  E=m_{p}c^{2}-\left(m_{r}c^{2}+m_{\alpha}c^{2}\right) \ \ \ \ \ (28)

From the numbers above, we get

\displaystyle   r_{1} \displaystyle  = \displaystyle  1.07\left(238\right)^{1/3}=6.63\mbox{ fm}\ \ \ \ \ (29)
\displaystyle  E \displaystyle  = \displaystyle  \left(238.050782-234.043595-4.002603\right)\times931.5\ \ \ \ \ (30)
\displaystyle  \displaystyle  = \displaystyle  4.27\mbox{ MeV} \ \ \ \ \ (31)

To calculate {\gamma} we use 15 with {Z=90} for the decay product (thorium), and we get

\displaystyle  \gamma=50.06 \ \ \ \ \ (32)

To get a prediction of the lifetime of a {^{238}\mbox{U}} nucleus, we can use a crude argument to get a rough estimate. Suppose that, before its escape, the alpha particle is moving within the nucleus the some velocity {v}, giving it a kinetic energy of

\displaystyle  E_{k}=\frac{1}{2}m_{\alpha}v^{2} \ \ \ \ \ (33)

If we’re interested only in an order of magnitude estimate, we can ignore the nuclear binding energy {-V_{0}} within the nucleus and take {E\approx E_{k}}. (In reality, {E=E_{k}-V_{0}}, so making this assumption underestimates the velocity.) On average, the alpha particle will collide with the boundary (that is, it will reach {r=r_{1}}) with a frequency of {v/2r_{1}} times per second (since it travels the diameter of the nucleus between collisions), and from 4, the probability of tunnelling through the barrier is {T\approx e^{-2\gamma}} on each collision. Thus the probability of emission per second is about {ve^{-2\gamma}/2r_{1}} making the approximate lifetime {\tau} of the parent nucleus equal to the reciprocal of this or

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

For {^{238}\mbox{U}} we get

\displaystyle   v \displaystyle  = \displaystyle  \sqrt{\frac{2E}{m_{\alpha}}}\ \ \ \ \ (35)
\displaystyle  \displaystyle  = \displaystyle  \sqrt{\frac{2E}{m_{\alpha}c^{2}}}c\ \ \ \ \ (36)
\displaystyle  \displaystyle  = \displaystyle  \sqrt{\frac{2\times4.27}{4.002603\times931.5}}c\ \ \ \ \ (37)
\displaystyle  \displaystyle  = \displaystyle  1.436\times10^{7}\mbox{ m s}^{-1} \ \ \ \ \ (38)

This is about 4% of the speed of light, so we’re justified in using a non-relativistic approximation. The lifetime is thus

\displaystyle  \tau=2.8\times10^{22}\mbox{ s}=8.87\times10^{14}\mbox{ years} \ \ \ \ \ (39)

The actual half-life of {^{238}\mbox{U}} is around {4.5\times10^{9}\mbox{ years}}, so this estimate is pretty wide of the mark, but at least it agrees that the lifetime is very long.

Example 2 We can do the same calculation for {^{212}\mbox{Po}} which decays to {^{208}\mbox{Pb}} and we find the atomic masses are

\displaystyle   ^{212}\mbox{Po}:\;211.988851
\displaystyle  ^{208}\mbox{Pb}:\;207.976635

Plugging in the numbers as before, we get

\displaystyle   r_{1} \displaystyle  = \displaystyle  6.34\mbox{ fm}\ \ \ \ \ (40)
\displaystyle  E \displaystyle  = \displaystyle  8.954\mbox{ MeV}\ \ \ \ \ (41)
\displaystyle  v \displaystyle  = \displaystyle  2.079\times10^{7}\mbox{ m s}^{-1}\ \ \ \ \ (42)
\displaystyle  \gamma \displaystyle  = \displaystyle  20.399\ \ \ \ \ (43)
\displaystyle  \tau \displaystyle  = \displaystyle  3.19\times10^{-4}\mbox{ s} \ \ \ \ \ (44)

The experimental value is {3\times10^{-7}\mbox{ s}} so again the WKB calculation is quite a bit longer, but at least it gives good qualitative agreement. Because the dependence on {\gamma} (and thus on {E}, {Z} and {r_{1}}) is exponential, relatively small changes in these quantities translate into a huge difference in lifetime.

Pulsar planets

References: edX MOOC on Exoplanets; Week 1

An exoplanet is a planet that orbits a star other than the Sun. Recently, many such exoplanets have been discovered. We’ll look at some of the techniques used to discover them.

One of the more unlikely places in which exoplanets have been found is in orbit around pulsars. A pulsar is a neutron star (a star where the matter has become so compressed it is essentially composed of nuclear matter without the empty space that is present in normal atoms) that emits pulses of radiation at periodic intervals. A millisecond pulsar emits pulses hundreds of times a second.

One way of detecting that a planet is orbiting a pulsar is to measure variations in the pulse rate. Suppose a planet is in an orbit that is edge-on as viewed from Earth. Since the planet-pulsar system orbits around its centre of mass, the pulsar’s distance from Earth varies with a period equal to the period of the orbit. When the pulsar is furthest away from Earth, its pulses must cover an extra distance equal to the diameter of the orbit compared to when the pulsar is closest to the Earth, meaning that its pulses are delayed by a time equal to {2r_{*}/c} where {r_{*}} is the radius of the orbit and {c} is the speed of light.

To get more information about the planet, we use a bit of Newtonian physics. In what follows, we assume that the mass {M_{*}} of the pulsar is much larger than {M_{p}}, the mass of the planet. This in turn leads to the assumption that the radius of the planet’s orbit about the centre of mass is much larger the radius of the star’s orbit: {r_{p}\gg r_{*}} and so the distance between the star and planet is {r_{p}+r_{*}\approx r_{p}}.

Now for a bit of force balancing. The centripetal force holding the planet in its orbit is provided by the gravitational force, so in our approximation

\displaystyle  \frac{GM_{p}M_{*}}{r_{p}^{2}}=\frac{M_{p}v_{p}^{2}}{r_{p}} \ \ \ \ \ (1)

where {v_{p}} is the speed of the planet in its orbit. We can’t measure {v_{p}} directly, but we know that (assuming the orbit is circular)

\displaystyle  v_{p}=\frac{2\pi r_{p}}{T} \ \ \ \ \ (2)

where {T} is the period of the orbit, which we can measure by observing the variation in the timing of the pulsar’s pulses. Combining these two equations, we get

\displaystyle  r_{p}=\left(\frac{GM_{*}T^{2}}{4\pi^{2}}\right)^{1/3} \ \ \ \ \ (3)

This is Kepler’s third law of planetary motion: the square of the orbital period is proportional to the cube of the radius (or more generally, semi-major axis for elliptical orbits). To use this formula we do need to know {M_{*}} but this can usually be obtained by other means (from the type of star, for example).

To find the mass {M_{p}} of the planet, we can use the fact that because planet and star both orbit the centre of mass:

\displaystyle  M_{*}r_{*}=M_{p}r_{p} \ \ \ \ \ (4)

We don’t know either {r_{*}} or {M_{p}} at this point, but if the star moves along the line of sight (that is, its orbital plane isn’t perpendicular to the line of sight), there is a Doppler shift of the star’s light (or, more generally, in the radiation emitted by the star). We can measure this and determine the component of the speed of the star along the line of sight. Unless the orbit is edge-on, this speed is less than the actual speed {v_{*}} of the star in its orbit. Suppose the angle between the axis of the orbit and the line of sight from Earth is {i} (the inclination); then the line of sight component of the speed that we measure is {v_{*}\sin i}. Thus we get

\displaystyle   v_{*} \displaystyle  = \displaystyle  \frac{2\pi r_{*}}{T}\ \ \ \ \ (5)
\displaystyle  \displaystyle  = \displaystyle  \frac{2\pi}{T}\frac{M_{p}}{M_{*}}r_{p}\ \ \ \ \ (6)
\displaystyle  v_{*}\sin i \displaystyle  = \displaystyle  \frac{2\pi}{T}\frac{M_{p}}{M_{*}}r_{p}\sin i\ \ \ \ \ (7)
\displaystyle  M_{p}\sin i \displaystyle  = \displaystyle  \frac{M_{*}T}{2\pi r_{p}}\left(v_{*}\sin i\right)\ \ \ \ \ (8)
\displaystyle  M_{p} \displaystyle  \ge \displaystyle  \frac{M_{*}T}{2\pi r_{p}}v_{*measured} \ \ \ \ \ (9)

Thus the best we can do is get a lower limit on the mass of the planet.

Follow

Get every new post delivered to your Inbox.

Join 286 other followers