final wkb

docs/figures/wkb_currents.py

+import kwant
+import numpy as np
+import matplotlib.pyplot as plt
+"""
+Code used to generate the plot of the probability currents found in wkb_tunel.md
+"""
+
+# Create system
+a = 1
+L = 200
+W = 40
+t = 1
+V = 1
+
+lat = kwant.lattice.square(a=1, norbs=1)
+
+syst = kwant.Builder()
+
+syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
+syst[(lat(x, y) for x in range(80, L-80) for y in range(W))] = 4 * t + V
+syst[lat.neighbors()] = -t
+lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
+lead[(lat(0, j) for j in range(W))] = 4 * t
+lead[lat.neighbors()] = -t
+syst.attach_lead(lead)
+syst.attach_lead(lead.reversed())
+
+# Define scattering states
+E0 = V - 0.7
+E1 = V - 0.15
+scattering_states_0 = kwant.wave_function(syst.finalized(), energy=E0)
+scattering_states_1 = kwant.wave_function(syst.finalized(), energy=E1)
+density = kwant.operator.Current(syst.finalized())
+
+# Plot probability currents and potential profile
+fig, ax = plt.subplots(ncols=1, nrows=3, figsize=(10, 10))
+x = np.arange(0, 200)
+y = np.heaviside(x-80, 0)*np.heaviside(-x+120, 0)
+ax[0].plot(x, y, color='black', label=r'V(x)')
+ax[0].hlines(y=E0, xmin=0, xmax=200, color='darkblue', linestyle='dashed', label='$E_0$')
+ax[0].hlines(y=E1, xmin=0, xmax=200, color='red', linestyle='dashed', label='$E_1$')
+ax[0].legend(fontsize=15)
+ax[0].set_xlabel(r'$x$', fontsize=15)
+ax[0].set_title(r'Potential profile', fontsize=15)
+ax[1].set_title(r'Transmitted wave at $E=E_1$', fontsize=15)
+ax[2].set_title(r'Reflected wave at $E=E_0$', fontsize=15)
+i = 0
+for axx in ax:
+    if i > 0:
+        axx.set_xlabel(r'$x$', fontsize=15)
+        axx.set_ylabel(r'$y$', fontsize=15)
+    axx.set_xticks([]);
+    axx.set_yticks([]);
+    i += 1
+kwant.plotter.current(syst.finalized(), density((scattering_states_1(0)[0])), ax=ax[1])
+kwant.plotter.current(syst.finalized(), density((scattering_states_0(0)[0])), ax=ax[2])
+plt.savefig(fname='wkb_current.png')
+
+"""
+Previous plot
+
+```python inline
+import matplotlib.pyplot as plt
+import numpy as np
+
+E = 6
+amplitude_left = 4.5
+amplitude_right = 0.5
+
+x = np.linspace(0, 5, 200)
+left_side = np.linspace(0, 1, 200)
+middle_side = np.linspace(1, 4, 200)
+right_side = np.linspace(4, 5, 200)
+
+y = 7*np.heaviside(x-1, 0)*np.heaviside(-(x-4), 0)*(2-(x - 2.5)**2/8)
+
+prop_wave = E + amplitude_left*np.cos(4*left_side*np.pi)
+transmitted_wave = E + amplitude_right*np.cos(4*left_side*np.pi)
+
+one_over_l = -(1/3)*np.log((E+amplitude_left)/(E+amplitude_right))
+decaying_wave = (E+amplitude_left)*np.exp((middle_side-1)*one_over_l)
+
+fig, ax = plt.subplots(figsize=(10, 5))
+ax.plot(x, y, color='black', label=r'V(x)')
+ax.plot(left_side, prop_wave, color='darkblue')
+ax.plot(middle_side, decaying_wave, color='darkblue', label=r'$\psi(x)$')
+ax.plot(right_side, transmitted_wave, color='darkblue')
+ax.hlines(y=E, xmin=0, xmax=5, color='red', linestyle='dashed', label='E')
+ax.set_xticks([]);
+ax.set_yticks([]);
+ax.set_xlabel(r'$x$', fontsize=15)
+ax.legend(fontsize=15)
+fig.show()
+```
+
+"""
\ No newline at end of file

