An example: Kronig-Penney

The Kronig-Penney model is a favorite of physics instructors across the world, because it nicely illustrates some quantum mechanical phenomena. The basic idea of the original model is the following: Suppose you want to model the behavior of a solid, which is basically an infinitely repeating arrangement of atoms. To simplify this 3-dimensional problem, we tackle the one-dimensional case. In this case, the potential is basically just a bunch of "potential peaks"\ evenly spaced, as shown in the following figure.

Kronig-Penney
Kronig-Penney potential.

Suppose we simulate this on a computer by repeating this structure something like 10 times: we put 10 potential peaks in some finite area, say, of length 10. These peaks have a particular length and width. If they're supposed to simulate atoms, we can say that the peaks are rather sharp and their height is great; for our purposes, let's set the height of the potential at \( V_0 = 100 \) Hartrees and the width at, say, 0.2 Bohrs.

Using the knowledge of the previous two sections, what should we do? I recommend that you actually think about this for a while before looking at the code. On the diagonal, we want the potential, which is zero in most places but 100 in others. Thus, our potential will be an array with as many points as there are points on the grid, with certain elements being equal to 100 and the rest to zero. The kinetic energy operator remains unchanged. Hence, the code we previously used is just slightly modified:


import numpy as np
import matplotlib.pyplot as plt

V0 = 100
npoints = 1000
V = np.zeros(npoints)
x = np.linspace(0,10,npoints)
dx = x[1]-x[0]
V[40:61] = V0 
V[140:161] = V0
V[240:261] = V0
V[340:361] = V0
V[440:461] = V0
V[540:561] = V0
V[640:661] = V0
V[740:761] = V0
V[840:861] = V0
V[940:961] = V0
kinetic_energy = -0.5/dx**2 \
                * (np.diagflat(-2*np.ones(1000))\
                +np.diagflat(np.ones(999),-1)\
                +np.diagflat(np.ones(999),1))
H = kinetic_energy + np.diagflat(V)

e,vec=np.linalg.eigh(H)


plt.rcParams.update({'font.size': 15})
plt.xlabel("k*L")
plt.ylabel("Energy")
plt.scatter(range(30),e[0:30]/np.pi**2)
plt.show()

And the output should be something like this:
Kronig-Penney band gaps
Kronig-Penney band gaps.

The gaps between the energy levels, where there are many close to each other after which there is a jump, are called "band gaps". The reason Kronig-Penney is used in learning is precisely because of these gaps: they explain why there are insulators. Electrons conduct electricity. In insulators, there is a band gap such that it is difficult for the electrons from the lower energy levels to jump to the higher ones to conduct electricity, as they need lots of energy to do it. In conductors, by contrast, this band gap is either very small or doesn't exist. In semiconductors, there is a gap and whether the electrons can make the jump depends on the conditions.

An interesting thing to test, by the way, is to add more potential peaks and see what happens (just make sure they're evenly spaced). What you should see is that the points get closer together. In the limit with \( N_{peaks}\rightarrow \infty\) and \( L\rightarrow \infty \), we get basically continuous strings of energy levels, plus the gaps. That's why they're called "band gaps" instead of "energy level gaps".

So far, we've been using a grid. That is a strategy that often works in computational problems. Nevertheless, we will next cover an important technique that is used in virtually every real electronic structure package: basis functions.

Chapter 1.3Index Chapter 1.5