## 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*.

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.

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.

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. https://doi.org/10.1007/BF03012263
- 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).