Kronig-Penney again

We once again do what we did last time, which is put 10 repetitions of the Kronig-Penney model in a box. As basis functions we use quite reasonably, the exact solutions of the particle-in-a-box problem:

\begin{align} \psi _n (x)&= \sqrt{\frac{2}{L}}\sin (\frac{n \pi x}{L} )\\ E_n &=\frac{n^2 \pi ^2}{2L ^2} \end{align}

Now, for these functions the kinetic energy part ( \(-1/2 d^2/dx^2\) ) of our Hamiltonian just gives us the energy levels of the particle in a box, \(E_n \) (can you see why?). As for the repeating constant potential, the integrals can be performed (and if you want to practice integration, you might want to do it by hand), but I will just give the result:

\begin{align} H _{ji} &= E_j\delta _{ji}+V_0\sum _{r=1}^{10}h_{ji}(x_r,b) \\ h_{ji}(s,b) &= \frac{2}{L}\int _{barriers}\sin \bigg( \frac{j\pi x}{L} \bigg)\sin\bigg( \frac{i\pi x}{L} \bigg) \end{align}

The last integral is just done over the barriers where the potential is not zero. This gives, if barrier width is b and the locations of barrier centers are \(x_r \), the result:

\begin{align} h_{ji}(s,b) &= F_{ji}(s+b/2)-F_{ji}(s-b/2)\label{kronigbasis1}\\ F_{ii}(x) &= \frac{x}{L}-\frac{\sin (2\pi i x/L)}{2\pi n}\\ F_{ji}(x) &= \frac{\sin [(j-i)\pi x/L]}{\pi (j-i)} - \frac{\sin [(i+j)\pi x/L]}{\pi (j+i)}, i\neq j \label{kronigbasis2} \end{align}

We are now ready to implement a code to calculate these matrix elements. This time, rather than setting certain areas of a matrix to a constant value, we implement precisely the functions \eqref{kronigbasis1}-\eqref{kronigbasis2}, and the result is:

import matplotlib.pyplot as plt
from scipy.linalg import eigh
import numpy as np
#Number of barriers
nbar = 10
a = 1
b = 1/6
V0 = 100
xi = np.zeros(nbar)
npoints = 100
H = np.zeros((npoints,npoints))
def fnm(x,n,m):
	if n==m:
		return x/L - \
		return np.sin((m-n)*np.pi*x/L)/\
		(np.pi*(m-n))- \
def hnm(x,n,m):
	return fnm(x+b/2,n,m) - fnm(x-b/2,n,m)        
for i in range(nbar):
	xi[i] = (-0.5)*a + (i+1) * a

for i in range(npoints):
	for j in range(i,npoints):
		for k in range(nbar):
			H[i,j] += V0\

		if i == j:
			H[i,j] += np.pi**2\
		H[j,i] = H[i,j]

e, vec = eigh(H)
x = np.zeros(30)
for i in range(30):
	x[i] = (i+1)*np.pi/L
plt.xlabel("Wave number k of solution")

The solution may look somewhat different than in the discrete case (due to numerics), but you should certainly see the band gaps.

Numerical integration (quadrature)

In the previous code, we calculated the Hamiltonian integrals by using the exact solution. However, you could also use numerical integration, which I'll now explain for completeness.

Let's remind ourselves of the definition of the (Riemannian) integral. Essentially, you can define it using functions of the type

\begin{align} \int _{a}^b f(x)\mathrm{d}x \approx \sum _{i=1}^N f(x_i)\Delta x \end{align}

when the interval \([a,b]\) is divided in to \(N\) partitions, as we have done in our examples here. Hence, for our purpose, a decent approximation to the integral numerically is just to take this sum for a sufficiently large \(N\). In code:

import numpy as np

theta = np.linspace(0,4,1000)
dtheta = x[1]-x[0]
f = np.sin(theta)

The analytical solution gives the integral of \(\sin(x)\) between 0 and 4 as about 1.65. Check that you get the same number with our numerical method!

We're now finished with going through our basic methods. Using these rather simple methods, we'll implement successively more complicated programs for doing calculations. So far, we have not been taking in to account the interactions between electrons; next, we turn to a method that allows us to solve one-particle equations while still taking interactions in to account.

Before proceeding, I recommend playing around with these methods by using e.g. different potentials. Mastering these simple methods is key to understanding what goes on in more complicated codes.

Chapter 1.5Index Chapter 2.1