Understand and calculate Bloch functions

In previous tutorials, we solved the Schrödinger equation for single or multiple potential wells. In solid crystalline materials, electrons move through a nearly infinite repetition of potential wells, which arise from atoms arranged in a crystal lattice. In such periodic potentials, the wave function takes the form of a Bloch function. In this tutorial, we will demonstrate how to calculate Bloch functions and explore some of their key properties.

A bit of theory

In periodic potentials, the wave function also becomes periodic and is modulated by plane wave, forming what is known as a Bloch function. Suppose the potential takes the following form:

\begin{equation} V(x) = V(x + n \cdot L) \quad \textrm{with} \quad L > 0 \quad \textrm{and} \quad n = ... -1, 0, 1, ... \end{equation}

This potential is periodic with period $ L $. According to Bloch's theorem, the corresponding Bloch function takes the form:

\begin{equation} \Psi(r) = e^{ik \cdot x} u(x) \end{equation}

Here, u(x) is a function with the same periodicity as the crystal lattice, satisfying $ u(x) = u(x + n\cdot L) $, where $ n $ is an integer and $ L $ is the lattice period. The wave vector $ k $ describes how the wave propagates through the crystal and is inversely proportional to the wavelength, given by $ \lambda = \frac{2 \pi}{|k|} $.

In essence, the wave function $\Psi$ is a plane wave ($ e^{ik \cdot x} $) which is modulated by a periodic function $ u(x) $.

In the following sections, we will use SigSpace to solve the Schrödinger equation for periodic potentials. Since the time-independent Schrödinger equation, $ {\cal H} \Psi = E \Psi $, contains no imaginary terms, its solution $ \Psi $ is typically real-valued. However, as we impose a bloch type wave with a complex exponential term $ e^{ik \cdot x} $, also $ u(x) $ and the overall solution $ \Psi $ become complex-valued.

First example

Let's calculate the wave function for a periodic quantum well. The periodicity is set to $ L = 3 $, with the well width of 0.2 centered in the periodic cell:

In [1]:
import sigspace as sg
import numpy as np

def potential(x):
    return -np.array((x >= 1.4) & (x <= 1.6), dtype=float)

H = sg.Hamiltonian(sdims=1) + potential

solver = sg.Solver(H)
solver.nodes.add(sg.BlochNode(k=1))

solver.grid(0, 3, 1000)
solver.print_info = True

sol = solver.run()

sg.quickplot(sol)
Iter: 346, energy 0.394, residual 0.007463, duration 0.088389s

As in previous examples, we define the potential using a callback function. However, instead of adding a PointNode (as in earlier tutorials), we introduce a different type of node object that instructs SigSpace to compute a Bloch function:

sg.BlochNode(k=1)

This tells SigSpace to impose a plane wave with the wave vector $ k = 1.0 $. This node type cannot be combined with others, and SigSpace automatically enforces the necessary periodic boundary conditions:

  • $ \psi(0) = \psi(a) e^{ik \cdot a} $
  • $ u(0) = u(a) $
  • $ u_x(0) = u_x(a) $

You'll notice that the resulting energy is now a complex value, although the imaginary part is very small, likely due to numerical errors. Increasing accuracy will cause the imaginary component to approach zero.

The plot displays the real and imaginary parts of the wave function within the periodic cell, from $ x=0 $ to $ 3 $. The wave function is not strictly periodic, consistent with Bloch's theorem, which states that it is a periodic function ($u(x)$) modulated by a plane wave $e^{ik \cdot x}$.

This becomes evident when we extrapolate the wave function over a larger spatial region:

In [2]:
x = np.linspace(-100, 100, 1000)

sg.quickplot((x, sol.wave(x)))

The wave function is typically obtained by calling sol.wave(). If you provide an array of coordinates as an argument, SigSpace will attempt to extrapolate the wave function. Thanks to the periodic nature of Bloch functions, this is entirely feasible, as the periodic cell contains all necessary information.

At larger scales, the modulation of the periodic function $ u(x) $ by the plane wave $ e^{ikx} $ becomes clearly visible.

Based on the definition given at the beginning of this tutorial, it is also possible to obtain the periodic function $ u(x) $ itself. Just divide the wave function by the term of the plane wave:

In [3]:
lx = sol.wave.x
k = sol.wave.nodes.k
phi = sol.wave()[:,0]

u = phi/np.exp(1j * k * lx)

sg.quickplot((lx, u), layout={"xlabel": "x", "ylabel": "u"})