3.2 Periodic systems

Exploiting the results of the previous section, we can now consider the motion of non-interacting particles in a static potential, which is described by the time-independent Schrödinger equation 3.20. In the study of bulk crystals, the system is infinite but periodic, and so it is necessary to be able to reduce this problem to the study of a finite system. This approach turns out to have several advantages so that it is often easiest to study even aperiodic systems by imposing some false periodicity. The system is contained within a supercell which is then replicated periodically throughout space (see figure 3.1). The supercell must be large enough so that the systems contained within each one, which in reality are isolated, do not interact significantly.

Figure 3.1: Using the supercell approximation, an isolated molecule can be studied using the same techniques which are usually applied to crystals.

3.2.1 Bloch's theorem

See [45] for a fuller discussion of the proof outlined here. We consider non-interacting particles moving in a static potential $V({\bf r})$, which may be the Kohn-Sham effective potential $V_{\mathrm{KS}}({\bf r})$ (3.18). In a perfect crystal, the nuclei are arranged in a regular periodic array described by a set of Bravais lattice vectors $\{ {\bf R} \}$. The system, being infinite, is invariant under translation by any of these lattice vectors, and in particular the potential is also periodic i.e.

$$V({\bf r}+{\bf R}) = V({\bf r})$$ (3.31)

for all Bravais lattice vectors ${\bf R}$.

The Schrödinger equation which describes the motion of a single particle in this potential is

$${\hat H} \vert \psi \rangle = \left[ -{\textstyle{1 \over 2}... ...}) \right] \vert \psi \rangle = \varepsilon \vert \psi \rangle$$ (3.32)

and we define translation operators ${\hat T}_{\bf R}$ for each lattice vector ${\bf R}$ which act in the following manner on any function of position $f({\bf r})$:

$${\hat T}_{\bf R} f({\bf r}) = f({\bf r}+{\bf R}) .$$ (3.33)

Since the potential and hence the Hamiltonian are periodic i.e. ${\hat H}({\bf r}+{\bf R}) = {\hat H}({\bf r})$, these operators commute with the translation operators:

$${\hat T}_{\bf R} {\hat H}({\bf r}) \psi({\bf r}) = {\hat H}(... ...r}+{\bf R}) = {\hat H}({\bf r}) {\hat T}_{\bf R} \psi({\bf r})$$ (3.34)

