Numerical Integration over Brillouin Zones: The idea behind Monkhorst-Pack

Numerical integration is the art of estimating integrals over domains by sampling a function at a set of points. Symbolically, this is \begin{align*} \int_{\Omega} f(\mathbf{x}) \, \mathrm{d}\mathbf{x} \approx \sum_{i = 0}^{n-1} w_{i}f(\mathbf{x}_i) \end{align*} where \(\{\mathbf{x}_{i}\}_{i=0}^{n-1} \subset \Omega\) are called quadrature nodes and \(\{w_{i}\}_{i=0}^{n-1}\) are called weights.

The Brillouin zone of the FCC lattice, shown with the quadrature nodes of the Monkhorst-Pack scheme. The reciprocal lattice vectors \(\mathbf{b}_1, \mathbf{b}_2\) and \(\mathbf{b}_3\) are shown in yellow.

If \(f \colon \mathbb{R} \to \mathbb{R}\), the question of how to numerically integrate \(f\) is largely settled, though the rules are complicated:

  • If \(f\) is periodic and you are integrating it over a period, use trapezoidal quadrature.
  • If \(f\) is holomorphic on the interior of a disc encompassing the interval of integration and you have high accuracy requirements, use tanh-sinh quadrature.
  • If you require a deterministic runtime and need an error estimate, use Gauss-Kronrod integration.
  • If you require a deterministic runtime and do not need an error estimate, use Gaussian quadrature.
The convergence rate of all these methods is dependent on the smoothness of \(f\). For example, if \(f\) is \(m\)-times differentiable, the error trapezoidal quadature on a periodic function is a consequence of the Euler-Maclaurin summation formula \begin{align*} T_{h}[f] := h\sum_{k=0}^{N-1} f(a+kh) = \int_{a}^{b} f(x) \, \mathrm{d}x - (-1)^{m}h^{m} \int_{a}^{b} \tilde{B}_{m}\left(\frac{x-a}{h} \right)f^{(m)}(x) \, \mathrm{d}x \end{align*} where \(h := \frac{b-a}{n}\) and \(\tilde{B}\) is the periodic extension of the \(m\)th Bernoulli polynomial. For infinitely differentiable periodic functions, this trapezoidal quadrature is fantastic. But if \(f\) is discontinuous at (say) a point \(c \in (a, b)\) the convergence is spoiled and the most reasonable approach is to split the integral at the discontinuity: \begin{align*} \int_a^b f(x) \, \mathrm{d}x = \int_a^c f(x) \, \mathrm{d}x + \int_c^b f(x) \, \mathrm{d}x \end{align*} Both these ideas will come in handy while designing quadrature in higher dimensions.

Desirable properties of quadrature formulae are

  • Error estimation: It should be possible to specify an accuracy goal rather than the number of function evaluations.
  • Nesting: Producing an error estimate requires calculation of the integral on a coarse set of nodes, and a fine set nodes-the difference of these estimates being the error estimate. If the coarse set of nodes is a subset of the fine set, then the quadrature rule is said to nest, resulting in a reduction of computational effort.
  • Interior nodes: The nodes should be strictly contained in \(]a, b[\). This is because we may wish to compute an integral where the integrand is singular on the boundary. If the quadrature nodes lie on the boundary, then this is impossible.
  • Positive weights: If a weight is negative, then we can construct a strictly positive function which the quadrature rule evaluates to a negative number.
All these goals generalize to multiple dimensions-though they are harder to fulfill.

In 1D, we can classify the domains into three types: Bounded, half-infinite, and infinite-the latter two being reducible to a bounded interval via variable transformations. In higher dimensions, the number of types of domains is infinite. For instance, Gaussian quadrature schemes in 2D must be worked out for the circle, the triangle, the square, the hexagon, and so on. Even on a single domain (say, the circle) the quadrature rules are far from unique-and entire encyclopedias have been made of the various ways in which node and weight choice can be made.

In solid state physics, we are interested in numerical integration of functions over Brillouin zones. For metals, these functions are discontinuous at the Fermi surface, which can have cusps and enclose disjoint regions. For insulators, this situation is better as the integrands are at least continuous. But before entering the veritable zoo of quadrature methods common in solid state physics, let's first try to understand our use case.

From basic quantum mechanics, we know that if an electron is in a state \(\psi\), then the expectation of an observable \(\left< X \right>\) is \begin{align*} \left< X \right> = \int \psi^{*}X\psi \, \mathrm{d} \mathbf{x} \end{align*} If \(\psi\) is a Bloch state with crystal momentum \(\mathbf{k}\) and band index \(n\), this becomes \begin{align*} \left< X_{n}(\mathbf{k}) \right> &= \int_{\mathrm{crystal}} \exp(-i\mathbf{k}\cdot \mathbf{x})u_{n,\mathbf{k}}(\mathbf{x}) X \exp(i\mathbf{k}\cdot \mathbf{x})u_{n,\mathbf{k}}(\mathbf{x}) \, \mathrm{d} \mathbf{x} \\ &= N_{\mathrm{cells}}\int_{\mathrm{unit\, cell}} \exp(-i\mathbf{k}\cdot \mathbf{x})u_{n,\mathbf{k}}(\mathbf{x}) X \exp(i\mathbf{k}\cdot \mathbf{x})u_{n,\mathbf{k}}(\mathbf{x}) \, \mathrm{d} \mathbf{x} \\ \end{align*} A couple examples are \(\left< \mathcal{H}_{\mathbf{k}} \right> = E_{n}(\mathbf{k})\) and \begin{align*} \left< \mathbf{v}_n(\mathbf{k}) \right> &= \frac{1}{\hbar} \nabla_{\mathbf{k}} E_{n}(\mathbf{k}) = \nabla_{\mathbf{k}} E_{n}(\mathbf{k}) \, (\mathrm{a.u.}) \end{align*}

