Calculate Brillouin diagram for a periodic potential

The dispersion relation links the wavelength of a particle to its frequency (or energy). This concept is particularly useful when examining the electrical properties of crystalline solid-state materials. In these materials, the periodic potential allows wave propagation for most energy levels. However, certain energy levels are prohibited, leading to energy gaps where no states are allowed states exist. These gaps are referred to as bandgaps. The arrangement of allowed and forbidden energy ranges, known as the band structure, determines wether a material behaves as an insulator), metal, or semiconductor. Consequently, the study of quantum states in periodic potentials is of great importance and has been extensively researched for many materials.

In this tutorial, we will calculate the dispersion relation for a one-dimensional potential. We begin by considering the free electron, introduce the concept of the Brillouin zone, and then explore how periodic potentials modify the energy levels and quantum states.

Before proceeding, please ensure you have read the tutorial on bloch functions, as they form the foundation for understanding this tutorial.

The free particle

Let's start with a simple example. We calculate the dispersion relation for a free electron:

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

La = 3          # Width of solution space

k_list = np.linspace(0, 5, 200)
E_list = np.array([])

for k in k_list:

    H = sg.Hamiltonian(sdims=1)
    
    solver = sg.Solver(H)
    solver.nodes.add(sg.BlochNode(k=k))
    solver.grid(0, La, 100)
    
    sol = solver.run()
    E_list = np.append(E_list, sol.energy.real)

sg.quickplot((k_list, E_list), layout={"xlabel": "Wave number [1/a0]", "ylabel": "Energy [Ha]"})

Let's walk through the code step by step:

  • In the dispersion relation, the energy of allowed states is plotted as a function of the wave number k. First, we define k_list, which contains the wave numbers we wish to calculate, and E_list to store the corresponding energies.The calculation is done in one dimension for a periodic unit cell extending from $ 0 $ to $ La $, where La defines the width of the solution space.
  • For each wave number k, we calculate the corresponding energy. To do this, we start by creating an empty Hamiltonian. As shown in the tutorial on Bloch functions, we add a BlochNode with the appropriate wave number, solve the system via solver.run(), and store the resulting energy in 'E_list'.
  • The quickplot command can also take a list of arrays for easy plotting.

The result is a quadratic curve, which matches the expected relation $ E = 0.5 \cdot k^2$, representing the energy of a free electron.

The Brillouin zone

In the example above, we calculated the full dispersion relation of a free particle. However, in most literature, dispersion relations are typically presented within the first Brillouin zone. Although the resulting graphs may appear different, they represent a different visualization of the same curve.

The Brillouin zone is the unit cell in reciprocal space, ranging from $ K_- = -\pi/La $ to $K_+ = +\pi/La $, where $La$ is the lattice constant. The boundary of this zone at $ K_\pm $ corresponds to the first standing wave in the periodic lattice. Its significance lies in the fact that, within a periodic potential, all wave behavior can be fully described within this first Brillouin zone. This makes it an invaluable tool to gain an overview of all possible wave modes that can propagate in the cristal.

In the graph, the dispersion relation (shown above) must be mapped to the first Brillouin zone within the range of $ K_\pm $. This is done by calculating $ k' = k + n \cdot 2\pi/a $, where $n$ is an integer. For instance, a wave with $k$ and $k′$ values that differ by $ 2\pi/a $ will occupy the same position in reciprocal space, though they may have different energy levels

SigSpace provides a function that simplifies plotting dispersion relations. It handles both the mapping and the calculation of the Brillouin zone boundaries, making visualization easier:

In [2]:
import sigspace.brillouin as br
br.plot_bzone(k_list, E_list, La, mirror_k=True)

The first two parameters in plot_bzone(k_list, E_list, La, mirror_k=True) represent the dispersion relation curve we previously calculated. The third parameter, La, specifies the width of the unit cell and is used to determine the location of the Brillouin zone boundaries. mirror_k is set to True to mirror the dispersion relation for negative wave numbers, reducing computational effort since the relation is symmetric.

The resulting graph displays three lines:

  • Dispersion relation: These are the mapped values from (k_list, E_list)
  • Extrapolated: These lines extends from the last calculated point to the zone boundary, helping to visualize potential bandgaps. If you zoom in on the intersections at the zone boundaries (the graph is interactive, just select the area you want to enlarge), you'll notice the extrapolated lines touch each other. As expected for a free electron, no forbidden energy levels are present.
  • Zone boundary: This line marks the edges of the first Brillouin zone