docs/wkb_wf.md

@@ -35,26 +35,16 @@ common.configure_plotting()
 
 ## Heuristic derivation
 
-<<<<<<< HEAD
 The WKB wavefunction can be derived following physical intuition. We will follow this approach before doing a formal derivation. Let us recall that:
-=======
-The WKB wavefunction can be derived as follows by considering the following assumptions:
->>>>>>> eed1e4215688d87d6774c8db51da6a2a1bc8b42c
 
 * A smooth potential can be decomposed into piecewise constant pads. 
 * There is no back reflection in a smooth potential.
 
-<<<<<<< HEAD
 A general ansatz for a wavefunction has an amplitude and a phase, that is,
-=======
-The later can be justified by considering sufficiently small pads. In general, a wavefunction has an amplitude and a phase, that is,
-
->>>>>>> eed1e4215688d87d6774c8db51da6a2a1bc8b42c
 $$
 \psi(x) = A(x)e^{i\phi(x)}.
 $$
 
-<<<<<<< HEAD
 ### Propagating waves: $E > V(x)$
 
 
@@ -113,30 +103,29 @@ fig.show()</pre>
 </div>
 
 
-First, let us focus on the phase as can be shown in the right plot above. When the wavefunction moves along a single potential step, the phase acquires a constant shift by moving a certain distance,
+First, let us focus on the phase as can be shown in the right plot above. When the wavefunction moves along a single potential step, the phase acquires a constant shift proportional to the wavevector and the travelled distance,
 $$
 \phi(x_{i+1}) = \phi(x_{i})+ \Delta \phi.
 $$ 