However, experimental methods struggle to detect the state of a single electron in a crystal, and often we must be content to observe bulk properties. For example, the conductivity of copper is not due to a single electron passing through the crystal, it is the cumulative movement of all the electrons. Hence to calculate a bulk property, we sum over all occupied Bloch states. For instance, the energy per unit volume of a crystal is \begin{align*} \frac{2}{\Omega_{\mathrm{crystal}}} \sum_{n}\sum_{\mathbf{k} \in \mathbf{BZ}} \mathbf{E}_{n}(\mathbf{k})f(E_{n}(\mathbf{k})) &= \frac{1}{4\pi^3} \sum_{n}\sum_{\mathbf{k} \in \mathbf{BZ}} \mathbf{E}_{n}(\mathbf{k})f(E_{n}(\mathbf{k})) \Delta \mathbf{k} \to \frac{1}{4\pi^3} \sum_{n}\int_{\mathbf{BZ}} \mathbf{E}_{n}(\mathbf{k})f(E_{n}(\mathbf{k})) \mathrm{d} \mathbf{k} \end{align*} where \begin{align*} f(E) := \frac{1}{\exp(\beta(E-\mu)) + 1} \end{align*} is the Fermi function which chemical potential \(\mu\). Let's do a few more examples:

Electric current density: An electron in a Bloch state carries a current \(\mathbf{j}_{e} = -e\mathbf{v}_n(\mathbf{k})\), hence the total current density in a band is \begin{align*} \mathbf{j}_{n} = -\frac{e}{4\pi^3\hbar} \int_{\mathrm{BZ}} \nabla_{\mathbf{k}}E_{n}(\mathbf{k}) f(E_{n}(\mathbf{k})) \, \mathrm{d}\mathbf{k} \end{align*}

Energy current density: An electron in a Bloch state carries energy current \(E_{n}(\mathbf{k})v_{n}(\mathbf{k})\), so the total energy current is \begin{align*} \mathbf{j}_{n,E} = \frac{1}{4\pi^3\hbar} \int_{\mathrm{BZ}} E_{n}(\mathbf{k})\mathbf{v}(\mathbf{k}) f(E(\mathbf{k})) \, \mathrm{d}\mathbf{k} \end{align*}

Semiclassical conductivity in the relaxation time approximation: The conductivity tensor relates the current density to an electric field via \(\mathbf{j} = \mathbf{\sigma}\mathbf{E}\). With a large number of hypotheses and approximations (discussed in Ashcroft and Mermin) we can write \begin{align*} \sigma = \frac{e^2 }{4\pi^3\hbar}\sum_{n} \tau_n(E_F)\int_{\mathrm{BZ}} \nabla_{\mathbf{k}}\mathbf{v}_{n}(\mathbf{k})f(E(\mathbf{k})) \, \mathrm{d}\mathbf{k} \end{align*} where \(\frac{1}{\hbar}\nabla_{\mathbf{k}}\mathbf{v}_{n}\) is the inverse of the effective mass tensor and \(\tau_n(E_F)\) is the scattering time in the \(n\)th band at the Fermi level.

Other examples of properties computable from integrals over Brillouin zones are AC electric conductivity, thermal current density, the Lindhard dielectric constant, and constant volume specific heats.

