Written by spxtr. Please feel free to email me with any questions or comments. This document assumes familiarity with quantum mechanics.
Hofstadter's butterfly is one of the most aesthetically pleasing results in all of physics. It is an intricate fractal, yet at the same time it is the solution to a simple physical system. It carefully toes the line between pathological mathematics and experimental physics. Let's learn how that works.
We start with a single-particle tight-binding model on a square lattice with lattice constant \(a\) [1]. Let \(x=ma\) and \(y=na.\) The basis states are then \(|m,n\rangle\) for \(m,n\in\mathbb{N}.\) We next assume the simple Hamiltonian $$\hat{H}=\sum_{m,n}|m+1,n\rangle\langle m,n| +|m,n+1\rangle\langle m,n| + \mathrm{h.c.}.$$ We use the Landau gauge \(\vec{A}=Bx\hat{y}.\) The Peierl's subsitution, whereby the translation operators pick up a phase of $$\theta = \frac{q}{\hbar}\int\vec{A}\cdot\mathrm{d}\vec{r},$$ does not affect hopping terms in \(x.\) Hopping terms in \(y,\) however, pick up a phase of \(eBxa/\hbar.\) Let \(\alpha=\Phi/\Phi_0=eBa^2/h,\) the ratio of the magnetic flux per unit cell to the magnetic flux quantum [2]. With this addition, the Hamiltonian becomes $$\begin{equation}\hat{H}=\sum_{m,n}|m+1,n\rangle\langle m,n| + e^{2\pi i \alpha m}|m,n+1\rangle\langle m,n|+\mathrm{h.c.}.\label{eq:twodeehamiltonian}\end{equation}$$
None of the hopping terms in the Hamiltonian depend on \(n,\) so we take as an ansatz plane wave behavior in the \(y\) direction, so \(|m,n\rangle\rightarrow e^{-i\nu n}|m\rangle .\) After expanding the Hermitian conjugate, we find $$\hat{H}=\sum_m|m+1\rangle\langle m| + |m-1\rangle\langle m| + 2\cos\left(2\pi\alpha m + \nu\right)|m\rangle\langle m|.$$ Plug it into the time-independent Schrödinger equation \(\hat{H}\Psi=E\Psi,\) and we can write it as a simple difference equation $$\begin{equation}g_{m+1} + g_{m-1} + 2\cos(2\pi m\alpha + \nu)g_m=Eg_m,\label{eq:harper}\end{equation}$$ where the full wave function is $$\Psi(x,y)\propto \sum_{m,n} g_m e^{-i\nu n}\Phi(x-ma,y-na)$$ for a single atomic orbital \(\Phi.\) Equation \(\eqref{eq:harper}\) is known as Harper's equation or the Almost Mathieu operator, depending on who you ask. There are multiple ways of computing the spectrum of this Hamiltonian. We will start with Hofstadter's original method [3].
We can rewrite equation \(\eqref{eq:harper}\) as $$\begin{pmatrix} g_{m+1} \\ g_m \end{pmatrix} = \begin{pmatrix}E - 2\cos(2\pi m\alpha+\nu) & -1 \\ 1 & 0\end{pmatrix}\begin{pmatrix} g_m \\ g_{m-1} \end{pmatrix}.$$ We call this \(2\times 2\) matrix \(A(m).\) In order to be physically realistic, \(g\) must remain bounded for all \(m.\) If \(A\) is periodic in \(m\) with period \(q,\) then there must be an integer \(p\) such that $$2\pi\alpha(m+q)+\nu=2\pi\alpha m+\nu+2\pi p,$$ which implies \(\alpha = p/q,\) and \(\alpha\) is rational. We will proceed with this rationality as an assumption. We call the product of \(q\) successive \(A\) matrices \(Q,\) and it is a function of \(E\) and \(\nu.\) For \(g\) to remain bounded, the eigenvalues of \(Q\) must have a magnitude of at most 1. This will be the case if the trace of \(Q(E,\nu)\) is less than or equal to 2 in magnitude. Hofstadter then argues that we can eliminate the factor of \(\nu\) by restricting the trace to 4, yielding the final inequality $$|\mathrm{Tr}Q(E)|\leq 4.$$ This is now straightforward to solve: let \(\alpha=p/q,\) compute \(Q(E),\) and plot its trace. Where it is less than 4 in magnitude, those are the allowed energies.
The above plots show \(\mathrm{Tr}Q(E)\) vs \(E\) for \(\alpha=p/q=1/2\) (left), \(1/3\) (middle), and \(4/7\) (right). The dashed horizontal lines are \(\pm 4\) and the thick red horizontal lines mark where the trace is within the lines. Note that for \(\alpha=1/2,\) the trace is equal to \(-4\) at \(E=0,\) so the two subbands kiss in the middle. The trace of \(Q\) is a degree-\(q\) polynomial, and so we expect it to have \(q\) roots and thus approximately \(q\) regions where its magnitude is under 4. Put succinctly, when the magnetic flux is a rational number \(p/q,\) the Bloch band splits into \(q\) subbands.
This statement is both alarming and intriguing. Small variations in \(\alpha\) can lead to arbitrarily large changes in \(q,\) which seems at first glance to render the entire calculation unphysical. Fortunately, there is an important aspect of continuity to the spectra of nearby values of \(\alpha\) that saves the day.
As \(q\) grows, the width of the subbands become thinner and thinner. Irrational numbers can be (very) loosely thought of as rational numbers with infinitely large numerator and denominator, so we might guess that the spectrum for irrational \(\alpha\) will consist of infinitely-many subbands with 0 length, perhaps somehow equivalent to a Cantor set. This guess turns out to be true [4].
By choosing all \(q\leq q_{\mathrm{max}},\) \(p\) such that \(0\lt p\lt q\) with \(p\) and \(q\) coprime, and then computing the trace, we find the full spectrum, shown above for \(q_{\mathrm{max}}=20.\) The y axis ranges from \(\alpha=0\) to \(1,\) and the x axis ranges from \(E=-4\) to \(4.\) A few notes on symmetries:
Returning to equation \(\eqref{eq:twodeehamiltonian}\): if \(\alpha\) is rational, then the matrix elements of the Hamiltonian repeat after \(q\) steps. This motivates us to define a larger magnetic unit cell, the same size in \(y\) but \(q\) times larger in \(x,\) in which we can apply Bloch's theorem.
On the left, the unit cell is shaded and some hopping terms are shown for \(\alpha=0\). On the right, \(\alpha=1/3\), each hopping term in \(y\) picks up a phase of \(\exp(2\pi i/3)\) compared to its neighbor on the left, and so the magnetic unit cell is tripled in size in \(x\) but not in \(y.\) We then solve the problem in reciprocal space in the magnetic Brillouin zone \(0\leq k_x \lt 2\pi/aq,\) \(0\leq k_y \lt 2\pi/a.\)
To solve this easily, we set \(q\) to some large prime number so that \(k_x\approx 0,\) and then sweep \(p\) from 1 to \(q-1,\) solving the associated eigenvalue equation for each value of \(p/q.\) The Hamiltonian for each \(p/q\) is a \(q\times q\) matrix with the form $$\hat{H}=\begin{pmatrix}V_0 & 1 & 0 & 0 & \cdots & 0 & 1\\ 1 & V_1 & 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & V_2 & 1 & \cdots & 0 & 0\\ 0 & 0 & 1 & V_3 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & V_{q-2} & 1 \\ 1 & 0 & 0 & 0 & \cdots & 1 & V_{q-1} \end{pmatrix}$$ for \(V_m=2\cos(2\pi mp/q + k_y).\) This matrix can be numerically solved by a computer in \(\mathcal{O}(q^3)\) time [5], shown below for \(q=997\) (subbands for \(q=1/3\) and \(2/3\) highlighted in red).
At \(\alpha=p/q,\) the Bloch band splits into \(q\) subbands. Each of these subbands has the same integrated density of states as the others [6], so if the total number of states in the band is equal to 1, then each subband has \(1/q\) states. This is in spite of the fact that their widths are not identical.
Curiously, the gaps between states have a distinctly continuous appearance. If we follow these gaps along and keep track of the number of states with lower energy, then we recover behavior that does not appear to depend on the rationality of the magnetic field. Instead, the gaps between subbands follow straight lines in density/field space. This was first pointed out by Gregory Wannier in 1978 [7].
The plot on the left is the butterfly as computed above. On the right, I show the inverse density of states \(\mathrm{d}\mu/\mathrm{d}n\) as a function of normalized density (\(x\)-axis) and magnetic field \(y\)-axis. Dark lines indicate a large gap from one state to the next. The two highlighted lines correspond to the gaps indicated on the left. These trajectories can be written as a Diophantine equation $$n=t\alpha + s$$ for integers \(s\) and \(t.\) In traditional field-effect nanodevices, the coplanar gates allow for direct control of electron density. In those devices, we therefore expect to see signatures of this Wannier plot in transport measurements. It is a beautiful thing that such a complicated fractal spectrum yields such a simple physical description: each gap is labelled by only two integers.
The triangular lattice (left) is topologically distinct from the square lattice in that each lattice point has six neighbors, rather than four. The calculation of the spectrum in a magnetic field follows the same course. We end up with $$\begin{eqnarray*}\hat{H}&=&2\sum_m[\cos(2\pi \alpha m)|m\rangle\langle m| \\&+& \cos(\pi\alpha(m+1/2))|m+1\rangle\langle m| \\&+& \cos(\pi\alpha(m-1/2))|m-1\rangle\langle m|].\end{eqnarray*}$$
The honeycomb lattice (right) can be thought of as either two overlaid triangular sublattices (here one shown unshaded, the other shaded red), or as a single triangular lattice with two atoms in its unit cell. Each site is adjacent to three sites from the opposite sublattice. We find $$\begin{eqnarray*}\hat{H}&=&\sum_m[\exp(2\pi i \alpha m / 3)|m,B\rangle\langle m,A| \\&+& \exp(-\pi i \alpha (m + 1/2) / 3)|m+1,B\rangle \langle m,A|\\&+&\exp(-\pi i \alpha (m-1/2)/3)|m-1,B\rangle\langle m,A|] + \mathrm{h.c.},\end{eqnarray*}$$ where \(A\) and \(B\) index the sublattice. The Hamiltonian has now doubled in size, but otherwise it is easy to compute.
The solutions for the triangular lattice (left) and honeycomb lattice (right) have the same butterfly-like appearance, but the details of each are different. The triangular lattice looks like it has only one wing, and the spectrum is not symmetric about 0 energy or 1/2 field. The honeycomb lattice looks like two butterflies side-by-side. The Wannier plots for the two lattices again have Diophantine gaps, however the presence of specific gaps and the relative sizes of the gaps are distinct from the square lattice and each other.
The gaps in the butterfly spectra follow straight lines given by the diophantine equation \(n=t\alpha+s\) for integers \(t\) and \(s\) [7]. The number \(t\) is, somewhat surprisingly, proportional to the Hall conductivity of the system when the Fermi energy is tuned into the respective gap: $$\sigma_H = \frac{te^2}{h}.$$ This fact, along with an algorithm for computing the numbers \(t\) and \(s\) for a given gap at a given field, along with a very important in-between equation relating the Hall conductivity to a particular integral around the Brillouin zone, was derived in a now-famous paper by Thouless, Kohmoto, Nightingale, and den Nijs, hereafter referred to as the TKNN paper [8]. This work was a key component in the 2016 Nobel Prize in Physics (pdf).
The algorithm for computing \(s\) and \(t\) is simple. For \(\alpha=p/q,\) there are \(q\) subbands with \(q-1\) gaps between them. For the \(r\)-th gap, solve the diophantine equation $$r=s_rq+t_rp$$ for integer \(r\) with \(1 \leq r < q\) and \(|t_r| \leq q/2.\) For instance, at a field of \(p/q=2/5\) and \(r = 1,\) the equation \(1=5s+2t\) has solution \(s=1\) and \(t=-2.\) So, the Hall conductivity within that gap is \(-2e^2/h.\) This is pure magic.
On the left I show the butterfly, with the Hall resistivity (\(1/t_r\)) for each gap indicated with color, ranging from \(-h/e^2\) (blue) to \(h/e^2\) (red). On the right, the associated Wannier plot.
If that's not a nice result then I don't know what is.