i.e. $\left[ {\hat H} , {\hat T}_{\bf R} \right] = 0$, and the translation operators commute with each other i.e. ${\hat T}_{\bf R}{\hat T}_{\bf R'} = {\hat T}_{\bf R'}{\hat T}_{\bf R} = {\hat T}_{{\bf R}+{\bf R'}}$.

There must, therefore, exist a good quantum number corresponding to each lattice vector ${\bf R}$, and it must also be possible to choose the eigenstates of the Hamiltonian to be simultaneous eigenstates of all the translation operators;

$${\hat H} \vert \psi \rangle = \varepsilon \vert \psi \rangle ,$$ (3.35)

$${\hat T}_{\bf R} \vert \psi \rangle = c({\bf R}) \vert \psi \rangle .$$ (3.36)

From the commutation relations of the translation vectors it follows that the eigenvalues must satisfy

$$c({\bf R}+{\bf R'}) = c({\bf R}) c({\bf R'}) .$$ (3.37)

We can define the eigenvalues for the three primitive lattice vectors $\{ {\bf a}_i \}$ in terms of three complex numbers $\{ x_i \}$ by

$$c({\bf a}_i) = \exp ( 2 \pi {\mathrm i} x_i ) .$$ (3.38)

Since all lattice vectors can be expressed in the form ${\bf R} = n_1 {\bf a}_1 + n_2 {\bf a}_2 + n_3 {\bf a}_3$, where the $n_i$ are integers, it follows from equation 3.37 that

$$c({\bf R}) = c({\bf a}_1)^{n_1} c({\bf a}_2)^{n_2} c({\bf a}_3)^{n_3}$$ (3.39)

which is equivalent to

$$c({\bf R}) = \exp ({\mathrm i} {\bf k} \cdot {\bf R}) ,$$ (3.40)

$${\bf k} = x_1 {\bf g}_1 + x_2 {\bf g}_2 + x_3 {\bf g}_3 ,$$ (3.41)

where the $\{ {\bf g}_i \}$ are the reciprocal lattice vectors satisfying ${\bf g}_i \cdot {\bf a}_j = 2 \pi \delta_{ij}$, and the $\{ x_i \}$ are complex numbers in general.

Thus we have shown that

$${\hat T}_{\bf R} \psi({\bf r}) = \psi({\bf r}+{\bf R}) = c({... ...f r}) = \exp ({\mathrm i} {\bf k} \cdot {\bf R}) \psi({\bf r})$$ (3.42)

which is one statement of Bloch's theorem. Consider the function $u({\bf r}) = \exp(-{\mathrm i} {\bf k} \cdot {\bf r}) \psi({\bf r})$.

$$u({\bf r}+{\bf R}) = \exp(-{\mathrm i} {\bf k} \cdot [{\bf r... ...-{\mathrm i} {\bf k} \cdot {\bf r}) \psi({\bf r}) = u({\bf r})$$ (3.43)

i.e. the function $u({\bf r})$ also has the periodicity of the lattice, and so the wave-function $\psi({\bf r})$ can also be expressed as

$$\psi({\bf r}) = \exp({\mathrm i}{\bf k} \cdot {\bf r}) u({\bf r}) ,$$ (3.44)

where $u({\bf r})$ is a strictly cell-periodic function i.e. $u({\bf r}+{\bf R}) = u({\bf r})$.

We thus label the eigenstates of the Hamiltonian and the translation operators $\vert \psi_{n{\bf k}} \rangle$ where $n$ is the good quantum number labelling different eigenstates of the Hamiltonian with the same good quantum vector ${\bf k}$, related to the translational symmetry.

At this point we note that a periodic function can always be expressed as a Fourier series i.e.

$$u({\bf r}) = \sum_{\bf G} {\tilde u}_{\bf G} \exp({\mathrm i} {\bf G} \cdot {\bf r})$$ (3.45)

where ${\bf G}$ is reciprocal lattice vector ${\bf G} = m_1 {\bf g}_1 + m_2 {\bf g}_2 + m_3 {\bf g}_3$ and the $m_i$ are integers. Thus the state $\vert \psi_{n{\bf k}} \rangle$ can be expressed as a linear combination of plane-waves:

$$\psi_{n{\bf k}}({\bf r}) = \exp({\mathrm i}{\bf k} \cdot {\bf r}) u_{n{\bf k}}({\bf r})$$ (3.46)

$$= \sum_{\bf G} c_{n{\bf k}}({\bf G}) \exp[{\mathrm i}({\bf k}+{\bf G}) \cdot {\bf r}].$$ (3.47)

Instead of having to solve for a wave-function over all of (infinite) space, the problem now becomes one of solving for wave-functions only within a single (super)cell, albeit with an infinite number of possible values for ${\bf k}$. In order to simplify the problem to manageable proportions, it is necessary to impose some boundary conditions on the wave-function, which restrict the allowed values of ${\bf k}$.

3.2.2 Brillouin zone sampling

We choose to model the infinite periodic system by a large number of primitive cells $N_{\mathrm{cells}} = N_1 N_2 N_3$ stacked together, with $N_i$ cells along the ${\bf a}_i$ direction, and we apply periodic or generalised Born-von Karman boundary conditions to the wave-functions, which can be interpreted by saying that a particle which leaves one surface of the crystal simultaneously enters the crystal at the opposite surface. In fact it can be shown [46] that the choice of boundary conditions does not affect the bulk properties of the system. This condition is expressed mathematically as

$$\psi({\bf r}+N_i{\bf a}_i) = \psi({\bf r}) , \qquad i=1,2,3 .$$ (3.48)

Applying Bloch's theorem (3.42) gives

$$\psi({\bf r}+N_i{\bf a}_i) = \exp({\mathrm i} N_i {\bf k} \cdot {\bf a}_i) \psi({\bf r})$$ (3.49)

so that the values of ${\bf k}$ are restricted such that

$$\exp({\mathrm i} N_i {\bf k} \cdot {\bf a}_i) = \exp(2 \pi {\mathrm i} N_i x_i) = 1, \qquad i=1,2,3$$ (3.50)

using equation 3.41. Therefore the values of the $\{ x_i \}$ are required to be real and equal to

$$x_i = \frac{l_i}{N_i} , \qquad i=1,2,3,$$ (3.51)

where the $\{ l_i \}$ are integers, so that the general allowed form for the Bloch wave-vectors ${\bf k}$ is

$${\bf k} = \sum_{i=1}^3 \frac{l_i}{N_i} {\bf g}_i .$$ (3.52)

Taking the limit to the true infinite perfect crystal ( $N_i \rightarrow \infty$) we see that there is still an infinite number of allowed ${\bf k}$-vectors, but that they are now members of a countably infinite set. Furthermore, we see that ${\bf k}$-vectors which differ only by a reciprocal lattice vector are in fact equivalent. Consider two such wave-vectors related by ${\bf k'} = {\bf k} + {\bf G}$, then the corresponding Bloch states are also related by

$$\psi_{n{\bf k'}}({\bf r}) = \exp({\mathrm i} {\bf k'} \cdot {\bf r}) u_{n{\bf k'}}({\bf r}) =... ...) \left[ u_{n{\bf k'}}({\bf r}) \exp({\mathrm i} {\bf G} \cdot {\bf r}) \right]$$

$$= \exp({\mathrm i} {\bf k} \cdot {\bf r}) {\tilde u}({\bf r}) = \psi_{{n'}{\bf k}}({\bf r}) .$$ (3.53)

Since the expression in square brackets on the first line is a cell-periodic function the whole expression is a valid Bloch wave-function with wave-vector ${\bf k}$. Thus we can restrict our attention to those ${\bf k}$-vectors which lie within the first Brillouin zone, that volume of reciprocal-space enclosing the origin which is bounded by the planes which perpendicularly bisect lines from the origin to surrounding lattice points.

The situation now is that for each allowed ${\bf k}$-vector within the first Brillouin zone we must calculate the occupied Hamiltonian eigenstates in order to construct the density. However, the wave-functions and other properties such as Hamiltonian eigenvalues vary smoothly over the Brillouin zone [47] so that in practice only a finite set of points need to be chosen, and methods for making efficient choices have been developed [48,49,50,51,52,53]. From the calculation of the wave-functions at a certain set of ${\bf k}$-points, ${\bf k} \cdot {\bf p}$ perturbation theory [54,55] can be used to approximate the wave-functions at other nearby ${\bf k}$-points.

In this work, we are interested in the behaviour of very large systems. The volume of the Brillouin zone $\Omega_{\mathrm{BZ}}$ is related to the volume of the supercell $\Omega_{\mathrm{cell}}$ by

$$\Omega_{\mathrm{BZ}} = \frac{(2 \pi)^3}{\Omega_{\mathrm{cell}}}$$ (3.54)

so that for large systems, the Brillouin zone volume is very small and only a few ${\bf k}$-points need to be considered to describe the variation across the Brillouin zone accurately. In this work we therefore only calculate the wave-functions at the centre of the Brillouin zone, ${\bf k}=0$, known as the $\Gamma$-point. This has the added advantage that at this ${\bf k}$-point the wave-functions can be chosen to be real (recall that there is always an arbitrary global phase factor) without loss of generality.