All our examples are of the form \begin{align*} \int_{\mathrm{BZ}} \phi(\mathbf{k}) \, \mathrm{d}\mathbf{k} \end{align*} where \(\phi(\mathbf{k}+ \mathbf{K}) = \phi(\mathbf{k})\) for all reciprocal lattice vectors \(\mathbf{K}\). Because \(\phi\) is periodic, we might suspect that like trapezoidal quadrature, equispaced sampling over the Brillouin zone should converge rapidly. However, to my knowledge, there is no extension of the Euler-Maclaurin summation formula to multiple dimensions, so that route to proving the intuition seems closed. But using Fourier analysis, we can get the same result as we obtained with Euler-Maclaurin, and this proof will extend to higher dimensions. Suppose \(f\colon \mathbb{R}\to \mathbb{R}\) is \(T\) periodic with a convergent Fourier series. Then \begin{align*} f(x) = \sum_{m \in \mathbb{Z}} c_{m}\exp(2\pi i x/T) \end{align*} and \begin{align*} \int_{0}^{T} f(x) \, \mathrm{d}x = Tc_{0} \end{align*} The trapezoidal rule applied to \(f\) is \begin{align*} T_{h}[f] = \sum_{k=0}^{N-1} f(kh)h \end{align*} where \(h:= T/N\). Plugging in the Fourier series representation of \(f\) into this gives \begin{align*} T_{h}[f] &= \sum_{m \in \mathbb{Z}} \frac{Tc_{m}}{N} \sum_{k=0}^{N-1} \exp(2\pi i km/N) \\ &= T\sum_{j \in \mathbb{Z}} c_{jN} \\ &= \int_{0}^{T} f(x) \, \mathrm{d}x + T\sum_{j \in \mathbb{Z}, j\ne 0} c_{jN} \end{align*} so the error of the trapezoidal rule is \begin{align*} E_{N} := \left|T_{h}[f] - \int_{0}^{T} f(x) \, \mathrm{d}x \right| = T\left| \sum_{j \in \mathbb{Z}, j\ne 0} c_{jN} \right| \le T\sum_{j=1}^{\infty}\left|c_{jN} + c_{-jN} \right| \end{align*} Now we just need to know how fast the coefficients of the Fourier series decays. There are lots of really tight error bounds known for Fourier series, but the following estimate will suffice: If \(f\) is \(p\)-times continuously differentiable, then \begin{align*} |c_{n}| \le \frac{\left\|f^{(p)} \right\|_{L_1}}{|n|^{p}} \end{align*} and hence \begin{align*} E_{N} \le \frac{2\left\|f^{(p)} \right\|_{L_1}}{N^{p}} \sum_{j=1}^{\infty} \frac{1}{j^{p}} = \frac{2\left\|f^{(p)} \right\|_{L_1}}{N^{p}}\zeta(p) \end{align*} In higher dimensions, things are a bit more tedious, but not fundamentally different. We write \begin{align*} \phi(\mathbf{k}) = \sum_{\mathfrak{t}\in \mathrm{Bravais}} c_{\mathfrak{t}}\exp(i\mathbf{k}\cdot \mathfrak{t}) \end{align*} where \(\mathfrak{t} = n_1\mathbf{a}_1 + n_2\mathbf{a}_2 + n_3\mathbf{a}_3\) is a Bravais lattice vector. With some suffering we can extend the argument of the 1D case to 3D, with one difference: If it takes us \(N\) function evaluations to reach a certain accuracy in 1D, then it will take us roughly \(N^3\) evaluations to reach the same accuracy in 3D. This is a very general idea in quadrature-going roughly under the heading of the curse of dimensionality-but nonetheless, if \(\phi\) is very smooth (\(p\) is large), then \(\mathcal{O}(1/N^{p/3})\) convergence is still quite fast. The quadrature nodes of trapezoidal quadrature over a primitive cell in reciprocal space can be taken as \begin{align*} \mathbf{k}_{ijk} = \frac{2i - N + 1}{2N}\mathbf{b}_1 + \frac{2j - N + 1}{2N}\mathbf{b}_2 + \frac{2k - N + 1}{2N}\mathbf{b}_3 \end{align*} for \(i,j,k \in \{0, 1, 2,\ldots, N-1\}\). The quadrature weights are all equal and given by \(\Omega_{\mathrm{BZ}}/N^3\). This generalization of trapezoidal quadrature to \(k\)-space integrals is called Monkhorst-Pack quadrature. However, although the quadrature nodes fill a primitive cell in reciprocal space (and hence should converge to the same value as an integral taken over the Brillouin zone), they are not strictly contained in the Brillouin zone, as can be seen in the figure. This can cause some suffering if you absolutely require the nodes to be contained in the BZ-some method of reducing points out of the BZ into it via reciprocal lattice vectors is required. In addition, the method doesn't nest as nicely as 1D trapezoidal quadrature. With 1D trapezoidal quadrature, if we evaluate \(T_{h}[f]\), then half the nodes of \(T_{h/2}[f]\) have already been computed. With 3D Monkhorst-Pack quadrature, halving the spacing (i.e, \(N \mapsto 2N\)) doesn't leave any quadrature nodes in common. We actually have to cut the node spacing by a factor of three to get nesting, and then only \(1/3^3\) nodes have already been computed. So error estimation via nesting is much more time consuming for Monkhorst-Pack integration than for 1D trapezoidal quadrature.

In solid state physics, many of the integrals that interest us require calculation of the band structure \(E_{n}(\mathbf{k})\). Calculation of the band structure requires solution of a PDE, so function evaluations over the BZ are expensive. However, we have already learned that many interesting crystals have a large point group in addition to translational symmetry. If \(\phi\) is invariant under the point group of the crystal, then \begin{align*} \int_{\mathrm{BZ}} \phi(\mathbf{k}) \, \mathrm{d}\mathbf{k} = |\mathbf{P}|\int_{\mathrm{IBZ}} \phi(\mathbf{k}) \, \mathrm{d}\mathbf{k} \end{align*} where \(|\mathbf{P}|\) is the cardinality of the point group and \(\mathbf{IBZ}\) is the irreducible Brillouin zone. This reduces the number of function evaluations by a factor of \(|\mathbf{P}|\), but requires we characterize the IBZ in a way we can readily filter the Monkhorst-Pack nodes contained in the IBZ. Since there are 73 topologically inequivalent Brillouin zones, this is a task of some difficulty.

