PHY.K02UF Molecular and Solid State Physics

Plane wave method

Consider an electron moving in a periodic potential. The Schrödinger equation for this case is,

$$-\dfrac{\hbar^2}{2m} \, \nabla^2 \psi_{\vec{k}} + U(\vec{r}) \psi_{\vec{k}} = E \psi_{\vec{k}}. $$

Since the potential $U(\vec{r})$ is periodic, it can be written as Fourier series,

$$U(\vec{r}) = \sum_{\vec{G}} U_{\vec{G}} e^{i \vec{G} \cdot \vec{r}},$$

Felix Bloch showed that the solutions to the Schrödinger equation for an electron in a periodic potential have the form,

$$\psi_{\vec{k}}(\vec{r}) = e^{i\vec{k}\cdot\vec{r}}u_{\vec{k}}(\vec{r}) = e^{i\vec{k}\cdot\vec{r}}\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}},$$

where the periodic function $u_{\vec{k}}=\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}}$ has been expressed as a Fourier series. The reciprocal lattice vectors have been relabeled as running over $\vec{G}'$ instead of $\vec{G}$. It does not matter how they are labeled since we sum over all of the reciprocal lattice vectors, $\sum\limits_{\vec{G}}C_{\vec{G}}e^{i\vec{G}\cdot\vec{r}}=\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}}$. This form for the wavefunction is substituted into the Schrödinger equation,

$$\sum\limits_{\vec{G}'}\dfrac{\hbar^2(\vec{k}+\vec{G}')^2}{2m} C_{\vec{G}'}e^{i(\vec{k}+\vec{G}')\cdot\vec{r}} + \sum_{\vec{G}}\sum\limits_{\vec{G}''} U_{\vec{G}} C_{\vec{G}''}e^{i(\vec{k}+\vec{G}+\vec{G}'')\cdot\vec{r}} = E \sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i(\vec{k}+\vec{G}')\cdot\vec{r}}. $$

In the middle term with the double sum, the sum over $\vec{G}'$ has been relabeled as a sum over $\vec{G}''$. Again it does not matter that the label has changed since the sum is over all of the states but we need a way to keep track of the product terms in the double sum. Next we collect like terms. The exponential factors can be written as $e^{i \vec{k} \cdot \vec{r}}= \cos(\vec{k} \cdot \vec{r})+ i\sin(\vec{k} \cdot \vec{r})$. Only terms that have the same wavelength can be equal to each other so only the terms where $\vec{k} +\vec{G}'= \vec{k} + \vec{G}+\vec{G}''$ can be equal to each other. This results in the condition $\vec{G}'' = \vec{G}' - \vec{G}$ and for each $\vec{G}'$ there is an equation,

$$\dfrac{\hbar^2(\vec{k}+\vec{G}')^2}{2m}C_{\vec{G}'} + \sum_{\vec{G}} U_{\vec{G}} C_{\vec{G}'-\vec{G}} = E C_{\vec{G}'}. $$

This set of equations are called the central equations. The Schrödinger equation, a differential equation for $\psi_{\vec{k}}$, has been replaced with the central equations which are algebraic equations for the coefficients $C_{\vec{G}'}$. The algebraic equations can be put in the form of an eigenvalue problem. These equations written out in the one-dimensional case for $-G_0, 0, G_0$ are,

$$\dfrac{\hbar^2(k-G_0)^2}{2m}C_{-G_0} + (\cdots + U_{-2G_0}C_{G_0} + U_{-G_0}C_{0} + U_{0}C_{-G_0} + \cdots ) = E C_{-G_0}. $$ $$\dfrac{\hbar^2k^2}{2m}C_{0} + (\cdots + U_{-G_0}C_{G_0} + U_{0}C_{0} + U_{G_0}C_{-G_0} + \cdots ) = E C_{0}. $$ $$\dfrac{\hbar^2(k+G_0)^2}{2m}C_{G_0} + (\cdots + U_{0}C_{G_0} + U_{G_0}C_{0} + U_{2G_0}C_{-G_0} + \cdots ) = E C_{G_0}. $$

If the terms indicated by $\cdots$ are negelected, this can be written in matrix form as an eigenvalue problem that can be solved for the energies.

$$\textbf{M} \vec{C} = E \vec{C}.$$

When more equations for $\vec{G}'$ are included, the matrix gets bigger and since the coefficients $U_{\vec{G}}$ typically get smaller as $\vec{G}$ gets larger, it becomes more justified to neglect the terms indicated by $\cdots$. An $N\times N$ matrix will have $N$ solutions for the energy $E$ at each value of $\vec{k}$. The different solutions correspond to different bands.

In three dimensions the reciprocal lattice vectors are typically indexed by three integers $h$, $k$, and $l$. However to construct the matrix, it is necessary to index the reciprocal lattice vectors by a single integer. This can be done by making a list of $[hkl]$ triples and assigning each triple to an integer $i$. The python code below defines $h$, $k$, and $l$ vectors and assigns them to an integer $i$ like this:

i = 0
for nx in range(5):  
    for ny in range(5):
        for nz in range(5):
            h[i] = nx-2 # h = -2,-1,0,1,2
            k[i] = ny-2 # k = -2,-1,0,1,2
            l[i] = nz-2 # l = -2,-1,0,1,2
            i= i+1 

The periodic potential could be a periodic Coulomb potential, a muffin tin potential, or a pseudopotential. To calculate the dispersion relation for a simple cube, bcc, or fcc lattice with a periodic potential not listed below, use the code for the Coulomb potential and it will only be necessary to change the line that specifies the off-diagonal elements of the matrix $\textbf{M}$ for that potential.

Comparisons with DFT calculations: Li, Na, Mg, Al, Cu.


1-D zero potential


1-D comb potential


1-D square-wave potential


1-D cosine potential


Body Centered Cubic with a zero potential


Body Centered Cubic with a Comb Potential


Body Centered Cubic with a Coulomb Potential


Face Centered Cubic with a zero potential


Face Centered Cubic with a Comb Potential


Face Centered Cubic with a Coulomb Potential


Simple Cubic with a zero potential


Simple Cubic with a Comb Potential


Simple Cubic with a Coulomb Potential