Jupyter notebook mode

4.3. Connection formulas

Expected prerequisites

Before the start of this lecture, you should be able to:

1. Write down the Schrödinger equation.
2. Do the Taylor expansion of a function and perform simple integrations.

Learning goals

After this lecture, you will be able to:

1. Describe the behaviour of the Airy functions and why they arise in the WBK approximation.
2. Write down the WKB wave function for a given potential in all regions.

4.3.1. On why turning points are problematic

So far, we have only derived the WKB wave function for cases where either $E>V(x)$ for all $x$ or $E<V(x)$ for all $x$. If we have a potential $V(x)$ and an energy $E$ such that both regions $E>V(x)$ and $E < V(x)$ are present, there will be a point $x_t$ where $E = V(x_t)$. In classical mechanics, $x_t$ is the point where the motion of a particle is reversed, and is hence referred to as turning point.

The WKB approximation inevitably breaks down as we approach the turning point $x_t$. On a mathematical level, we see that $p(x)\rightarrow 0$ as $x\rightarrow x_t$, and hence, $x_t$ appears in the middle region as can be observed in the plot below. At this point we have that, $$\psi_{WKB}(x \rightarrow x_t) \sim \frac{1}{\sqrt{p(x\rightarrow x_t)}} \rightarrow \infty.$$ As the result, the wave function becomes infinite and unphysical. The underlying reason for the failure of this model is the following: As we approach $x_t$, the wavelength $\lambda \sim \frac{1}{p(x)}$ becomes longer and longer. Thus, the condition that the potential must be smooth compared to the wavelength can no longer be satisfied and WKB breaks down.

Strategy to deal with turning points in WKB approximation

This is a fundamental problem of WKB - it will always fail at a turning point. We can however solve this problem by:

1. defining three well-behaved regions,
2. solving for the wave functions in each region,
3. and finally matching the solutions.

The regions are indicated in the plot below.

import common
common.configure_plotting()
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(0, np.pi, 100)
y = -np.cos(x)
z = x-np.pi/2

fig, ax = plt.subplots(figsize=(10, 5))
plt.plot(x, y, label=r'$V(x)$', color='black', linewidth=1)
plt.axhline(0, color='r', label=r'$E$')
plt.axvline(1, linestyle='dashed', color='gray')
plt.axvline(2.15, color='gray', linestyle='dashed')
plt.plot(x, z, label='Linear approx.', linewidth=3, color='darkblue')
plt.legend(fontsize=12)
plt.ylim(-1.1, 1.1);
plt.xlim(0, 3.15);
plt.yticks([]);
plt.xticks([]);
plt.xlabel(r'$x_t$', fontsize=15)
plt.text(0.3, 0.1, r'$\Psi_I$', fontsize=18)
plt.text(1.5, 0.8, r'$\Psi_{II}$', fontsize=18)
plt.text(2.7, 0.1, r'$\Psi_{III}$', fontsize=18)
plt.scatter(np.pi/2, 0, color='black', zorder=10)
plt.vlines(x=np.pi/2, ymin=-2, ymax=0, color='black', linestyle='dashed')
fig.show()

The solutions for each of the three regions that we will consider are:

$$\begin{align} & \Psi_I: &&E < V(x) \text{ and WKB works.} \\ & \Psi_{II}: &&E \sim V(x) \text{ and WKB doesn't work.} \\ & \Psi_{III}: &&E > V(x) \text{ and WKB works.} \end{align}$$ The idea for solving this problem is to use a different approximation in region II, and match it with the WKB solution at the boundaries with regions I and III.

4.3.2. Quantum solution for a linear potential - the Airy functions

We now consider a linear approximation for the potential around a turning point, i.e. Taylor expansion of the potential around $x_t$. To keep the notation simple, we will set $x_t=0$ in the following. We can get the general result by shifting the final result back by $x_t$.

In the linear approximation, the potential reads $$V(x) = E + V'(0)x.$$ We then define $$\alpha = \left( \frac{2m}{\hbar^2}V^{'}(0) \right)^{1/3}.$$ Around the turning point, the Schrödinger equation then becomes $$\begin{align} - &\frac{\hbar^2}{2m} \frac{d^2}{dx^2} \psi(x) + (E + V'(0) x)\, \psi(x) = E \psi(x) \\ \Rightarrow \quad & \psi^{''}(x) = \alpha^3 x \psi(x)\\ \Rightarrow \quad &\frac{d^2}{dz^2} \psi(z) = z\psi(z) \quad \text{with $z = \alpha x$}. \end{align}$$ The last equation corresponds to the Airy equation. The solution to this equation is given by Airy functions, and so the general solution to the Airy equation reads: $$\psi(x) = a\, \mathrm{Ai}(\alpha x) + b\, \mathrm{Bi}(\alpha x).$$ The Airy functions are special functions, but luckily enough they have been studied long enough that we know a lot about their properties. In particular, we will need the asymptotic form for $|z|\gg 0$, which for the first Airy function $\mathrm{Ai}(x)$ reads $$\mathrm{Ai}(z) = \left\{ \begin{array}{lc} \frac{1}{\sqrt{\pi}(-z)^{1/4}}\sin\left( \frac{2}{3}(-z)^{3/2} + \frac{\pi}{4} \right), & z <<0,\\ \frac{1}{2\sqrt{\pi}z^{1/4}}e^{-\frac{2}{3}z^{3/2}}, & z >>0. \end{array} \right.$$

The function $\mathrm{Ai}(x)$ is shown together with the asymptotic form in the following plot:

import common
common.configure_plotting()
import numpy as np
from scipy.special import airy
import matplotlib.pyplot as plt

def Ai_asymp_neg(z):
    return 1.0/(np.sqrt(np.pi) * (-z)**0.25) * np.sin(2.0/3.0 * (-z)**1.5 + np.pi/4.)
def Ai_asymp_pos(z):
    return 1.0/(2 * np.sqrt(np.pi) * (z)**0.25) * np.exp(-2.0/3.0 * z**1.5)

xs = np.linspace(-10, 10, 101)
xs_neg = np.linspace(-10, 0, 201, endpoint=False)
xs_pos = np.linspace(0.00001, 10, 101)

fig, ax = plt.subplots(figsize=(10,5))
ax.plot(xs, airy(xs)[0], lw=2.5, label="$\mathrm{Ai}(z)$")
ax.plot(xs_neg, Ai_asymp_neg(xs_neg), "r--", label="asymptote $z<0$")
ax.plot(xs_pos, Ai_asymp_pos(xs_pos), "g--", label="asymptote $z>0$")
ax.legend()
ax.set_xlim(-10,10);
ax.set_ylim(-0.6,2);
ax.set_xlabel("$z$")

fig.show()

$\mathrm{Ai}(x)$ shows oscillating behavior for $z<0$ and decaying behavior for $z>0$, as expected from physical intuition. We can observe that $Ai(x)$ smoothly connects these two regions (unlike the WKB wave function!). Observe that the asymptotic solutions match the expected behavior away from the turning point at $z=0$. In fact, the asymptotic forms agree very well with the Airy function already for small values of $z$, less than one oscillation length!

The asymptotic behavior of the second Airy function $\mathrm{Bi}(x)$ is given by $$\mathrm{Bi}(z) = \left\{ \begin{array}{lc} \frac{1}{\sqrt{\pi}(-z)^{1/4}}\cos\left( \frac{2}{3}(-z)^{3/2} + \frac{\pi}{4} \right), & z <<0,\\ \frac{1}{\sqrt{\pi}z^{1/4}}e^{\frac{2}{3}z^{3/2}}, & z >>0. \end{array} \right.$$

The function $\mathrm{Bi}(x)$ is shown in the following plot. The behaviour is similar as in the previous case, except that this function grows exponentially for $z>0$.

import common
common.configure_plotting()
import numpy as np
from scipy.special import airy
import matplotlib.pyplot as plt


def Bi_asymp_neg(z):
    return 1.0/(np.sqrt(np.pi) * (-z)**0.25) * np.cos(2.0/3.0 * (-z)**1.5 + np.pi/4.)
def Bi_asymp_pos(z):
    return 1.0/(np.sqrt(np.pi) * (z)**0.25) * np.exp(2.0/3.0 * z**1.5)

xs = np.linspace(-10, 10, 101)
xs_neg = np.linspace(-10, 0, 201, endpoint=False)
xs_pos = np.linspace(0.00001, 10, 101)

fig, ax = plt.subplots(figsize=(10,5))

xs = np.linspace(-10, 5, 101)
xs_neg = np.linspace(-10, 0, 1001, endpoint=False)
xs_pos = np.linspace(0.00001, 5, 101)

ax.plot(xs, airy(xs)[2], lw=2.5, label="$\mathrm{Bi}(z)$")
ax.plot(xs_neg, Bi_asymp_neg(xs_neg), "r--", label="asymptote $z<0$")
ax.plot(xs_pos, Bi_asymp_pos(xs_pos), "g--", label="asymptote $z>0$")
ax.legend()
ax.set_ylim(-0.6, 2)
ax.set_xlim(-10, 10)
ax.set_xlabel("$z$")

fig.show()

Airy functions - basis for a solution around the turning point

The Airy functions thus form a basis for a solution around the turning point, smoothly connecting the regions:

$$I: \ E>V(x) \quad \text{ and } \quad III: \ E<V(x).$$

We achieved this by using a specific approximation, namely, by approximating (linearising) the potential $V(x)$ around the turning point, and finding asymptotically exact solutions of the Schrödinger equation.

In the next step, we will use these solutions to "glue" our WKB wave functions in region I and III together.

4.3.3. Derivation of the connection formulas

Once we will have found the solution for region $II$, we will extend the approximation to the overlapping region $I-II$ and $II-III$.

1. Using the linear approximation for the potential in the WKB wave function, we find: $$p(x) \approx \sqrt{2m(E-(E+V'(0)x ))} = \hbar \sqrt{-\frac{2mV'(0)}{\hbar^2} x} = \hbar \alpha^{3/2} \sqrt{-x}$$
2. For region $III, \ x\gg 0$, we have that $E < V(x)$. Therefore, we expect that the solution behaves as exponentially decaying or blowing up. Computing the integral in the exponent of the WKB wave function we find: $$\int_{0}^x |p(x')| dx' = \frac{2}{3} \hbar (\alpha x)^{3/2}$$ By inserting this into the WKB wave function, we arrive at: $$\begin{split} \psi_{III}(x) &= \frac{C}{\sqrt{|p(x)|}}e^{-\frac{1}{\hbar} \int_0^x |p(x')| dx'} + \frac{D}{\sqrt{|p(x)|}}e^{+\frac{1}{\hbar} \int_0^x |p(x')| dx'}\\ &= \frac{C}{(\alpha x)^{1/4}\sqrt{\hbar \alpha}}e^{-\frac{2}{3}(\alpha x)^{3/2}} + \frac{D}{(\alpha x)^{1/4}\sqrt{\hbar \alpha}}e^{+\frac{2}{3}(\alpha x)^{3/2}} \end{split}$$
3. Now, let us compare this result with the result from the asymptotic forms of the Airy functions in region $II$, $$\psi_{II}(x) = \frac{a}{2\sqrt{\pi}z^{1/4}}e^{-\frac{2}{3}z^{2/3}} + \frac{b}{\sqrt{\pi}z^{1/4}}e^{\frac{2}{3}z^{2/3}}$$ By matching the coefficients, i.e. $\psi_{II}(x)=\psi_{III}(x)$, we find that, $$a = 2C \sqrt{\frac{\pi}{\alpha \hbar}},\qquad b = D \sqrt{\frac{\pi}{\alpha \hbar}} \tag{4.7}\label{eq:conn1}$$

4. For $x < 0$, in region $I$, we can apply a similar derivation. However, the behaviour of the solutions will be oscillatory and not decaying: $$\psi_{I}(x) = \frac{\tilde{A}}{\sqrt{p(x)}}e^{\frac{i}{\hbar} \int_x^0 p(x') dx'} + \frac{\tilde{B}}{\sqrt{p(x)}}e^{-\frac{i}{\hbar} \int_x^0 p(x') dx'}$$ It turns out that these formulas become simpler when we perform a change of basis. The wave function $\psi_{I}(x)$ above is written in terms of plane-wave-like solutions. We can redefine the coefficients to write $\psi_{I}(x)$ in terms of $\sin$ and $\cos$ instead: $$\begin{split} \psi_{I}(x) &= \frac{A}{\sqrt{p(x)}}\sin\left(\frac{1}{\hbar} \int_x^0 p(x') dx'+ \theta \right) &&+ \frac{B}{\sqrt{p(x)}}\cos\left(\frac{1}{\hbar} \int_x^0 p(x') dx'+ \theta \right)\\ &= \frac{A}{(-\alpha x)^{1/4}\sqrt{\hbar \alpha}}\sin\left(\frac{2}{3}(-\alpha x)^{3/2} + \theta \right) &&+ \frac{B}{(-\alpha x)^{1/4}\sqrt{\hbar \alpha}}\cos\left(\frac{2}{3}(-\alpha x)^{3/2} + \theta \right) \end{split}$$ Here we included a phase shift $\theta$ for convenience.

5. Let us now compare this solution with the asymptotic solution from region II, $$\psi_{II}(x) = \frac{a}{\sqrt{\pi}(-z)^{1/4}}\sin\left( \frac{2}{3}(-z)^{3/2} + \frac{\pi}{4} \right) + \frac{b}{\sqrt{\pi}(-z)^{1/4}}\cos\left( \frac{2}{3}(-z)^{3/2} + \frac{\pi}{4} \right)$$ By matching the coefficients again, $\psi_{I}(x)=\psi_{II}(x)$, we find $$\theta = \frac{\pi}{4},\quad a = A \sqrt{\frac{\pi}{\hbar \alpha}},\quad b = B \sqrt{\frac{\pi}{\hbar \alpha}} \tag{4.8}\label{eq:conn2}$$

6. As a final step, we now eliminate $a$ and $b$ from the system of equations \eqref{eq:conn1} and \eqref{eq:conn2}, giving $$2C = A,\qquad D =B.$$ We thus found a relation between the coefficients in region $I$ and those in region $III$!

We can now relate wave functions at both sides of the potential, and write down the solution of the system. In fact, these connection formulas dictate that if we know the form of the WKB wave function in e.g. region $I$ (i.e. the coefficients $A$ and $B$), the WKB wave function in region $III$ is uniquely determined. This is very similar to the wave function matching method which you know from quantum mechanics, namely, matching $\psi$ and $\psi'$ at a potential step.

4.3.4. Summary of the connection formulas

We now give the full form of the connection formulas, for a general turning point $x_t$ and potentials of either positive and negative slope.

Positive slope

import common
common.configure_plotting()
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(-1, 1, 100)
y1 = x
E = 0
fig, ax = plt.subplots(figsize=(4, 4))
ax.plot(x, y1, color='black', label=r'V(x)')
ax.hlines(y=E, xmin=-1, xmax=1, color='red', label='E')
ax.vlines(x=0, ymin=-2, ymax=0, linestyle='dashed', color='black')
ax.scatter(0, 0, color='black', zorder=100)
ax.set_ylim(-1, 1);
ax.legend(fontsize=15)
ax.set_xlabel(r'$x_t$', fontsize=15)
ax.set_xticks([]);
ax.set_yticks([]);
fig.show()

$E > V(x)$	$E < V(x)$
$$\frac{2C}{\sqrt{p(x)}} \sin\left( \frac{1}{\hbar} \int_x^{x_t} p(x') dx' + \frac{\pi}{4} \right)$$	$$\frac{C}{\sqrt{\|p(x)\|}} e^{-\frac{1}{\hbar} \int_{x_t}^x \|p(x')\|dx'} \quad \text{ (decaying)}$$
$$\frac{D}{\sqrt{p(x)}} \cos\left( \frac{1}{\hbar} \int_x^{x_t} p(x') dx' + \frac{\pi}{4} \right)$$	$$\frac{D}{\sqrt{\|p(x)\|}} e^{\frac{1}{\hbar} \int_{x_t}^x \|p(x')\|dx'} \quad \text{ (blowing up)}$$

Negative slope

import common
common.configure_plotting()
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(-1, 1, 100)
y1 = -x
E = 0
fig, ax = plt.subplots(figsize=(4, 4))
ax.plot(x, y1, color='black', label=r'V(x)')
ax.hlines(y=E, xmin=-1, xmax=1, color='red', label='E')
ax.vlines(x=0, ymin=-2, ymax=0, linestyle='dashed', color='black')
ax.scatter(0, 0, color='black', zorder=100)
ax.set_ylim(-1, 1);
ax.legend(fontsize=15)
ax.set_xlabel(r'$x_t$', fontsize=15)
ax.set_xticks([]);
ax.set_yticks([]);
fig.show()

$E < V(x)$	$E > V(x)$
$$\text{(decaying) } \quad \frac{C}{\sqrt{\|p(x)\|}} e^{-\frac{1}{\hbar} \int_{x}^{x_t} \|p(x')\|dx'}$$	$$\frac{2C}{\sqrt{p(x)}} \sin\left( \frac{1}{\hbar} \int_{x_t}^{x} p(x') dx' + \frac{\pi}{4} \right)$$
$$\text{(blowing up) } \quad \frac{D}{\sqrt{\|p(x)\|}} e^{\frac{1}{\hbar} \int_{x}^{x_t} \|p(x')\|dx'}$$	$$\frac{D}{\sqrt{p(x)}} \cos\left( \frac{1}{\hbar} \int_{x_t}^{x} p(x') dx' + \frac{\pi}{4} \right)$$

Warning

In these tables, we also indicated which wave functions are decaying and which are blowing up as you move away from the turning point.

When using these formulas, be careful to use exactly the forms given here, and pay attention to the different integral boundaries!

Integration tricks

To deal with WKB problems, one often requires to change the limits of integrations. The simplest case is to flip the integrations limits as, $$\int_{x_0} ^{x_1} f(x) dx = - \int_{x_1} ^{x_0} f(x) dx.$$ Another useful trick when $x_{0}<x<x_1$ is $$\int_{x_0} ^{x} f(x') dx' = \int_{x_0} ^{x_1} f(x') dx' - \int_{x} ^{x_1} f(x') dx'.$$

4.3.5. Validity of the connection formulas

We did a rather involved derivation above taking several approximations, but when are the connection formulas actually valid?

From our derivation, it is clear that there should be a region where the Airy solution $\psi_{II}(x)$ has at least some overlap with the WKB wave functions $\psi_{I}(x)$ and $\psi_{III}(x)$ - otherwise it does not "connect" those solutions well. In the interactive plot below, we show the WKB wave functions as well as the Airy function around the turning point for different values of the energy:

import common
from IPython.display import HTML, display
common.configure_plotting()
import numpy as np
from wkb import plot_patching_region_plt

def V(x):
    return 0.01 * (x+10)**2 -2  
anim = plot_patching_region_plt(V, Erange=np.round(np.linspace(-1.98, 2.4, 101), 2),
                     start=-11, stop=11, m=1, wf_scaling_factor=0.5, y_range=(-5,5))
display(HTML(anim.to_jshtml(default_mode='reflect')))                   

Once Loop Reflect

Validity of the connection formulas

From our derivation, we can deduce the essential conditions:

1. The range of validity of the asymptotic form of the Airy functions must overlap with the region where the linear potential approximation is valid.
2. The range, where the linear potential approximation is valid, must overlap with the range where WKB is a good approximation, i.e. the potential smooth compared to the wave length.

In the interactive plot above, you can see that this model breaks down when the energy is small; then the wavelength becomes longer and the curvature of the potential has a larger influence on the wave functions.

It is good to be aware of the limitations. In practice however, the connection formulas tend to work sufficiently well.

4.3.6. Summary

Take-home message

For the connection formulas, one can always refer to the tables presented above. However, it might turn useful to always remember that:

1. A sine wave function at one side will be connected to a decaying wave function.
2. A cosine wave function at one side will be connected to a blowing up wave function.