Monkhorst-Pack quadrature is unlikely to be improved on for smooth periodic functions. However, the Fermi surface of metals represents a discontinuity which spoils its convergence. In 1D, our strategy for integrating discontinuous functions was to break the integral into pieces where the integrand was continuous. Similarly, in 3D, the traditional method of integrating functions which are discontinuous at the Fermi surface is to tesselate the interior(s) of the Fermi surface with tetrahedra. The function of interest is evaluated at the corners of the tetrahedron, and linearly interpolated inside it. Since there are analytic formulae for integrals of linear functions over tetrahedra, we can evaluate the integral on each tetrahedron rapidly, and the entire integral is the sum of the integrals over the tetrahedra. This "linear tetrahedron integration" is easy-once the Fermi surface has been elucidated and the tesselation has been computed! But tesselation and Fermi surface characterization is a task of sufficient difficulty that we will analyze it another day.

  • Cools, R. & Kim, K.J. Korean J. Comput. & Appl. Math. (2000) 7: 477.
  • Blöchl, Peter E., Ove Jepsen, and Ole Krogh Andersen. "Improved tetrahedron method for Brillouin-zone integrations." Physical Review B 49.23 (1994): 16223.
  • Methfessel, M. P. A. T., and A. T. Paxton. "High-precision sampling for Brillouin-zone integration in metals." Physical Review B 40.6 (1989): 3616.
  • Monkhorst, Hendrik J., and James D. Pack. "Special points for Brillouin-zone integrations." Physical review B 13.12 (1976): 5188.
  • Ashcroft, Neil W., N. David Mermin, and Sergio Rodriguez. "Solid state physics." (1998).

Irreducible Brillouin Zones and Band Structures

In Bloch's Theorem via representation theory we discussed how the discrete translational symmetry of a crystal leads to solutions of Schrodinger's equations which are highly structured, in the sense that translations of the wavefunctions can be computed via phase shifts. The goal of this post is to obtain a greater understanding of Brillouin zones and band structures, but first we need to review a few aspects of Bloch's theorem that weren't relevant to the aforementioned post. Our fundamental conclusion was that if \(\mathcal{H}\psi = E\psi\) then \(\psi(\mathbf{r} + \mathfrak{t}) = \exp(i\mathbf{k}\cdot\mathfrak{t})\psi(\mathbf{r})\) where \(\mathbf{k} := \frac{n_1}{N_1}\mathbf{b}_1 + \frac{n_2}{N_2}\mathbf{b}_2 + \frac{n_3}{N_3}\mathbf{b}_3\) and \(\mathfrak{t} := m_1\mathbf{a}_1 + m_2\mathbf{a}_2 + m_3\mathbf{a}_3\) is a translation vector of the crystal. Let \(u_{\mathbf{k}}(\mathbf{r}) := \exp(-i\mathbf{k}\cdot \mathbf{r})\psi(\mathbf{r})\). A quick calculation shows that \(u_{\mathbf{k}}\) has the same periodicity of the lattice, i.e., \(u_{\mathbf{k}}(\mathbf{r} + \mathfrak{t}) = u_{\mathbf{k}}(\mathbf{r})\). By turning this around, we see that \(|\psi|^{2}\) takes the same value at any two points that differs by a lattice vector, and hence the probability of an electron in a given unit cell is the same as it being in any other unit cell. Hence we can use translation symmetry to reduce computation. To wit, our original goal was to solve \begin{align}\label{periodic_schrod} -\frac{1}{2}\nabla^2 \psi + V\psi = E\psi, \end{align} within the entire crystal. But since we know that \(\psi(\mathbf{r}) = e^{i\mathbf{k}\cdot \mathbf{r}}u_{\mathbf{k}}(\mathbf{r})\), we have \begin{align*} -\frac{1}{2}\nabla^2 \psi &= -\frac{1}{2}\nabla\cdot\left(i\mathbf{k}e^{i\mathbf{k}\cdot \mathbf{r}}u_{\mathbf{k}}(\mathbf{r}) + e^{i\mathbf{k}\cdot \mathbf{r}} \nabla u_{\mathbf{k}} \right) \\ &= \frac{1}{2}\left(k^2e^{i\mathbf{k}\cdot \mathbf{r}}u_{\mathbf{k}}(\mathbf{r}) - 2ie^{i\mathbf{k}\cdot \mathbf{r}}\mathbf{k}\cdot \nabla u_{\mathbf{k}} - e^{i\mathbf{k}\cdot \mathbf{r}} \nabla^2u_{\mathbf{k}}\right) \\ &= -\frac{1}{2}e^{i\mathbf{k}\cdot\mathbf{r}}(\nabla + i\mathbf{k})^{2}u_{\mathbf{k}} \end{align*} which shows that instead of solving (\ref{periodic_schrod}) we can instead solve \begin{align}\label{schrod_unit_cell} \mathcal{H}_{\mathbf{k}}u_{\mathbf{k}} := -\frac{1}{2}(\nabla + i\mathbf{k})^2u_{\mathbf{k}} + Vu_{\mathbf{k}} = E_{\mathbf{k}}u_{\mathbf{k}} \end{align} subject to periodic boundary conditions \(u_{\mathbf{k}}(\mathbf{r} + \mathfrak{t}) = u_{\mathbf{k}}(\mathbf{r})\).

This is a huge win. Without Bloch's theorem, we need to solve a Sturm-Liouville problem over a domain the size of \({\sim}10^{23}\) unit cells. With Bloch's theorem, we only need to sample some values of \(\mathbf{k}\) in the Brillouin zone and solve (\ref{schrod_unit_cell}) in one unit cell. But we can do even better than this. Since most crystals have point group symmetries in addition to translational symmetries, the number of \(\mathbf{k}\) values for which we need to solve (\ref{schrod_unit_cell}) is reduced dramatically. In addition, the point group symmetries will give some justification to our traditional presentation of band structure-requiring in general a 4D map-but by exploitation of symmetry can be presentated as a multivalued 2D graph of \(E\) vs \(\mathbf{k}(t)\) where \(t \mapsto \mathbf{k}(t)\) is known as a k-path.