-Over a constant potential, the acquired phase is $\Delta \phi = p(x_i) \Delta x / \hbar$, where $p(x)=\sqrt{2m(E-V(x))}$. The total phase can be obtaining by summing the contributions of all the steps, that is,
+Over a constant potential, the acquired phase is $\Delta \phi = p(x_i) \Delta x / \hbar$, where 
+$$
+p(x)=\sqrt{2m(E-V(x))}.
+$$
+The phase can be obtaining by summing the contributions of all the steps, that is,
 $$
 \psi(x) \sim \exp\left( \frac{i}{\hbar} \sum_{j=0}^N p(x_j) \Delta x\right) =\exp\left( \frac{\pm i}{\hbar} \int_{x_0}^x p(x') d x'\right).
 $$
-In the last equality, we take the limi $\Delta x \rightarrow 0$. Second, let us consider the amplitude. In quantum mechanics, probability currents are conserved. The current is given by, 
-=======
-First, let us focus on the phase. By moving along a single pad, the phase evolve as $\Delta \phi = \phi(x_{i+1})-\phi(x_{i})$. Over a constant potential, the acquired phase is $\Delta \phi = p(x_i) \Delta x / \hbar$, where $p(x)=\sqrt{2m(E-V(x))}$. Therefore, the total phase can be obtaining by summing the contributions of all the pads, that is,
-
->>>>>>> eed1e4215688d87d6774c8db51da6a2a1bc8b42c
+In the last equality, we take the limit $\Delta x \rightarrow 0$. Second, let us consider the amplitude. In quantum mechanics, probability currents are conserved. The current is given by, 
 $$
 j(x) = |\psi(x)^2| v(x) = |A(x)|^2 \frac{p(x)}{m},
 $$
-<<<<<<< HEAD
 where the velocity is $v(x) = p(x)/m$. Since the current is constant, we find that the amplitude of the wavefunction goes as 
 $$
 A(x) \sim 1/\sqrt{p(x)}.
 $$
 Finally, the WKB wavefunction will be given by,
 $$
-\psi_{WKB}(x) = \frac{1}{\sqrt{p(x)}} \exp\left( \frac{i}{\hbar} \int_{x_0}^x p(x') d x'\right).
+\psi_{WKB}(x) = \frac{1}{\sqrt{p(x)}} \exp\left( \pm \frac{i}{\hbar} \int_{x_0}^x p(x') d x'\right).
 $$
 ### Evanescent waves: $E < V(x)$
 We assumed that $E>V(x)$, but this heuristic derivation holds for $E < V(x)$ as well. In this case, the wavefunction does not accumulate a phase, but accumulates a decaying amplitude. On the formal level, 
@@ -150,19 +139,21 @@ $$
 
 !!! summary "Energy dependence and turning points"
 
-    - Observe that for the WKB approximation, the energy enter as a parameter rather than being an outcome of calculations.
+    - Observe that the energy enter as a parameter instead of being an outcome of calculations.
     - When $V(x_0)=E$, the amplitude $A$ of $\psi$ diverges. This means that the WKB approximation is no longer valid. Observe that these points are a function of energy, i.e. $x_0 = x_0(E)$.
 
 ## Formal derivation
-Starting from the Schr$\"{o}$dinger equation,
+Starting from the Schrodinger equation,
 $$
+\begin{equation}\label{eq:schrod}
 -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} \psi(x) + V(x) \psi(x) = E\psi(x) \longrightarrow  \psi^{''}(x) + \frac{p^2(x)}{\hbar^2} \psi(x) = 0.
+\end{equation}
 $$
 Consider $E>V(x)$. Then, $p(x)$ is real. Consider a general wavefunction as ansatz, i.e. $\psi(x) = A(x) \exp(i\varphi(x))$, where $A$ and $\varphi$ are real. Then, the second spatial derivative of $\psi$ is,
 $$ 
 \psi^{''} = A^{''}e^{i\varphi}+2i A^{'} \varphi^{'} e^{i\varphi}- A (\varphi^{'})^2 e^{i\varphi} + i A \varphi^{''} e^{i\varphi}. 
 $$
-By replacing this into the Schr\"odinger equation, one finds,
+By replacing this into the original equation. Then, one finds,
 $$
 A^{''} + 2i A^{'} \varphi^{'} - A ( \varphi^{'} )^2  + i A \varphi^{''}  + \frac{p^2}{\hbar^2} A = 0.
 $$
@@ -183,30 +174,21 @@ Then, we can conclude that the WKB wavefunction is given as,
 $$
 \psi(x) = \frac{1}{\sqrt{p(x)}} \exp\left( \mp \frac{i}{\hbar} \int_{x_0}^x p(x') d x'\right).
 $$
-For the case $E < V(x)$, we use the same replacement as before.
+For the case $E < V(x)$, we use the same argument as before to find exponential decaying and exponential growing solutions.
+
+!!! summary "Smooth potential approximation"
+    
+    We have assumed that $A^{''} << A(\varphi^{'})^2$. This implies that the length scale over wich $V(x)$ changes, $\Delta L$, is much larger than the wavefunction's wavelength, $\lambda$. Note that
+    $$
+    A^{''} \sim \frac{1}{\Delta L} << A(\varphi^{'})^2 \sim \frac{1}{\lambda^2} \longrightarrow \lambda << \Delta L.
+    $$
 
 ## Summary
 
 The WKB wavefunction can be derived from two assumptions: a smooth potential, and a slowly varying wavefunction. The general solution is,
-=======
-
-Note that, for the moment, it is implicitly assumed that $E>V(x)$. In the last equality, the sum is taken to the continuum limit. Second, consider the amplitude. The current is given by $j(x) = |\psi(x)^2| v(x)$ where the velocity is $v(x) = p(x)/m$. Since there is no back reflection, the current is constant. Therefore, $|\psi(x)^2| \sim 1/p(x)$. From here, we find that the amplitude of the wavefunction goes as $A(x) \sim 1/\sqrt{p(x)}$. Then, the WKB wavefunction will be given by,
-
-$$
-\psi_{WKB}(x) = \frac{1}{\sqrt{p(x)}} \exp\left( \frac{i}{\hbar} \int_{x_0}^x p(x') d x'\right).
-$$
-
-### WKB for evanescent waves
-We assumed that $E>V(x)$, but this heuristic derivation holds for $E < V(x)$ as well. In this case, the wavefunction does not accumulate a phase, but accumulates a decaying amplitude. On the formal level, $p(x)=\sqrt{2m(E-V(x))}=i\sqrt{2m(V(x)-E)}=i|p(x)|$. Therefore, the WKB function in a region where $E < V(x)$ is,
->>>>>>> eed1e4215688d87d6774c8db51da6a2a1bc8b42c
 $$
 \begin{split}
 \psi(x)_{E > V(x)} &=  \frac{A}{\sqrt{p(x)}}e^{\frac{i}{\hbar} \int_x^{x_1} p(x') dx'} + \frac{B}{\sqrt{p(x)}}e^{-\frac{i}{\hbar} \int_x^{x_1} p(x') dx'},\\
 \psi(x)_{E < V(x)} &=  \frac{C}{\sqrt{|p(x)}|}e^{\frac{1}{\hbar} \int_x^{x_1} |p(x')| dx'} + \frac{D}{\sqrt{|p(x)|}}e^{-\frac{1}{\hbar} \int_x^{x_1} |p(x')|dx'}.
 \end{split}
 $$
-<<<<<<< HEAD
-As can be seen, its phase sums over all the contributions of $p(x)$, and the amplitude is inversly proportional to it. The WKB approximation breaks down at the turning points $x_0 = x_0(E)$, where $p(x_0)=0$.
-=======
-## Formal derivation
->>>>>>> eed1e4215688d87d6774c8db51da6a2a1bc8b42c

src/adiabatic_theorem.md

@@ -85,7 +85,9 @@ On the other hand, for a smooth change in the potential, the adiabatic theorem
 tells us to expect that the ground state remains in the ground
 state even as the wall moves, so long as it moves slowly enough.
 
+
 ```python inline
+
 from adiabatic import make_adiabatic_potential_animation
 
 L = 200
@@ -96,4 +98,7 @@ def pot_slow(i, t):
 
 anim = make_adiabatic_potential_animation(potential=pot_slow)
 display(HTML(common.to_svg_jshtml(anim)))
-```
\ No newline at end of file
+
+```
+
+## Summary

src/figures/wkb_currents.py

+import kwant
+import numpy as np
+import matplotlib.pyplot as plt
+"""
+Code used to generate the plot of the probability currents found in wkb_tunel.md
+"""
+
+# Create system
+a = 1
+L = 200
+W = 40
+t = 1
+V = 1
+
+lat = kwant.lattice.square(a=1, norbs=1)
+
+syst = kwant.Builder()
+
+syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
+syst[(lat(x, y) for x in range(80, L-80) for y in range(W))] = 4 * t + V
+syst[lat.neighbors()] = -t
+lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
+lead[(lat(0, j) for j in range(W))] = 4 * t
+lead[lat.neighbors()] = -t
+syst.attach_lead(lead)
+syst.attach_lead(lead.reversed())
+
+# Define scattering states
+E0 = V - 0.7
+E1 = V - 0.15
+scattering_states_0 = kwant.wave_function(syst.finalized(), energy=E0)
+scattering_states_1 = kwant.wave_function(syst.finalized(), energy=E1)
+density = kwant.operator.Current(syst.finalized())
+
+# Plot probability currents and potential profile
+fig, ax = plt.subplots(ncols=1, nrows=3, figsize=(10, 10))
+x = np.arange(0, 200)
+y = np.heaviside(x-80, 0)*np.heaviside(-x+120, 0)
+ax[0].plot(x, y, color='black', label=r'V(x)')
+ax[0].hlines(y=E0, xmin=0, xmax=200, color='darkblue', linestyle='dashed', label='$E_0$')
+ax[0].hlines(y=E1, xmin=0, xmax=200, color='red', linestyle='dashed', label='$E_1$')
+ax[0].legend(fontsize=15)
+ax[0].set_xlabel(r'$x$', fontsize=15)
+ax[0].set_title(r'Potential profile', fontsize=15)
+ax[1].set_title(r'Transmitted wave at $E=E_1$', fontsize=15)
+ax[2].set_title(r'Reflected wave at $E=E_0$', fontsize=15)
+i = 0
+for axx in ax:
+    if i > 0:
+        axx.set_xlabel(r'$x$', fontsize=15)
+        axx.set_ylabel(r'$y$', fontsize=15)
+    axx.set_xticks([]);
+    axx.set_yticks([]);
+    i += 1
+kwant.plotter.current(syst.finalized(), density((scattering_states_1(0)[0])), ax=ax[1])
+kwant.plotter.current(syst.finalized(), density((scattering_states_0(0)[0])), ax=ax[2])
+plt.savefig(fname='wkb_current.png')
+
+"""
+Previous plot
+
+```python inline
+import matplotlib.pyplot as plt
+import numpy as np
+
+E = 6
+amplitude_left = 4.5
+amplitude_right = 0.5
+
+x = np.linspace(0, 5, 200)
+left_side = np.linspace(0, 1, 200)
+middle_side = np.linspace(1, 4, 200)
+right_side = np.linspace(4, 5, 200)
+
+y = 7*np.heaviside(x-1, 0)*np.heaviside(-(x-4), 0)*(2-(x - 2.5)**2/8)
+
+prop_wave = E + amplitude_left*np.cos(4*left_side*np.pi)
+transmitted_wave = E + amplitude_right*np.cos(4*left_side*np.pi)
+
+one_over_l = -(1/3)*np.log((E+amplitude_left)/(E+amplitude_right))
+decaying_wave = (E+amplitude_left)*np.exp((middle_side-1)*one_over_l)
+
+fig, ax = plt.subplots(figsize=(10, 5))
+ax.plot(x, y, color='black', label=r'V(x)')
+ax.plot(left_side, prop_wave, color='darkblue')
+ax.plot(middle_side, decaying_wave, color='darkblue', label=r'$\psi(x)$')
+ax.plot(right_side, transmitted_wave, color='darkblue')
+ax.hlines(y=E, xmin=0, xmax=5, color='red', linestyle='dashed', label='E')
+ax.set_xticks([]);
+ax.set_yticks([]);
+ax.set_xlabel(r'$x$', fontsize=15)
+ax.legend(fontsize=15)
+fig.show()
+```
+
+"""
\ No newline at end of file

src/wkb_connection.md

@@ -29,11 +29,11 @@ common.configure_plotting()
 
 ## Turning points
 
-If both regions $E>V(x)$ and $E < V(x)$ are present, the a turning point $E = V(x_0)$ appears in the middle as can be observed in the plot above. At this point we have that $p(x_0) \rightarrow 0$, and then,
+If both regions $E>V(x)$ and $E < V(x)$ are present, the a turning point $E = V(x_0)$ appears in the middle as can be observed in the plot above. At this point we have that,
 $$
 \psi_{WKB}(x_0) \sim \frac{1}{\sqrt{p(x_0)}} \rightarrow \infty.
 $$
-That is, the wavelength becomes infinite and the WKB is no longer valid. However, we can solve this problem by defining three regions, solving for the wavefunctions in each region, and finally matching the solutions. The regions are indicated in the plot below
+That is, the wavelength becomes infinite and the WKB is no longer valid. However, we can solve this problem by defining three well-behaved regions, solving for the wavefunctions in each region, and finally matching the solutions. The regions are indicated in the plot below
 
 ```python inline
 import matplotlib.pyplot as plt
@@ -62,10 +62,11 @@ plt.vlines(x=np.pi/2, ymin=-2, ymax=0, color='black', linestyle='dashed')
 fig.show()
 ```
 
-* Region I: $E < V(x)$ and WKB works.
-* Region II: $E \sim V(x)$ and WKB doesn't work.
-* Region III: $E > V(x)$ and WKB works.
+The solutions for each of the three regions that we'll consider are:
 
+* $\Psi_I$: $E < V(x)$ and WKB works.
+* $\Psi_{II}$: $E \sim V(x)$ and WKB doesn't work.
+* $\Psi_{III}$: $E > V(x)$ and WKB works.
 
 The idea to solve this problem is by using a different approximation in region II, and matching it with the WKB solution at the boundaries with regions I and III.
 
@@ -221,30 +222,29 @@ $$
 $$
 We have found a relation between the coefficients at region I with those at region III. Therefore, we can relate wavefunctions at both sides of the potential, and write down the solution of the system. Let us now observe what terms are connected at each side of the potential for the following two cases plotted below.
 
+
+For a potential with a positive slope,
+
 ```python inline
 import matplotlib.pyplot as plt
 import numpy as np
 
 x = np.linspace(-1, 1, 100)
 y1 = x
-y2 = -x
 E = 0
-fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
-axes[0].plot(x, y1, color='black', label=r'V(x)')
-axes[1].plot(x, y2, color='black', label=r'V(x)')
-for ax in axes:
-    ax.hlines(y=E, xmin=-1, xmax=1, color='red', label='E')
-    ax.set_xticks([]);
-    ax.set_yticks([]);
-    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);
-axes[0].legend(fontsize=15)
-axes[0].set_xlabel(r'$x_t$', fontsize=15)
-axes[1].set_xlabel(r'$x_t$', fontsize=15)
+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()
 ```
-For the left plot,
+
 $$
 \begin{split}
 \text{term at }E >  V(x) &\leftrightarrow \text{term at }E < V(x)\\
@@ -253,7 +253,28 @@ $$
 \end{split}
 $$
 
-For the right plot, 
+For a potential with a negative slope,
+
+```python inline
+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()
+```
+
 $$
 \begin{split}
 \text{term at }E <  V(x) &\leftrightarrow \text{term at }E > V(x)\\
@@ -296,6 +317,9 @@ plot_patching_region(V, Erange=np.round(np.linspace(-1.9, 2, 51), 2),
 
 ## Summary
 
-From the connection formulas one can always look back at the table presented above. However, one can remember as well that a sine wavefunction at one side will always be connected to a decaying wavefuncion at the other side, and a cosine wavefunction with a blowing up wavefunction.
+From the connection formulas one can always look back at the table presented above. However, one can always remember that:
+
+* A sine wavefunction at one side will be connected to a decaying wavefuncion.
+* A cosine wavefunction at one side will be connected to a blowing up wavefunction.
 
 

src/wkb_smooth.md

@@ -2,18 +2,29 @@
 title: The WKB approximation
 ---
 
+# A particle moving in a potential
+
+!!! success "Expected prerequisites"
+
+    Before the start of this lecture, you should be able to:
+
+    - Describe the motion of a classical particle in a potenial.
+
 !!! summary "Learning goals"
 
     After this lecture you will be able to:
 
-    - Conceptually describe the motion of a particle in a smooth or abrupt potential for the classical and quantum cases.
+    - Describe the motion of a quantum particle in a smooth or abrupt potential.
+
+## Classical intuition
 
-# Quantum mechanics in a smooth potential
+The WKB approximation is used to describe a particle of mass $m$ in a one-dimensional potential. To call some intuition, let us recall the classical case: consider a particle of energy $E$ moving in a potential $V(x)$. If $E > V$, then the particle will always move above the potential. There will always be some non-zero kinetic energy irrespective of the shape of the potential. If $E < V$, will reach a maximum position $x_0$ such that $E = V(x_0)$. Then, the particle will be reflected. 
 
-The WKB approximation is used to describe a particle of mass $m$ in a one-dimensional potential. To call some intuition, let us recall the classical case. Consider a particle of energy $E$ moving in a potential $V(x)$. If $E > V$, then the particle will always move about the potential. It will have always some non-zero kinetic energy irrespective of the shape of the potential. If $E < V$, the particle will be reflected. 
+## Quantum case
 
+### Smooth potential
 
-Let us consider the quantum case where we talk about wave packets. In the following animation we describe the cas where $V(x)$ is a smooth potential describing a dip. One can observe that the intuition from classical mechanics holds because the particle moves through the potential with a sightly increase in velocity along the dip.
+Let us consider the quantum case where we talk about wave packets. In the following animation we describe the case where $V(x)$ is a smooth potential describing a dip. One can observe that the intuition from classical mechanics holds because the particle moves through the potential with a sightly increase in velocity along the dip.
 
 ```python inline
 import common
@@ -36,12 +47,13 @@ anim = make_wave_packet_animation(L=3500, pot_func=pot_gauss,
 display(HTML(common.to_svg_jshtml(anim)))
 ```
 
+### Abrupt potential
 
-One the other hand, if $V(x)$ describes an abrupt potential, one can observe in the following animation that there is some backreflection everytime that the wavepacket crosses an abrupt step. This is an effect that is not present in the classical case.
-
+One the other hand, if $V(x)$ describes an abrupt potential, one can observe in the following animation that the wavefunction reflects back whenever it crosses an abrupt step.
 
 
 ```python inline
+
 from wkb import make_wave_packet_animation
 import math
 
@@ -56,4 +68,13 @@ anim = make_wave_packet_animation(L=3500, pot_func=pot_step,
                                   energy=0.01)
 
 display(HTML(common.to_svg_jshtml(anim)))
-```
\ No newline at end of file
+```
+
+This is an effect that is not present in the classical case.
+
+## Summary 
+
+For a quantum particle with $E > V(x)$:
+
+* When a particle moves along a smooth potential, it resembles the classical case.
+* When a particle moves along an abrupt potential, there's always some backreflection.
\ No newline at end of file

src/wkb_tunel.md

@@ -31,41 +31,11 @@ common.configure_plotting()
 
 ## Physical intuition
 
-```python inline
-import matplotlib.pyplot as plt
-import numpy as np
-
-E = 6
-amplitude_left = 4.5
-amplitude_right = 0.5
-
-x = np.linspace(0, 5, 200)
-left_side = np.linspace(0, 1, 200)
-middle_side = np.linspace(1, 4, 200)
-right_side = np.linspace(4, 5, 200)
-
-y = 7*np.heaviside(x-1, 0)*np.heaviside(-(x-4), 0)*(2-(x - 2.5)**2/8)
-
-prop_wave = E + amplitude_left*np.cos(4*left_side*np.pi)
-transmitted_wave = E + amplitude_right*np.cos(4*left_side*np.pi)
+![Probability current](figures/wkb_current.svg)
 
-one_over_l = -(1/3)*np.log((E+amplitude_left)/(E+amplitude_right))
-decaying_wave = (E+amplitude_left)*np.exp((middle_side-1)*one_over_l)
-
-fig, ax = plt.subplots(figsize=(10, 5))
-ax.plot(x, y, color='black', label=r'V(x)')
-ax.plot(left_side, prop_wave, color='darkblue')
-ax.plot(middle_side, decaying_wave, color='darkblue', label=r'$\psi(x)$')
-ax.plot(right_side, transmitted_wave, color='darkblue')
-ax.hlines(y=E, xmin=0, xmax=5, color='red', linestyle='dashed', label='E')
-ax.set_xticks([]);
-ax.set_yticks([]);
-ax.set_xlabel(r'$x$', fontsize=15)
-ax.legend(fontsize=15)
-fig.show()
-```
+A tunneling problem describes the situation found in the picture above. That is, a particle approaching a barrier from the left side, and being transmitted to the right side. On the top plot we observe a potential profile. There are two energies indicated by dashed lines. On the lower two plots we show the probability currents associated to each energy. Note that the system is two-dimensional. For $E=E_1$, on the left side of the barrier there are reflected and incident components, while on the right side there is a transmitted wave that has been attenuated due to the potential barrier. On the lower plot, the one with $E=E_0$, we observe that a particle with lower energy does not propagate to the other side, but is completely reflected. 
 
-A tunneling problem describes a situation where outside of a tunneling region, $0< x < L$, the particle propagates freely, i.e. $E>V(x)$. In such regions, we know that the wave can be incident, transmitted or reflected waves. That is,
+Now, let us describe this problem from a mathematical perspective. Outside of a tunneling region, $0< x < L$, the particle propagates freely, i.e. $E>V(x)$. The solution in this region can be an incident, transmitted or reflected wave. That is,
 $$
 \psi(x) = \left\{ \begin{array}{cc}
 A e^{ikx} + B e^{-ikx}, & x <0,\\
@@ -77,7 +47,7 @@ From quantum mechanics, recall that the tunneling probability is given as,
 $$
 t = \left| \frac{T}{A} \right|^2 \frac{v_R}{v_L}.
 $$
-The factor involving the velocities comes the probability current that is tunneling at the other side of the potential. It does not appear if the energy at both sides is the same. To solve for the region $0 < x < L$, we will use the WKB approximation. Observe that if the potential would be flat, then there would be a exponential decaying factor that depends on the distance. However, to account for an arbitrary potential shape, the WKB wavefunction is,
+The factor involving the velocities comes the conservation of probability current. It does not appear if the energy at both sides is the same. To solve for the region $0 < x < L$, we will use the WKB approximation. For a constant potential, there would be a exponential decaying factor proportional to the distance. However, to account for an arbitrary potential shape, the WKB wavefunction is,
 $$
 \psi_{WKB}(x) = \frac{1}{\sqrt{|p(x)|}} \left( C e^{ -\frac{1}{\hbar} \int_{0}^x |p(x')| d x'} + D e^{ \frac{1}{\hbar} \int_{0}^x |p(x')| d x'} \right).
 $$
@@ -90,6 +60,8 @@ $$
 t = \left| \frac{T}{A} \right|^2 \frac{v_R}{v_L} \sim e^{-2\gamma},\qquad \gamma = \frac{1}{\hbar}\int_0^L |p(x')|dx'.
 $$
 
+This result can be used to describe the plots on top of this section. Observe that $|E_0 - V(x)| >> |E_1 - V(x)| $. This means that the wave associated with $E_0$ will have a exponential factor that is much larger than the one with $E_1$. Therefore, the one with lower energy will have an exponentially smaller transmitted component than the one with higher energy.
+
 ## WKB for tunneling
 
 The tunneling amplitude can be found using the WKB approximation, and the result corresponds to the one derived above. In particular, we'll use the transfer matrix method. Let us recall that the WKB wavefunction for propagating waves is,
@@ -104,25 +76,20 @@ import numpy as np
 
 x = np.linspace(-1, 1, 100)
 y1 = x
-y2 = -x
 E = 0
-fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
-axes[0].plot(x, y1, color='black', label=r'V(x)')
-axes[1].plot(x, y2, color='black', label=r'V(x)')
-for ax in axes:
-    ax.hlines(y=E, xmin=-1, xmax=1, color='red', label='E')
-    ax.set_xticks([]);
-    ax.set_yticks([]);
-    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);