The energy levels shift when a potential is added to the Hamiltonian:

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

La = 3

k_list = np.linspace(0, 5, 200)
E_list = np.array([])

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

for k in k_list:

    H = sg.Hamiltonian(sdims=1) + potential
    
    solver = sg.Solver(H)
    solver.nodes.add(sg.BlochNode(k=k))
    
    solver.grid(0, La, 1000)
    solver.max_iterations = 1000
    
    sol = solver.run()

    if sol.residual < 1e-2:
        E_list = np.append(E_list, sol.energy.real)
    else:
        E_list = np.append(E_list, float("nan"))


sg.quickplot((k_list, E_list), layout={"xlabel": "Wave number [1/a0]", "ylabel": "Energy [Ha]"})

Several changes are made to the code compared to the free electron case:

  • We define a potential(x) callback function and add it to the Hamiltonian. This function is evaluated within the periodic unit cell.
  • The maximum number of iterations is limited to 1000. For every result we check if the solver has converged before adding the result. Since certain energy levels may be forbidden, the solver's convergence can significantly decrease near these forbidden energy levels. Without a limit on iterations, the calculation would enter an infinite loop. Additionally, to avoid including incorrect results in the dispersion relation, we insert a NaN where the solver fails to converge, which helps visualize these issues.

The resulting graph resembles the one from the free electron case, but now steps and gaps appear. While some gaps are due to solver limitations, noticeable steps occur near $ k = \pi/L_a = 1.04 $ or $ k = 2\pi/L_a = 2.09 $. These steps create energy gaps at the zone boundaries.

The appearance of these forbidden energies becomes even clearer when plotting the Brillouin zone:

In [4]:
import sigspace.brillouin as br

br.plot_bzone(k_list, E_list, La, mirror_k=True)

The commands for plotting are identical to those used in the free electron case.

However, when you zoom in on the zone boundaries or at the center, you'll notice small gaps forming. These are the bandgaps where no states exist for wave propagation. For example, the lowest bandgap at $ k = \pm\pi/a$ is approximately 113 mHa.

To understand the origin of these gaps, we perform another calculation. We solve the Hamiltonian for the same potential just above and below the zone boundary and plot the results:

In [5]:
import sigspace as sg
import numpy as np
import math

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

La = 3

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

solver = sg.Solver(H)
solver.nodes.add(sg.BlochNode(k=math.pi/La-0.1))
solver.grid(0, La, 1000)
solver.print_info = True
sol1 = solver.run()

solver = sg.Solver(H)
solver.nodes.add(sg.BlochNode(k=math.pi/La+0.1))
solver.grid(0, La, 1000)
solver.print_info = True
sol2 = solver.run()

x = sol1.wave.x
sg.quickplot((x, sol1.wave.density()), (x, sol2.wave.density()), (x, potential(x)), 
             layout={"xlabel": "x", "ylabel": "|Ψ|^2", "legend": ["k1=𝛱/a-0.1", "k2=𝛱/a+0.1", "V"]})
Iter: 215, energy 0.495, residual 0.007857, duration 0.021422s
Iter: 268, energy 0.741, residual 0.006903, duration 0.026508s

The figure shows the absolute squared wave function for $k_1$ and $k_2$. Despite the small change in the wave number, there is a jump in the wave function. As a result, the wave with $k_2$ has a much larger overlap with the potential function $V(x)$ compared to the wave with $k_1$. This increased overlap affects the energy of the wave function, causing a sudden shift in energy. Since no other propagation modes are available at this energy level, no electrons can exist here.

As a final note, Bloch waves also support operators. The following example demonstrates how the momentum can be obtained. Because a specific wave number is defined, the momentum is non-zero, in contrast to the static solutions in previous examples.

In [6]:
from sigspace.H.braket import *

solver = sg.Solver(H)
solver.nodes.add(sg.BlochNode(k=3))
solver.grid(0, La, 1000)
sol1 = solver.run()

psi = sol1.state

impulse = bra <psi|P|psi> ket
impulse
Out[6]:
np.complex128(2.9693418181916953+3.9847826866896974e-05j)