Buckle up kids because this is gonna be hard. To understand these ideas we will need to understand

  • point groups and space groups
  • Seitz symbols
  • irreducible representations of some space groups
  • the orbit and stabilizer of the crystal momentum \(\mathbf{k}\)
  • high-symmetry points of the Brillouin zone
Let's get goin':

Define the orthogonal group in three dimensions by \begin{align*} \mathbf{O}(3) := \{ \mathbf{A} \in \mathbf{GL}(3, \mathbb{R}) \colon \mathbf{A}^{T} = \mathbf{A}^{-1}\} \end{align*} and the translation group \(\mathbf{T}(3) := \{\mathfrak{t} \in \mathbb{R}^3\}\) with the group operation being vector addition. We can combine \(\mathbf{O}(3)\) and \(\mathbf{T}(3)\) into the Euclidean group \(\mathbf{E}(3)\), which acts on \(\mathbb{R}^3\) by \begin{align*} \{\mathbf{A}|\mathfrak{t}\} \mathbf{r} := \mathbf{A}\mathbf{r} + \mathfrak{t} \end{align*} The notation \(\{\mathbf{A}|\mathfrak{t}\}\) is called the Seitz symbol. The multiplication of Seitz symbols is given by \begin{align*} \{\mathbf{A}_1|\mathfrak{t}_1\}\{\mathbf{A}_2|\mathfrak{t}_2\} = \{\mathbf{A}_1\mathbf{A}_2|\mathbf{A}_1\mathfrak{t}_2 + \mathfrak{t}_1\} \end{align*} Proving that \(\mathbf{E}(3)\) is a group is an easy calculation: The identity of \(\mathbf{E}(3)\) is \(\{\mathbf{I}|0\}\), \(\{\mathbf{A}|\mathfrak{t}\}^{-1} = \{\mathbf{A}^{-1}| -\mathbf{A}^{-1}\mathfrak{t}\}\), closure exhibited in the definition of the multiplication and associativity inherited. \(\{\mathbf{O}(3)|\mathbf{0}\}\) is a subgroup of \(\mathbf{E}(3)\), and \(\{\mathbf{I}|\mathbf{T}(3)\}\) is a normal subgroup: \begin{align*} \{\mathbf{A}|\tau\}\{\mathbf{I}|\mathfrak{t}\}\{\mathbf{A}|\tau\}^{-1} &= \{\mathbf{A}|\mathbf{A}\mathfrak{t} + \tau\}\{\mathbf{A}^{-1}|-\mathbf{A}^{-1}\tau\} = \{\mathbf{I}|\mathbf{A}\mathfrak{t}\} \end{align*} Because

  • \(\mathbf{O}(3)\) is a subgroup
  • \(\mathbf{T}(3)\) is a normal subgroup
  • \(\mathbf{O}(3)\cap \mathbf{T}(3) = \{\mathbf{I}| \mathbf{0}\}\)
  • and \(\mathbf{E}(3) = \{\mathbf{I}|\mathbf{T}(3)\}\{\mathbf{O}(3)|\mathbf{0}\}\)
we say that the Euclidean group is the semidirect product of \(\mathbf{O}(3)\) and \(\mathbf{T}(3)\) and write \(\mathbf{E}(3) = \mathbf{O}(3)\ltimes \mathbf{T}(3)\). (It is certainly not the direct product of \(\mathbf{O}(3)\) and \(\mathbf{T}(3)\) because the direct product group operation is \((\mathbf{A}_1, \mathfrak{t}_1)(\mathbf{A}_2, \mathfrak{t}_2) = (\mathbf{A}_1\mathbf{A}_2, \mathfrak{t}_1+\mathbf{t}_2)\).)

The Euclidean group is a continuous group. Insertion of a crystal lattice into \(\mathbb{R}^3\) breaks the continuous symmetry and gives us a discrete symmetry group as a subgroup of the Euclidean group. The set of all \(\{\mathbf{A}|\mathfrak{t}\} \in \mathbf{E}(3)\) that are also symmetries of the crystal gives us the space group \(\mathbf{S}\). Let \(\{\mathbf{A}|\mathfrak{t}\} \in \mathbf{S}\) be associated with a symmetry operator \(P_{\{\mathbf{A}|\mathfrak{t}\}}\) in the standard way by defining \(P_{\{\mathbf{A}|\mathfrak{t}\}}\psi(\mathbf{r}) := \psi(\{\mathbf{A}|\mathfrak{t}\}^{-1}\mathbf{r})\). Then the condition that \(\{\mathbf{A}|\mathfrak{t}\}\) is a symmetry of the crystal implies that \([\mathcal{H}, P_{\{\mathbf{A}|\mathfrak{t}\}}] = 0\) and \begin{align*} \mathcal{H} \psi = E\psi \iff \mathcal{H}P_{\{\mathbf{A}|\mathfrak{t}\}}\psi = EP_{\{\mathbf{A}|\mathfrak{t}\}}\psi \end{align*} The most general way to proceed is to study the representations of the space group directly. But this group is massive, containing \({\sim}10^{23}\) elements, and the representations become unwieldy. So instead it would be nice to extend our existing understanding of the translational symmetry of the crystal by studying the action of the rotations on the Bloch waves. If we could write \(\{\mathbf{A}|\mathfrak{t}\} = \{\mathbf{I}|\mathfrak{t}\}\{\mathbf{A}|\mathbf{0}\}\), this goal could be realized, but in general, this factorization cannot be achieved due to non-symmorphic operations. We must formalize to proceed further.

