## Introduction

The WBK method is an approximation technique for approximating solutions to the time-independent Schrödinger equation
with non-constant potential.
It is named after Wentzel, Kramers, and Brillouin. In 1923, three years prior to the development of the WKB method, a
British mathematician Harold Jeffreys also developed it. But it was not very popular back then. In England, people
also call it the JWKB method for the contribution made by Jeffreys.


## Derivation of WKB

The Schrödinger equation with non-constant potential is,

\begin{equation}
    -\frac{\hbar^2}{2m} \frac{d^2 \psi}{dx^2} + V(x)\psi = E\psi.
\end{equation}

We assume $E>V(x)$ and rewrite it as

\begin{equation}
    \frac{d^2 \psi}{dx^2}= -\frac{p^2}{\hbar^2} \psi,
\end{equation}

where $p(x)=\sqrt{2m[E-V(x)]}$ is the classical momentum of a particle with total energy E and potential $V(x)$.

The general solution of this equation has the form

\begin{equation}
    \psi(x) = A(x)e^{i\phi(x)},
\end{equation}

where $A(x)$ is the amplitude and $\phi(x)$ is the phase, both as a function of x. Then we can calculate the first and
the second derivative of $\psi$ and plug them back to the Schrödinger equation to get

\begin{equation}
    A'' + 2i A' \phi' + iA\phi'' - A(\phi')^2 = - \frac{p^2}{\hbar^2}A.
\end{equation}

The real and imaginary part each gives us an equation.

The imaginary part is

\begin{equation}
    2A'\phi' + A\phi'' = 0.
\end{equation}

Its solution is

\begin{equation}
    A = \frac{C}{\sqrt{\phi'}},
\end{equation}

where $C$ is a real number.

The real part is

\begin{equation}
    A'' - A (\phi')^2 = -\frac{p^2}{\hbar^2}A.
\end{equation}

Unfortunately this equation cannot be solved in most cases. Here, the approximation comes in. By assuming $A(x)$ varies
slowly as a function of $x$, we can make $A''=0$. Then we can solve it to get

\begin{equation}
    \phi(x) = \pm \frac{1}{\hbar} \int p(x)dx.
\end{equation}

Therefore, the general solution is given by

\begin{equation}
    \psi(x) = \frac{C}{\sqrt{p(x)}} e^{\pm \frac{i}{\hbar} \int p(x)dx}.
\end{equation}

For $E<V$, we need to take $p(x)$ to absolute value to keep it a real function.

Please see Chapter 8 of Griffiths' Introduction to Quantum Mechanics for a more detailed discussion. Professor
Barton Zwiebach's lecture presented another derivation. It's very interesting as well.


## Validity of WKB

The argument we used in the derivation for the approximation is rather mathematical. A more physical way to look it is
that the probability of finding a particle at location $x$ is inversely proportional to its momentum. It agrees with our
intuition. If a particle moves fast at one location, it's less likely to be found by us. Professor Zwiebach derived some
more formal inequalities to justify the validity of the WKB approximation. Dropping the second derivative term leads
to the inequality.

\begin{equation}
    \lambda \frac{d \lambda(x)}{dx} << \lambda.
\end{equation}

There are several others, but this one makes the most sense physically. It says that the wavelength change over one
de Broglie wavelength is much smaller than the de Broglie wavelength itself. This is the physical assumption of this
approximation.

For a potential barrier with sloping walls, if the energy of the incoming particle is less than the potential of the
barrier in one portion of it. The points where $E=V$ is called the turning points. The WKB method is only valid away from
those turning points. As we move close to the turning points, the derivative of the local de Broglie wavelength blows up.

## The Connection Formulas

To approximate the solution near the turning points, people developed the connection formulas. Its derivation is rather
involved. So I decide to not talk about it here.

The idea is that we use a straight line to approximate the potential near the turning points. Then find a "patching"
wave function to connect both sides. The connection formulas that joining the WKB solutions at a turning point $x_1$ is

\begin{equation}
    \psi(x) =
    \begin{cases}
    & \frac{2D}{\sqrt{p(x)}} sin[ \frac{1}{\hbar} \int_{x}^{x_1}p(x')dx' + \frac{\pi}{4} ], \quad \text{if } x<x_1; \\
    & \frac{D}{\sqrt{|p(x)|}} exp[ -\frac{1}{\hbar} \int_{x}^{x_1}p(x')dx' + \frac{\pi}{4} ], \quad \text{if } x>x_1. \\
    \end{cases}
\end{equation}

Using WKB with connection formulas, we can calculate the bound state energies and transmission coefficient without actually
solving the Schrödinger equation.


## Analytical Solution of a Finite Rectangular Potential Barrier

We derived the expressions for the transmission coefficient the conditions that the wave function and its first
derivative must be continues everywhere. We wrote down the boundary conditions

\begin{aligned}
\psi_1 (0) &= \psi_2 (0) \\
\psi_1' (0) &= \psi_2' (0) \\
\psi_2 (L) &= \psi_3 (L) \\
\psi_2' (L) &= \psi_3' (L).
\end{aligned}

Then we solved for the transmission coefficient $T$ for three cases: $E<V_0$, $E>V_0$, and $E=V_0$. We have done one
of those three cases in homework 4. The other two are similar.

### Case 1: $E<V_0$
We have
\begin{equation}
T = \bigg[ 1 + \frac{V_0^2 sinh^2(k L)}{4E(V_0-E)} \bigg]^{-1},
\end{equation}

where $k=\sqrt{2m(V_0 - E)} / \hbar$. This is the quantum tunneling case. Classically, particles with energy less than
the potential barrier have zero chance going through. But in a quantum system, the particle has some probability
tunneling through. This probability is exponentially decayed as the width of of the barrier increases.

### Case 2: $E>V_0$
\begin{equation}
T = \bigg[ 1 + \frac{V_0^2 sin^2(k L)}{4E(E-V_0)} \bigg]^{-1},
\end{equation}

where $k=\sqrt{2m(E - V_0)} / \hbar$. Classically, particles with energy greater than
the potential barrier would go through certainly. But the quantum scenario, the transmission coefficient actually
oscillates as a function of the width and potential of the barrier and the energy of the particle.

### Case 3: $E=V_0$
\begin{equation}
T = \bigg[ 1 + mL^2 V_0 / 2\hbar^2 \bigg]^{-1}.
\end{equation}

In the case $E=V_0$ the transmission coefficient depends on the width and potential of the barrier.

### Plot of the Transmission Coefficient T

In the following graph, I plotted the three cases separately and combined them to form a single graph. We can see the
exponential decay of the transmission coefficient as the width of the barrier increases. And in the region that $E>V_0$,
$T$ oscillates and trends to 1.


```python
import numpy as np
import scipy.constants as sc
import matplotlib.pyplot as plt

# suppress divide by zero warnings
import warnings; warnings.simplefilter('ignore')
```


```python
# constants are imported from the scipy physical constants library
hbar = sc.physical_constants['reduced Planck constant in eV s'][0]
m_e = sc.physical_constants['electron mass energy equivalent'][0]
```


```python
T_l = lambda e_, l_, v_: (1 + (v_**2 * np.sinh(l_*np.sqrt(2*m_e*(v_-e_))/hbar)**2)/(4*e_*(v_-e_)))**-1
T_eq = lambda e_, l_, v_: (1 + m_e*l_**2 * v_ / (2*hbar**2))**-1
T_g = lambda e_, l_, v_: (1 + (v_**2 * np.sin(l_*np.sqrt(2*m_e*(e_-v_))/hbar)**2)/(4*e_*(e_-v_)))**-1
```


```python
n = 100
l = 10E-9 # m
vs = [0.05, 0.1, 0.5, 1] # eV

plt.figure()
for v in vs:
    es_lv = np.linspace(0, v*((n-1)/n), n-1)
    es_gv = np.linspace(v*((n+1)/n), v*3, n*2)
    es = np.append(es_lv, v)
    es = np.concatenate((es, es_gv))
    leq = np.append(T_l(es_lv, l, v), T_eq(v, l, v))
    Ts = np.concatenate((leq, T_g(es_gv, l, v)))
    es_v = es / v
    plt.plot(es_v, Ts, label=f'v={v} eV')

plt.vlines(1, 0, 1, linestyles='dashed')
plt.xlabel('E / V_0')
plt.ylabel('T')
plt.legend()
plt.show()

```

## The WKB approximation of a Finite Rectangular Potential Barrier

### WKB inside of the Barrier
Since the wave function of an electron passing through a finite rectangular potential barrier is known. We can compare
the exact analytical solution with the one getting from the WKB approximation. The first thing I tried is to solve for the
transmission coefficient using the WKB wave function inside the barrier.

\begin{equation}
    \psi_2(x) = \frac{C}{\sqrt{p}} e^{\frac{p}{\hbar}x} + \frac{D}{\sqrt{p}} e^{-\frac{p}{\hbar}x}, \\
\end{equation}

where $p=\sqrt{2m(V_0 - E)}$. The phase part is the result of the integral $\int_{0}^{L}|p(x')|dx'$. Then we use the
boundary conditions to solve for the transmission coefficient T.

\begin{align}
    T = \bigg\vert \frac{F}{A} \bigg\vert^2
      = \bigg\vert \bigg[
      \frac{1}{2} e^{ikL-\frac{pL}{\hbar}} \big(1+\frac{p}{ik\hbar} \big) \big(1+\frac{ik\hbar}{p} \big)
         + \frac{1}{2} e^{ikL+\frac{pL}{\hbar}} \big(1-\frac{p}{ik\hbar} \big) \big(1-\frac{ik\hbar}{p} \big)
         \bigg]^{-1} \bigg\vert^2.
\end{align}

This equation could be farther simplified. But since we are plotting it with a program. This form works perfectly fine.

For a rectangular potential barrier with a flat top. The second derivative term of the amplitude of the wave equation is
zero. So dropping it in the approximation would not affect the result. Any of its deviation from the exact analytical
solution should be because WKB doesn't work very well near the boundaries.

### WKB for Potential with Sloping Walls

The way of getting the exponential decay part of the transmission coefficient was presented in both homework 4 and
the Griffiths book. It's give by

\begin{equation}
    T \approx e^{-2\gamma}, \\
\end{equation}
where $\gamma = \frac{1}{\hbar}\int_0^L |p(x)|dx$. The problem 8.10 of the Griffiths book also talks about the WKB solution of a scattering from
barrier with sloping walls. It uses the connection formula. The expression is

\begin{equation}
    T = \frac{e^{-2\gamma}}{ [1-(e^{-\gamma /2})^2]^2 }, \\
\end{equation}
where $\gamma = \frac{1}{\hbar}\int_0^L |p(x)|dx$.
I'm not sure if it applies to a finite rectangular potential barrier with vertical walls. But it would be interesting to
see how it looks.



```python
gamma = lambda e_, l_, v_: np.exp(l_/hbar*np.sqrt(2*m_e*(v_-e_)))
T_wkb = lambda e_, l_, v_: np.exp(-2*gamma(e_, l_, v_)) / np.power(1 - np.power(np.exp(-gamma(e_, l_, v_)/2), 2), 2)

p = lambda e_, v_: np.sqrt(2*m_e*(v_ - e_))
k = lambda e_: np.sqrt(2*m_e*e_)/hbar
T_test = lambda e_, v_, l_: np.power(np.absolute( ((1+p(e_, v_)/(k(e_)*hbar*1j))*(0.5*np.exp(k(e_)*l_*1j - p(e_, v_)*l_/hbar)*(1+(k(e_)*hbar*1j)/p(e_, v_))) + (1-p(e_, v_)/(k(e_)*hbar*1j))*(0.5*np.exp(k(e_)*l_*1j + p(e_, v_)*l_/hbar)*(1-(k(e_)*hbar*1j)/p(e_, v_))))**-1 ) , 2)
```


```python
l = 10E-9
v = 1
es_lv = np.linspace(0, v*((n-1)/n), n-1)
es_vv = es_lv / v
plt.plot(es_vv, T_test(es_lv, v, l), label='WKB')
plt.plot(es_vv, T_wkb(es_lv, l, v), label='WKB-connection')
plt.plot(es_vv, T_l(es_lv, l, v), label='Analytical')
plt.xlabel('E / V_0')
plt.ylabel('T')
plt.legend()
plt.show()
```

The results agree with the approximation. As the energy of the particle approaches the potential of the barrier, the
transmission coefficient starts to deviate from the exact solution.


```python
l = 5E-9
v = 1
es_lv = np.linspace(0, v*((n-1)/n), n-1)
es_vv = es_lv / v
plt.plot(es_vv, T_test(es_lv, v, l), label='WKB')
plt.plot(es_vv, T_wkb(es_lv, l, v), label='WKB-connection')
plt.plot(es_vv, T_l(es_lv, l, v), label='Analytical')
plt.xlabel('E / V_0')
plt.ylabel('T')
plt.legend()
plt.show()
```

Another assumption of the WKB approximation is that the size of the aperture is much larger than the de Broglie wavelength of the incoming particle. So I reduced the width of the potential barrier by half and tested how it affects the
transmission coefficient T. For the exact solution, T increases overall. It agrees with our intuition since a particle
is more likely to tunnel through a narrower potential barrier. But in this case, the WKB solution deviates even more
from the exact solution.

## References

- Griffiths, D. J. (2004). Introduction to Quantum Mechanics (2nd Edition). Pearson Prentice Hall. ISBN: 0131118927
- Barton Zwiebach. 8.06 Quantum Physics III. Spring 2018. Massachusetts Institute of Technology: MIT OpenCourseWare,
https://ocw.mit.edu. License: Creative Commons BY-NC-SA.


64c61d24dfd6fafd62df340dcbe26821396bb0d4
ipynb
Jupyter Notebook
wkb/wkb.ipynb
VirosaLi/quantum-electronics
90d3188583d9b6172c4597f74640ae534b4b4579
MIT
wkb/wkb.ipynb
VirosaLi/quantum-electronics
90d3188583d9b6172c4597f74640ae534b4b4579
MIT
wkb/wkb.ipynb
VirosaLi/quantum-electronics
90d3188583d9b6172c4597f74640ae534b4b4579
MIT
2020-10-03T13:35:15.000Z
2020-10-03T13:35:15.000Z
Qwen/Qwen-72B
1. YES
2. YES