Let \(\mathbf{S}\) be the space group of some crystal and let \(\phi\colon \mathbf{S} \to \mathrm{Im}(\phi) \subset \mathbf{O}(3)\) be defined by \begin{align*} \phi(\{\mathbf{A}|\mathfrak{t}\}) = \mathbf{A} \end{align*} We call the image of \(\phi\) the point group of the crystal, which we will denote \(\mathbf{P}\). To clarify this definition, it might be of use to the reader to prove that \(\phi\) is a homomorphism and that \(\mathbf{P} =\mathrm{Im}(\phi)\) is indeed a group. Note that \(\{\mathbf{P}|\mathbf{0}\}\) is not necessarily a subset of \(\mathbf{S}\). There can exist non-symmorphic operations such that \(\{\mathbf{A}|\tau\} \in \mathbf{S}\), but \(\{\mathbf{A}|\mathbf{0}\} \not \in \mathbf{S}\) and \(\{\mathbf{I}|\tau\}\not\in \mathbf{S}\). Examples of non-symmorphic space group elements are screw axes and glide planes. Because the natural embedding of \(\mathbf{P}\) into \(\mathbf{O}(3)\) via \(\mathbf{A}\mapsto \{\mathbf{A}|\mathbf{0}\}\) is not necessarily even a subset of the space group \(\mathbf{S}\), the space group is no longer the semidirect product of the point group and translation group. The insertion of a lattice into empty space breaks not only the continuous symmetry of \(\mathbf{E}(3)\), but also its structure as a semidirect product.

The kernel of \(\phi\) are the translations \(\{\mathbf{I}|\mathfrak{t}\}\) where \(\mathfrak{t} = n_1\mathbf{a}_1 + n_2\mathbf{a}_2 + n_3\mathbf{a}_3\). (That translations are the only elements of the kernel is harder to prove, and I know of no slick way of demonstrating it. In the case that \(\{\mathbf{A}|\tau\}\) is a screw axis, the proof that \(\mathbf{A}\mathfrak{t} \in \mathbf{T} \, \forall \mathfrak{t} \in \mathbf{T}\) is shown graphically in Space Groups for Solid State Scientists. Please hit me up if you know a good proof.) Even though the discrete lattice has broken the semidirect product stucture, since translations are still a normal subgroup of the space group, we can form the factor group \(\mathbf{S}\setminus \mathbf{T}\). The point group \(\mathbf{P}\) is isomorphic to the factor group \(\mathbf{S}\setminus\mathbf{T}\). A representative of the factor group has the form \(\{\mathbf{A}|\tau\}\{\mathbf{I}|\mathbf{t}_{n}\}\), and under the image of \(\phi\) we obtain \(\mathbf{A}\), showing \(\phi\) covers \(\mathbf{P}\). Now assume that \(\phi(\{\mathbf{A}|\tau_1\}\{\mathbf{I}|\mathbf{t}_{n}\}) = \phi(\{\mathbf{A}|\tau_2\}\{\mathbf{I}|\mathbf{t}_{m}\})\). Because \(\phi\) is a homomorphism with translational kernel, \(\phi(\{\mathbf{A}|\tau_1\}) = \phi(\{\mathbf{A}|\tau_2\})\). But both \(\{\mathbf{A}|\tau_1\}\) and \(\{\mathbf{A}|\tau_1\}\) are elements of \(\mathbf{S}\), which is a group, so \begin{align*} \{\mathbf{A}|\tau_1\}\{\mathbf{A}|\tau_2\}^{-1} = \{\mathbf{I}|\tau_1 -\tau_2\} \in \mathbf{S} \end{align*} and hence \(\tau_1 - \tau_2 \in \mathbf{T}\) showing that \(\{\mathbf{A}|\tau_1\}\) and \(\{\mathbf{A}|\tau_2\}\) are elements of the same coset. Hence \(\phi\) is an isomorphism from \(\mathbf{S}\setminus \mathbf{T}\) to \(\mathbf{P}\).

The purpose of studying the isomorphism between the point group and the factor group \(\mathbf{S}\setminus \mathbf{T}\) becomes more clear in reciprocal space. The action of a space group operation on a Bloch function is \begin{align*} P_{\{\mathbf{A}|\tau\}}e^{i\mathbf{k}\cdot \mathbf{r}}u_{\mathbf{k}} = e^{i\mathbf{A}\mathbf{k}\cdot\tau}e^{i\mathbf{A}\mathbf{k}\cdot \mathbf{r}}P_{\{\mathbf{A}|\tau\}}u_{\mathbf{k}} \sim e^{i\mathbf{A}\mathbf{k}\cdot \mathbf{r}}P_{\{\mathbf{A}|\tau\}}u_{\mathbf{k}} \end{align*} where the \(\sim\) indicates the fact that wavefunctions are only defined up to a global phase. \(P_{\{\mathbf{A}|\tau\}}u_{\mathbf{k}}\) is invariant under lattice translations just like \(u_{\mathbf{k}}\), so \(P_{\{\mathbf{A}|\tau\}}e^{i\mathbf{k}\cdot \mathbf{r}}u_{\mathbf{k}}\) is also a Bloch wave. Because the effect of the nonsymmorphic translation \(\tau\) is a global phase shift, we see that the point group is the natural group to study in reciprocal space.

So clarify this idea, let \(\{\mathbf{A}|\tau\} \in \mathbf{S}\). Then \(\mathcal{H}_{\mathbf{A}\mathbf{k}}\) and \(\mathcal{H}_{\mathbf{k}}\) are isospectral. The idea is that \(\mathcal{H}_{\mathbf{A}\mathbf{k}}\) and \(P_{\{\mathbf{A}|\tau\}}\) almost commute, in that \begin{align}\label{almost_commute} \mathcal{H}_{\mathbf{A}\mathbf{k}}P_{\{\mathbf{A}|\tau\}} = P_{\{\mathbf{A}|\tau\}}\mathcal{H}_{\mathbf{k}} \end{align} Using standard manipulations \([\nabla^2, P_{\{\mathbf{A}|\tau\}}] = 0\), \([V, P_{\{\mathbf{A}|\tau\}}] = 0\), and \(\mathbf{A}\mathbf{k}\cdot \mathbf{A}\mathbf{k} = \mathbf{k}\cdot\mathbf{k}\) reduces the claim of (\ref{almost_commute}) to \begin{align*} \mathbf{A}\mathbf{k}\cdot \nabla (P_{\{\mathbf{A}|\tau\}}\phi) = P_{\{\mathbf{A}|\tau\}} \mathbf{k}\cdot\nabla\phi \end{align*} This says that the derivative of the rotated function in the direction of the rotated vector is the same as rotating the directional derivative, which is easily visualized. I find this mental picture of rotating the graph of a function fairly intuitive, but the same thing can be proved using the chain rule for the gradient. The claim is now easy: \begin{align*} \mathcal{H}_{\mathbf{k}}u_{\mathbf{k}} &= Eu_{\mathbf{k}} &\implies &\mathcal{H}_{\mathbf{A}\mathbf{k}}P_{\{\mathbf{A}|\tau\}}u_{\mathbf{k}} = P_{\{\mathbf{A}|\tau\}}\mathcal{H}_{\mathbf{k}}u_{\mathbf{k}} = EP_{\{\mathbf{A}|\tau\}}u_{\mathbf{k}} \\ \mathcal{H}_{\mathbf{A}\mathbf{k}}u_{\mathbf{A}\mathbf{k}} &= Eu_{\mathbf{A}\mathbf{k}} &\implies& \mathcal{H}_{\mathbf{k}}P_{\{\mathbf{A}|0\}}^{-1}u_{\mathbf{A}\mathbf{k}} = P_{\{\mathbf{A}|\tau\}}^{-1}\mathcal{H}_{\mathbf{A}\mathbf{k}}u_{\mathbf{A}\mathbf{k}} = EP_{\{\mathbf{A}|\tau\}}^{-1}u_{\mathbf{A}\mathbf{k}} \end{align*} Since \(\mathcal{H}_{\mathbf{k}}\) and \(\mathcal{H}_{\mathbf{A}\mathbf{k}}\) are isospectral, so are \(\mathcal{H}_{\mathbf{A}\mathbf{k}}\) and \(\mathcal{H}_{\mathbf{A}^2\mathbf{k}}\). This suggests that we study the orbit of \(\mathbf{k}\) under the action of the point group \(\mathbf{P} \cong \mathbf{S}\setminus \mathbf{T}\): \begin{align*} \mathbf{P}\cdot \mathbf{k} := \{\mathbf{A}\mathbf{k} \colon \mathbf{A} \in \mathbf{P}\} \end{align*} The orbits of the crystal momenta partition the Brillouin zone into equivalence classes, because the relation "\(\mathbf{k}_1\) is in the orbit of \(\mathbf{k}_2\)" is an equivalence relation. Hence we only need to solve (\ref{schrod_unit_cell}) once for each equivalence class, because all other representatives of the class have the same energy levels. The smallest convex subset of the Brillouin zone that contains a representative from each orbit is called an irreducible Brillouin zone. (Note: If the crystal lacks inversion symmetry, then \(E(\mathbf{k})\) and \(E(-\mathbf{k})\) must be calculated independently, and the irreducible Brillouin zone consists of two convex components which meet at \(\mathbf{k} = 0\).)

Although we \(\mathcal{H}_{\mathbf{A}\mathbf{k}}\) and \(\mathcal{H}_{\mathbf{k}}\) are isospectral, we do not say that \(u_{\mathbf{A}\mathbf{k}}\) and \(u_{\mathbf{k}}\) are degenerate--since we reserve that term for linearly independent wavefunctions solving the same Hamiltonian. There are circumstances where \(u_{\mathbf{A}\mathbf{k}}\) and \(u_{\mathbf{k}}\) are in fact degenerate: When \(\mathbf{A}\cdot \mathbf{k} = \mathbf{k} + \mathbf{K}\) where \(\mathbf{K}\) is a reciprocal lattice vector. This suggests we study the stabilizer group of \(\mathbf{k}\): \begin{align*} \mathrm{Stab}(\mathbf{k}) := \{ \mathbf{A} \in \mathbf{P}: \mathbf{A}\mathbf{k} = \mathbf{k} + \mathbf{K}\} \end{align*} It is a standard result that the stabilizer is a subgroup of the point group. If \(|\mathrm{Stab}(\mathbf{k})| > 1\), then we say that \(\mathbf{k}\) is a high-symmetry point.

High-symmetry points are used to create band structure diagrams. The band structure the information contained in the multivalued function \(E(\mathbf{k})\). Since \(E(\mathbf{k} + \mathbf{K}) = E(\mathbf{k})\) for any \(\mathbf{K}\) in the reciprocal lattice, then we understand the band structure if we understand it in a single Brillouin zone. In addition, since \(E(\mathbf{A}\mathbf{k}) = E(\mathbf{k})\), we understand the bandstructure if we understand it in the irreducible Brillouin zone (\(\mathrm{IBZ}\)). Unfortunately, even with this simplification, \begin{align*} E\colon \mathrm{IBZ}\subset \mathbb{R}^{3} \to \mathbb{R} \end{align*} so we need 4-dimensional graph paper to plot our bands. The way that physicists have wriggled out of this dilemma is to create k-paths, which are maps \(t \mapsto \mathbf{k}(t)\). Then the multi-valued 2D graph of \(E(\mathbf{k}(t))\) vs \(t\) is displayed. The problem with this is that there is no obvious choice for the k-path of a given \(\mathrm{IBZ}\). Band structure diagrams computed by use k-paths defined by Hinuma, et. al, but other reasonable choices can be made. One guiding principle of choosing a k-path is that we should connect high-symmetry points by traveling along high-symmetry lines. Since high-symmetry points exhibit more degeneracy than generic points in the \(\mathrm{IBZ}\), there are fewer bands along these lines, making the task of interpretation of the band structure diagram a little easier. There is no guarantee that all the information of interest is contained along a chosen k-path, and war stories related to confusion from this problem abound. (Note: Hopefully this problem will be solved soon: I'm working on the code to present surfaces above slices through the IBZ. The user will be able to evolve these surfaces by manipulation of the pitch and roll of the slice. Sign up for an account to be notified when this feature drops.)

To illustrate these ideas, we'll consider graphene. Graphene is a great example, because it's a 2D material with a symmorphic space group. In common three-dimensional materials, there are so many group operations that we get worn out long before we can extract some insight from the problem. The point group of graphene is \(D_{6h} := D_{6}\times \{E, \sigma_{h}\}\), where \(\{E, \sigma_{h}\}\) is the reflection in the plane containing the carbon atoms. This group is incredibly simple so we can understand its representations immediately. By the dimension formula we know that all its finite-dimensional irreps are 1-dimensional, and there are two of them: The trivial representation \(\Gamma_0(E) = \Gamma_0(\sigma_h) = 1\) and \(\Gamma_1(\sigma_h) = -1\). Hence for any solution \(\psi\) of (\ref{periodic_schrod}) in graphene which transforms as an irreducible representation of \(\{E, \sigma_h\}\) \begin{align*} P_{\sigma_h}\psi(x,y,z) = \psi(x, y, -z) = \pm \psi(x,y,z) \end{align*} Hence if we restrict ourselves to consideration of solutions which are either even or odd under reflection in the \(z\)-plane, we can understand the representations of \(D_{6h}\) by understanding representations of \(D_{6}\). \(D_6\) is the symmetry group of the regular hexagon. The symmetry operations are the six rotations by multiples of \(2\pi/6\), the three reflections in the planes which bisect the faces of the hexagon, and the three reflections in the planes that intersect the vertices of the hexagon. The entire multiplication table of \(D_{6}\) can be generated by two elements \(a, b\) which obey \begin{align}\label{d6_mults} a^6 = e,\quad b^2 = e,\quad ba = a^{-1}b. \end{align} Using these multiplication rules we find 6 conjugacy classes \begin{align*} \{e\}, \{a^3\}, \{a, a^5\}, \{a^2, a^4\}, \{b, a^2b, a^4b\}, \{ab, a^3b,, a^5b\} \end{align*} The group is order 12, so by the dimension formula \begin{align*} 12 &= 3^2 + 3\cdot 1^2 \\ &= 2\cdot 2^2 + 4\cdot 1^2 \\ &= 2^2 + 8\cdot 1^2 \\ &= 12\cdot 1^2 \end{align*} we have the possible dimensions of the irreducible representations. But since we have 6 conjugacy classes (=#irreps), then we know that in fact the only acceptable solution is \(12 = 2\cdot 2^2 + 4\cdot 1^2\) (two 2d irreps and four 1d irreps). Hence the possible degeneracies of solutions to (\ref{periodic_schrod}) in graphene are either 1 or 2, with no other possibilities.

The Brillouin zone of graphene. The irreducible wedge is shaded light blue. High symmetry points on the interior of the IBZ are labelled with Greek letters, those on the faces are labelled with Roman letters. Red and purple lines indicate reflection planes, red lines are in one conjugacy class, purple in another. The \(\Gamma\) point has the entire point group as its stabilizer, the M point has the identity, vertical and horizontal reflection as its stabilizer. The blue lines are the orbit of a random \(\mathbf{k}\in \mathrm{IBZ}\), generated by six applications of \(a = \) rotation by \(2\pi/6\) and \(b = \) reflection over horizontal line.
  • Seitz, Frederick. “On the Reduction of Space Groups.” Annals of Mathematics, vol. 37, no. 1, 1936, pp. 17–28. JSTOR, JSTOR,
  • Glazer, Michael, Gerald Burns, and Alexander N. Glazer. Space groups for solid state scientists. Elsevier, 2012.
  • Hinuma, Yoyo, et al. "Band structure diagram paths based on crystallography." Computational Materials Science 128 (2017): 140-184.
  • Dresselhaus, Mildred S., Gene Dresselhaus, and Ado Jorio. Group theory: application to the physics of condensed matter. Springer Science & Business Media, 2007.
  • Tinkham, Michael. Group theory and quantum mechanics. Courier Corporation, 2003.