4 Tight binding
A third perspective, standing at the opposite limit of the NFE model, starts from localized atomic orbitals. This is called tight binding (TB) or, in a chemical context, linear combination of atomic orbitals (LCAO).
4.1 Basic example
Many concepts have to be introduced, so we start from a minimalistic example to grasp them one by one. A basis built out of atomic orbitals has a major problem: orbitals are localized and using them as a basis does not take advantage of the great simplification of Bloch’s theorem… we thus first of all need a basis change.
Definition 4.1: Bloch sums (simplified)
Let us assume we have a 1D crystal with \(N\) cells and lattice spacing \(a\), in every cell we have a single atom with a single orbital \(\phi_a(x)\). The Bloch sum of quasi-momentum \(k\) is defined as
\[\psi_k(x) = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} e^{ikt_n}\, \phi_a(x-t_n), \tag{4.1}\]
where \(t_n=na\) and the sum runs over all \(N\) lattice sites of the crystal. As already noticed in previous chapters, in a finite system \(k\) values are discrete, with \(k=2\pi m/Na\) and \(m\) running from \(0\) to \(N-1\).
Exercise. Show that the Bloch sum \(\psi_k\) is a Bloch state with quasi-momentum \(k\).
CautionIs this a good new basis?
Let us run a few basic checks:
Dimension. There are \(N\) inequivalent \(k\) choices so this check is passed: the basis still contains \(N\) elements.
Orthogonality. This will get more complicated with more orbitals, but here it is easy to conclude that the new basis is orthogonal. Each element is a Bloch state with an inequivalent \(k\), and Bloch states with inequivalent \(k\) are orthogonal by definition: they are not mixed by any periodic Hamiltonian, including the identity \(\mathcal{H}=\mathbb{1}\).
Normalization. This is less obvious, and a proper calculation should be performed. The bra-ket
\[\langle \psi_k|\psi_k\rangle = \frac{1}{N}\sum_{n,n'} e^{ik(na-n'a)}\,\langle\phi_a(x-n'a)|\phi_a(x-na)\rangle,\]
can be greatly simplified if we express it in terms of \(m=n-n'\) and notice how it contains \(N\) identical copies
\[ \begin{aligned} \langle \psi_k|\psi_k\rangle &= \frac{1}{N}\sum_{n',m} e^{ikma}\,\langle\phi_a(x-n'a)|\phi_a(x-n'a-ma)\rangle\\ &= \sum_m e^{ikma}\,\langle\phi_a(x)|\phi_a(x-ma)\rangle \end{aligned} \]
In the TB spirit, we expect orbitals \(\phi_a(x)\) to be localized. While orbitals on different lattice sites are not orthogonal in general, their overlap is vanishingly small and it is reasonable to assume \(\langle\phi_a(x)|\phi_a(x-ma)\rangle\approx \delta_{m0}\). This finally leads to \(\langle \psi_k|\psi_k\rangle=1\), thus this is true in the limit of zero overlap between orbitals in nearby sites.
Anticipating future issues. Note that an inexact normalization (and, even worse, an inexact orthogonality when we have multiple orbitals or basis atoms) of the basis does not necessarily mean Bloch sums cannot be used to solve the problem, but additional care must be taken in the diagonalization problem.
Before moving to this ultra-simplified TB diagonalization, let us recap the hypotheses:
- the crystal is 1D, with lattice spacing \(a\) and one atom per cell (trivial basis);
- a single orbital \(\phi_a(x)\) per atom is included in the calculation;
- we have \(N\) cells with periodic boundary conditions;
- the Hamiltonian only enables hopping between the first neighboring atoms.
Given the localized basis \(|\phi_a(x-t_n)\rangle\) with \(t_n=na\), the last condition is encoded in the hopping integrals \(\langle\phi_a(x)|\mathcal{H}|\phi_a(x-ma)\rangle\) measuring the coupling between orbitals at sites separated by \(m\) lattice spacings. Here we only keep the nearest neighbor terms:
- \(\langle\phi_a(x)|\mathcal{H}|\phi_a(x)\rangle=\varepsilon_0\): the on-site energy;
- \(\langle\phi_a(x)|\mathcal{H}|\phi_a(x\pm a)\rangle = \gamma\) for nearest neighbors (\(\gamma < 0\));
- zero hopping otherwise.
In this very special limit exercise, the problem is so trivial that \(\mathcal{H}\) is already diagonal over the Bloch sums (basically, having restricted the problem to a single atom and single orbital, the \(S_k\) subspace is \(1\times 1\)). The band dispersion can be calculated as the expectation value \(E(k) = \langle\psi_k|\mathcal{H}|\psi_k\rangle\). Note we are assuming \(\langle\psi_k|\psi_k\rangle = 1\), otherwise a renormalization factor would be needed:
\[E(k) = \frac{1}{N}\sum_{n,n'} e^{ik(na-n'a)}\,\langle\phi_a(x-n'a)|\mathcal{H}|\phi_a(x-na)\rangle,\]
Substituting again \(m = n-n'\) and summing over \(n'\) (which gives a factor \(N\)) we obtain the dispersion
\[ \begin{aligned} E(k) &= \sum_{m=-\infty}^{+\infty} e^{ikma}\langle\phi_a(x)|\mathcal{H}|\phi_a(x-ma)\rangle\\ &= \sum_{m = 0,\pm 1} e^{ikma}\langle\phi_a(x)|\mathcal{H}|\phi_a(x-ma)\rangle\\ &= \varepsilon_0 + \gamma e^{-ika} + \gamma e^{+ika} = \varepsilon_0 + 2\gamma\cos(ka). \end{aligned} \tag{4.2}\]
This dispersion is displayed in Figure 4.1.
CautionBenzene as a TB textbook example
The molecule of benzene (C\(_6\)H\(_6\)) offers a textbook implementation of the physics discussed so far, where periodic boundary conditions are not a mathematical abstraction but rather a physical reality. Benzene has a skeleton consituted by six carbon atoms and hydrogen terminations arranged in an hexagonal network of \(sp^2\) bonds. This structure saturates most of available electrons except for six last ones residing in the carbon \(p_z\) orbitals. Organic chemistry typically describes this in terms of resonance hybrid between two competing \(\pi\) bond configurations. Our simplified TB allows a much more profound and accurate description, which is known in chemistry as HĂĽckel model.
In the language of this section, each carbon atom has a localized \(p_z\) orbital with some hopping probability towards the nearby atoms: this reproduces well the situation in Section 4.1 with \(N = 6\) and, given this is an actual ring, true periodic boundary conditions. In the limit of nearest neighbor hopping and neglecting the other electrons in the \(sp^2\) bonds, the Bloch sums
\[\psi_k = \frac{1}{\sqrt{6}}\sum_{n=0}^{5} e^{ikt_n}\,\phi_{p_z}(\mathbf{r}-\mathbf{t}_n), \qquad ka = \frac{2\pi m}{6}, \quad m = 0,\pm 1,\pm 2, 3,\]
diagonalize the problem exactly. The six molecular orbitals are Bloch states running on the benzene ring, with an energy sitting on the cosine band \(E(k) = \varepsilon_0 + 2\gamma\cos(ka)\). The actual orbital energy derive from \(k\) quantization: a single level at \(\varepsilon_0+2\gamma\), a degenerate pair at \(\varepsilon_0+\gamma\), another pair at \(\varepsilon_0-\gamma\) and a single level at \(\varepsilon_0-2\gamma\) (recall \(\gamma<0\)). The six \(\pi\) electrons fill the three lowest states, the closed-shell “aromatic sextet” behind the exceptional stability of the molecule.
Orbitals can be visually explored in Figure 4.2. Note that real-valued wavefunctions are typically preferred since they are easier to manage: consistenly, Bloch states are combined into even and odd superpositions \(\psi_{+k}\pm\psi_{-k}\), this is completely equivalent to using \(p_x\) and \(p_y\) orbitals instead of the eigenstates of the angular momentum with \(m=\pm 1\).
4.2 Semi-empirical tight binding
Let us now add one missing ingredient with respect to Section 4.1 — the expansion over multiple orbitals — while keeping things simple with a trivial basis (one atom per cell). The extra complication of a genuine multi-atom basis is postponed to the graphene case study of Section 4.3, where it is really needed.
Before diving into the formalism, note here we focus on semi-empirical TB meaning we will make relatively drastic simplifications with the explicit goal of reducing the band problem to an analytic model containing only a handful of parameters. These are then not computed from first principles but fitted to data. This is at the opposite end of the spectrum with respect to ab initio methods, where one computes everything from first principles, with little/no feedback from experimental data.
Definition 4.2: Bloch sums (multiple orbitals)
For a crystal with one atom per cell carrying several atomic orbitals \(\phi_i\) (e.g. \(s, p_x, p_y, p_z, \dots\)), one Bloch sum is constructed per orbital:
\[\psi_{i}(\mathbf{k},\mathbf{r}) = \frac{1}{\sqrt{N}}\sum_{\mathbf{R}\in\mathcal{BL}} e^{i\mathbf{k}\cdot\mathbf{R}}\,\phi_i(\mathbf{r}-\mathbf{R}). \tag{4.3}\]
This is the same construction as in Section 4.1, now repeated once per orbital rather than for a single one.
Definition 4.3: Tight binding secular equation
The key idea is that crystal eigenstates at fixed \(\mathbf{k}\) are linear combinations \(\psi = \sum_{i}c_{i}\psi_{i}\). After projection the classic eigenvalue equation \(\mathcal{H}\psi = E\psi\) will lead to a generalized eigenvalue problem
\[\sum_{j}\bigl[\mathcal{M}_{ij}(\mathbf{k}) - E\,\mathcal{S}_{ij}(\mathbf{k})\bigr]\,c_{j} = 0, \tag{4.4}\]
with an interaction matrix \(\mathcal{M}_{ij} = \langle\psi_{i}|\mathcal{H}|\psi_{j}\rangle\) and overlap matrix \(\mathcal{S}_{ij} = \langle\psi_{i}|\psi_{j}\rangle\).
The band energies are the roots of \(\det[\mathcal{M}(\mathbf{k}) - E\,\mathcal{S}(\mathbf{k})] = 0\).
CautionWhat if \(\mathcal{S}\) is different from identity?!
This is why the problem is called a “generalized” eigenvalue problem. The common and easy-to-grasp approach (Löwdin method) consists in going back to an orthogonal basis by noticing that the overlap matrix \(\mathcal{S}\), being Hermitian and positive-definite, can be diagonalized and has a well-defined Hermitian square root \(\sqrt{\mathcal{S}}\); let us call its inverse \(\mathcal{Q}\), so that \(\mathcal{S}=\mathcal{Q}^{-2}\). By multiplying the secular equation from the left by \(\mathcal{Q}\) and inserting the identity \(\mathcal{Q}\mathcal{Q}^{-1}\), the eigenvalue problem becomes a standard one:
\[\sum_{j}\bigl[\tilde{\mathcal{M}}_{ij}(\mathbf{k}) - E\,\delta_{ij}\bigr]\,\tilde{c}_{j} = 0,\]
where we have a new Hermitian interaction matrix \(\tilde{\mathcal{M}} = \mathcal{Q}^\dagger\mathcal{M}\mathcal{Q}\) and transformed eigenvectors \(\tilde{c} = \mathcal{Q}^{-1}c\).
This is not the only possibility, and in addition the subject would deserve a more profound discussion on numerical stability (e.g., using Cholesky factorization in actual computational implementations).
4.2.1 Key approximations
We now approximate Equation 4.4 to make it tractable without losing the essential physics.
1. Neglect of orbital overlap. We assume \(\langle\phi_i(\mathbf{r})|\phi_j(\mathbf{r}-\mathbf{t})\rangle \approx 0\) for \(\mathbf{t}\neq\mathbf{0}\), so \(\mathcal{S} = \mathbb{1}\) and the problem reduces to ordinary diagonalization of \(\mathcal{M}(\mathbf{k})\).
2. Decomposition of \(V(\mathbf{r})\) and the two-center approximation. Here we decompose
\[V(\mathbf{r}) = V_a(\mathbf{r}) + V'(\mathbf{r})\]
with \(V_a\) as atomic potential in the cell located at the origin, while \(V'(\mathbf{r})\) accounts for all the rest. Consider now a generic matrix element between an orbital at the origin and an orbital at site \(\mathbf{t}\) (a lattice vector, since the basis is trivial). Since atomic orbitals solve the atomic problem, \([\mathbf{p}^2/2m + V_a(\mathbf{r})]\,\phi_i(\mathbf{r}) = E_i\,\phi_i(\mathbf{r})\), easy algebra leads to
\[\langle\phi_i(\mathbf{r})|\mathcal{H}|\phi_j(\mathbf{r}-\mathbf{t})\rangle = E_i\,\langle\phi_i(\mathbf{r})|\phi_j(\mathbf{r}-\mathbf{t})\rangle + \langle\phi_i(\mathbf{r})|V'(\mathbf{r})|\phi_j(\mathbf{r}-\mathbf{t})\rangle. \tag{4.5}\]
The first term survives on-site only (\(E_i\,\delta_{ij}\) at \(\mathbf{t}=\mathbf{0}\)): same-site orbitals are orthonormal, while overlaps at \(\mathbf{t}\neq\mathbf{0}\) vanish by the previous approximation. In the second term, \(V'\) is itself a sum of atomic potentials sitting on all sites but the origin, and each contribution can be classified by its number of distinct spatial centers, as illustrated in Figure 4.3. Assuming here \(\mathbf{t}\neq\mathbf{0}\), we can only retain the potential centered on the site of the second orbital — the two-center integral \(\langle\phi_i(\mathbf{r})|V_a(\mathbf{r}-\mathbf{t})|\phi_j(\mathbf{r}-\mathbf{t})\rangle\) — and neglect all three-center integrals (orbitals on two distinct sites, potential on a third one): these are expected to be smaller since at the each center we have a the product of two tails.
3. Diagonal crystal field. Let us moved to the skipped case \(\mathbf{t}=\mathbf{0}\), which can yield a non-zero term in the second \(V'(\mathbf{r})\) term of Equation 4.5. This evaluates ho much the all the other atoms in the crystal - this is called crystal field - perturb the orbitals located at the origin \(I_{ij} = \langle\phi_i|V'|\phi_j\rangle\) (right panel of Figure 4.3). The term can couple different orbitals or lift degeneracies, but here it is assumed diagonal, \(I_{ij} \approx I_i\,\delta_{ij}\): this merely renormalises the on-site energies \(E_i \to E_i + I_i\) and neglects off-diagonal symmetry-breaking contributions.
4. Nearest-neighbor hopping. The lattice sum is restricted to first neighbors; further shells can be added for higher accuracy.
Under all four approximations, the matrix elements take the compact form
\[\mathcal{M}_{ij}(\mathbf{k}) = E_i\,\delta_{ij} + \sum_{\mathbf{t}_I\in\mathcal{N\!N}}\langle\phi_i(\mathbf{r})|V_a(\mathbf{r}-\mathbf{t}_I)|\phi_j(\mathbf{r}-\mathbf{t}_I)\rangle\,e^{i\mathbf{k}\cdot\mathbf{t}_I},\]
with the on-site energies \(E_i\) now including the crystal-field shifts \(I_i\), and the sum running over the nearest-neighbor lattice vectors \(\mathbf{t}_I\) only.
CautionSlater–Koster two-center integrals
The aim of empirical TB is to produce an analytical model with a small number of parameters, to be then determined by data fitting. Given crystals often have high coordination numbers, a simplification comes from the rotational properties of atomic orbitals. Here, we assume a spherically symmetric atomic potential \(V_a(r)\) and orbitals expressed in terms of spherical harmonics and \(s\), \(p\), etc. orbitals. \(\mathcal{M}\) will typically contain multiple two-center integrals at fixed nearest neighbor distance \(R\), involving the same orbitals, but with different bond orientations \(\hat{\mathbf{n}} = (n_x,n_y,n_z)\), where \(\mathbf{R}=R\hat{\mathbf{n}}\).
Slater–Koster parameters. Consider a two-center integral involving \(s\) and \(p_x\) orbitals: given \(p_x\) is odd upon \(x\) mirror, any bond oriented perpendicular to \(x\) yields a zero overlap integral due to parity, while a non-zero result is generally obtained when the bond develops parallel to the \(x\) direction. The former limit case is named a \(s\)-\(p\)-\(\sigma\) integral, with value \(V(sp\sigma)\), since it involves an \(s\) orbital, a \(p\) orbital, and both wavefunctions are rotationally invariant (\(m=0\)) along the bond axis.The key observation here is that \(n_x\phi_{p_x}+n_y\phi_{p_y}+n_z\phi_{p_z}\) is a \(p\) orbital oriented in the \(\hat{\mathbf{n}}\) direction, thus an overlap integral in a generic orientation \(\hat{\mathbf{n}}\) can be easily calculated by projecting \(p_x\) in the bond direction, leading to the final result \(n_xV(sp\sigma)\). A similar reasoning can be applied to any possible combiantion of \(s\) and \(p\) orbitals, using the four independent parameters (Figure 4.4) corresponding to configurations in Figure 4.4:
\[V(ss\sigma), \quad V(sp\sigma), \quad V(pp\sigma), \quad V(pp\pi).\]
Note the parameters are classified based on the two overlapping orbitals and on their angular momentum along the bond axis, with \(\sigma\) for \(m=0\) (rotationally invariant), \(\pi\) for \(m=\pm 1\) and, less commonly, \(\delta\) for \(m=\pm 2\).
Using a similar approach, the hopping integrals at generic orientation satisfy
\[\begin{aligned} \langle s|V_a|s\rangle &= \int\phi^*_s(\mathbf{r})V_a(\mathbf{r}-\mathbf{R})\phi_s(\mathbf{r}-\mathbf{R})d^3r = V(ss\sigma), \\ \langle s|V_a|p_x\rangle &= \int\phi^*_s(\mathbf{r})V_a(\mathbf{r}-\mathbf{R})\phi_{p_x}(\mathbf{r}-\mathbf{R})d^3r = n_x\,V(sp\sigma), \\ \langle p_x|V_a|p_x\rangle &= \int\phi^*_{p_x}(\mathbf{r})V_a(\mathbf{r}-\mathbf{R})\phi_{p_x}(\mathbf{r}-\mathbf{R})d^3r = n_x^2\,V(pp\sigma) + (1-n_x^2)\,V(pp\pi), \\ \langle p_x|V_a|p_y\rangle &= \int\phi^*_{p_x}(\mathbf{r})V_a(\mathbf{r}-\mathbf{R})\phi_{p_y}(\mathbf{r}-\mathbf{R})d^3r = n_x n_y\,(V(pp\sigma) - V(pp\pi)). \end{aligned} \tag{4.6}\]
The angular factors of Equation 4.6 can be explored interactively in Figure 4.5, where for simplicity \(V_a(\mathbf{r})=1\) was assumed. As the second site moves around the first, each pair of \(s\), \(p_x\), \(p_y\) orbitals traces the characteristic \(n_x\), \(n_y\) combination of the table.
4.3 Graphene
Graphene is a classic TB textbook example. It incorporates all the key ingredients (a basis of two atoms, multiple orbitals, symmetry-enforced band crossings) while remaining analytically tractable in 2D.
4.3.1 Crystal structure and Brillouin zone
The honeycomb lattice is not a Bravais lattice. It is described by a triangular Bravais lattice with a two-atom basis (\(A\) and \(B\) sublattices). The primitive vectors are
\[\mathbf{t}_{1,2} = \frac{a}{2}\left(\pm 1,\sqrt{3}\right), \tag{4.7}\]
with \(a = 2.46\,\text{Å}\) (C–C bond is \(a_0 = a/\sqrt{3} = 1.42\,\text{Å}\)). The basis vectors are \(\mathbf{d}_A = (0,0)\) and \(\mathbf{d}_B = (0, a/\sqrt{3})\). The reciprocal lattice vectors are \(\mathbf{g}_{1,2} = (2\pi/a)(\pm 1, 1/\sqrt{3})\), generating a hexagonal Brillouin zone. The two physically inequivalent corner points are
\[K, K' = \frac{4\pi}{3a}\left(\pm 1, 0\right), \tag{4.8}\]
known as the Dirac points, appearing in three copies at the corners of the BZ. The midpoints of the BZ edges are saddle points \(Q\). The direct and reciprocal lattices are collected in Figure 4.6.
4.3.2 \(\sigma\) and \(\pi\) bands
Graphene has a non-trivial basis and we have to update our approach: we ned to build one Bloch sum per orbital and per basis atom or, as we call it here, per sublattice.
Definition 4.4: Bloch sums (with two sublattices)
For every atomic orbital \(\phi_i\) we build two Bloch sums, one anchored on sublattice \(A\) and one on \(B\):
\[ \begin{aligned} \psi_{i,A}(\mathbf{k},\mathbf{r}) &= \frac{1}{\sqrt N}\sum_{\mathbf{R}\in\mathcal{BL}} e^{i\mathbf{k}\cdot\mathbf{R}}\,\phi_i(\mathbf{r}-\mathbf{R}-\mathbf{d}_A), \\ \psi_{i,B}(\mathbf{k},\mathbf{r}) &= \frac{1}{\sqrt N}\sum_{\mathbf{R}\in\mathcal{BL}} e^{i\mathbf{k}\cdot\mathbf{R}}\,\phi_i(\mathbf{r}-\mathbf{R}-\mathbf{d}_B), \end{aligned} \tag{4.9}\]
with \(\mathbf{d}_A=(0,0)\), \(\mathbf{d}_B=(0,a/\sqrt3)\) and the sum \(\mathbf{R}\) running over the Bravais lattice only.
CautionThe two-center approximation with a basis
The full extension of the TB formalism to non-trivial basis is left to specialized courses/books. Just note that the two-center approximation and subsequent potential decomposition \(V(\mathbf{r})=V_a(\mathbf{r})+V'(\mathbf{r})\) of Section 4.2.1 is now sublattice-dependent and can leave some ambiguity on the off-diagonal elements of \(\mathcal{M}\) (should we decompose according to the bra or to the ket?).
This is not a problem in the specific case of graphene studied here, where the two basis atoms are identical and symmetry naturally resolves the ambiguity. Furthermore, in the spirit of the semi-empirical TB, the issue is bypassed entirely: the method aims at building an analytical model where ambiguous integrals are simply replaced by phenomenological parameters fitted to experimental or ab-initio data.
Each carbon atom has two core electrons in the \(1s\) orbitals, plus four valence electrons in the \(2s\) and \(2p\) orbitals. Core electrons are not included in the model since they form very flat and deep bands, with a small mixing. The other four orbitals on the two sublattices give \(4\times 2 = 8\) distinct Bloch sums, so \(\mathcal{M}(\mathbf{k})\) is \(8\times 8\). However, the graphene plane is a mirror symmetry plane: \(s\), \(p_x\), \(p_y\) are even under reflection, while \(p_z\) is odd. Since the Hamiltonian respects this symmetry, the matrix decouples into:
- A \(6\times 6\) block for the even (\(\sigma\)) orbitals: these form deep \(\sigma\) bands (the \(sp^2\) backbone of graphene) and their antibonding \(\sigma^*\) counterparts. All \(\sigma\) bands are fully occupied.
- A \(2\times 2\) block for the odd (\(\pi\)) orbitals: the \(p_z\) Bloch sums on sublattices \(A\) and \(B\) produce the \(\pi\) (bonding) and \(\pi^*\) (antibonding) bands that control charge conduction.
Figure 4.7 shows this block-diagonal structure explicitly. The parity argument is illustrated in the bottom-left panel of Figure 4.7: since a Bloch sum inherits the \(z\)-parity of its atomic orbital, and the Hamiltonian commutes with \(\sigma_h\), matrix elements between even and odd Bloch sums vanish identically. As a result, the solutions of the \(6\times 6\) block (\(\sigma\) bands) and \(2\times 2\) block (\(\pi\) bands) are independent with no anticrossing and can be studied independently.
We now turn to the filling of this band dispersion: there are eight valence electrons per unit cell, which can be hosted in the bottom four bands highlighted in Figure 4.7. The \(\sigma\) bands are fully occupied and are mainly responsible for the mechanical backbone of graphene. Differnetly the \(\pi\) block is half-filled and is responsible for the electronic properties of this material.
4.3.3 \(\pi\)-band dispersion
Electrons in the \(\pi\) bands in TB are described in terms of a superposition of two Bloch sums
\[\psi(\mathbf{k},\mathbf{r}) = c_A\,\psi_{p_z,A}(\mathbf{k},\mathbf{r}) + c_B\,\psi_{p_z,B}(\mathbf{k},\mathbf{r}), \tag{4.10}\] p where \(c_{A/B}\) are the complex amplitudes on the two sublattices, satisfying \(|c_A|^2 + |c_B|^2 = 1\). Restricting hopping to the nearest neighbors, as visible in Figure 4.8 (see optional sections for further details about next order corrections), the interaction matrix is relatively simple
\[\begin{aligned} \mathcal{M}_{AA}(\mathbf{k}) &= \varepsilon_{p_z}\\ \mathcal{M}_{AB}(\mathbf{k}) &= \langle\psi_{p_z,A}|\mathcal{H}|\psi_{p_z,B} \rangle = sum_{\mathbf{t}_I\in\mathcal{N\!N}}\langle\phi_{p_z}(\mathbf{r})|V_a(\mathbf{r}-\mathbf{t}_I)|\phi_{p_z}(\mathbf{r}-\mathbf{t}_I)\rangle\,e^{i\mathbf{k}\cdot\mathbf{t}_I} \\ \end{aligned}\]
where \(\varepsilon_{p_z}\) is the on-site energy of the \(p_z\) orbital, which will be taken as the energy zero. The missing matrix elements are easy to derive from symmetry: \(\mathcal{M}_{BA}=\mathcal{M}_{AB}^*\) and \(\mathcal{M}_{BB}=\mathcal{M}_{AA}\). To complete the calculate, we note that each \(A\) atom is bonded to three \(B\) atoms through the vectors (see Figure 4.8)
\[\boldsymbol{\tau}_1 = \mathbf{d}_B,\qquad \boldsymbol{\tau}_2 = \mathbf{d}_B - \mathbf{t}_1,\qquad \boldsymbol{\tau}_3 = \mathbf{d}_B - \mathbf{t}_2,\qquad |\boldsymbol{\tau}_j| = \frac{a}{\sqrt3}, \tag{4.11}\]
which set \(\mathbf{t}_I = \mathbf{0},\,-\mathbf{t}_1,\,-\mathbf{t}_2\). Given \(p_z\) orbitals are rotationally symmetric in the graphene plane and given all A-B bonds have the same length, the overlap integral is always the same
\[\big\langle\phi_{p_z}(\mathbf{r})\big|V_a(\mathbf{r}-\boldsymbol{\tau}_j)\big|\phi_{p_z}(\mathbf{r}-\boldsymbol{\tau}_j)\big\rangle \;=\; V(pp\pi)\;\equiv\;t , \tag{4.12}\]
and the interaction matrix term is
\[\mathcal{M}_{AB}(\mathbf{k}) =\; V(pp\pi)\sum_{\mathbf{t}_I\in\mathcal{N\!N}} e^{i\mathbf{k}\cdot\mathbf{t}_I} = V(pp\pi)\Bigl(1 + e^{-i\mathbf{k}\cdot\mathbf{t}_1} + e^{-i\mathbf{k}\cdot\mathbf{t}_2}\Bigr)\;\equiv\;V(pp\pi)\,F(\mathbf{k}), \tag{4.13}\]
where, using the primitive vectors in Equation 4.7, the geometric structure factor \(F(\mathbf{k})\) is
\[F(\mathbf{k}) = 1 + 2\cos\!\left(\frac{k_x a}{2}\right)\exp\left(-i\frac{k_y a\sqrt{3}}{2}\right). \tag{4.14}\]
Finally, the secular equation for the \(\pi\) block is
\[ V(pp\pi) \begin{bmatrix} 0 & F(\mathbf{k}) \\ F^*(\mathbf{k}) & 0 \end{bmatrix} \begin{bmatrix} c_A\\c_B \end{bmatrix} = E(\mathbf{k}) \begin{bmatrix} c_A\\c_B \end{bmatrix}, \tag{4.15}\]
with eigenvalues \(E_\pm(\mathbf{k}) = \pm|V(pp\pi)|\,|F(\mathbf{k})|\). This symmetric dispersion is visible in Figure 4.8, where the valence/conduction bands are also called \(\pi\)/\(\pi^*\) bands.
CautionParticle–hole asymmetry in graphene
In the nearest-neighbor model just discussed, the spectrum is exactly particle–hole symmetric, \(E_\pm = \pm|V(pp\pi)||F|\). In real graphene the conduction band has a broader energy span with respect to the valence one. There are two equivalent ways to account for it.
(i) Orthogonality of Bloch sums. This is the origin of the asymmetry in Bassani and Pastori Parravicini (1967), whose bands are reproduced in Figure 4.7: the \(p_z\) orbitals of neighboring atoms do overlap, and \(\langle\phi_{p_z}(\mathbf{r})|\phi_{p_z}(\mathbf{r}-\boldsymbol{\tau})\rangle \neq 0\), thus \(\mathcal{S}\neq\mathbb{1}\) and this leads to corrections in the \(2\times2\) generalized problem Equation 4.4. As a rough explanation, the overlap in the valence band is larger than in the conduction one, due to the more consistent phase of nearby \(A\) and \(B\) sites (see also Figure 4.9). In the generalized diagonalization, this enhances the conduction band eigenenergies and reduces the valence ones.
(ii) Second neighbors. In more modern descriptions, asymmetry is ascribed to second-order hopping. The six second neighbors of an atom sit on the same sublattice, so they lead to a diagonal correction of the \(\pi\) block, adding the same \(\mathbf{k}\)-dependent term to both bands:
\[E_\pm(\mathbf{k}) = \pm|t|\,|F(\mathbf{k})| + t'\bigl(|F(\mathbf{k})|^2 - 3\bigr), \tag{4.16}\]
where the constant \(-3t'\) keeps the Dirac point at \(E=0\). This is what the NN/NNN selector of Figure 4.8 switches on: the \(\pi^*\) bandwidth grows while the \(\pi\) one shrinks, and the second-order term also adds an isotropic \(q^2\) correction to the cone.
Show proof
The six second neighbors of an \(A\) site are the six \(A\) sites of the surrounding cells, reached by the lattice vectors
\[\pm\mathbf{t}_1,\qquad \pm\mathbf{t}_2,\qquad \pm(\mathbf{t}_1-\mathbf{t}_2),\]
all at the same distance \(a\), hence all carrying the same hopping integral \(t'\). Since they sit on the same sublattice, they contribute to the diagonal entries of the \(\pi\) block — with the same lattice sum on \(A\) and on \(B\), because the two sublattices see identical second-neighbor shells:
\[\mathcal{M}_{AA}(\mathbf{k}) = \mathcal{M}_{BB}(\mathbf{k}) = \varepsilon_{2p} + t'\,f(\mathbf{k}), \qquad f(\mathbf{k}) = 2\Bigl[\cos(\mathbf{k}\!\cdot\!\mathbf{t}_1) + \cos(\mathbf{k}\!\cdot\!\mathbf{t}_2) + \cos\bigl(\mathbf{k}\!\cdot\!(\mathbf{t}_1-\mathbf{t}_2)\bigr)\Bigr],\]
where opposite vectors have been paired into cosines. The off-diagonal elements are untouched, so the block reads
\[ \begin{bmatrix} \varepsilon_{2p} + t'f(\mathbf{k}) & t\,F(\mathbf{k}) \\ t\,F^*(\mathbf{k}) & \varepsilon_{2p} + t'f(\mathbf{k}) \end{bmatrix} \begin{bmatrix} c_A\\ c_B \end{bmatrix} = E(\mathbf{k}) \begin{bmatrix} c_A\\ c_B \end{bmatrix}, \]
a matrix of the form \((\varepsilon_{2p}+t'f)\,\mathbb{1} + t\,(\dots)\): the new term is proportional to the identity, so it does not mix the two sublattices and simply shifts both eigenvalues by the same amount,
\[E_\pm(\mathbf{k}) = \varepsilon_{2p} + t'f(\mathbf{k}) \pm |t|\,|F(\mathbf{k})|.\]
The lattice sum \(f\) is not independent of \(F\). Expanding the modulus of \(F(\mathbf{k}) = 1 + e^{-i\mathbf{k}\cdot\mathbf{t}_1} + e^{-i\mathbf{k}\cdot\mathbf{t}_2}\),
\[|F|^2 = 3 + 2\cos(\mathbf{k}\!\cdot\!\mathbf{t}_1) + 2\cos(\mathbf{k}\!\cdot\!\mathbf{t}_2) + 2\cos\bigl(\mathbf{k}\!\cdot\!(\mathbf{t}_1-\mathbf{t}_2)\bigr) = 3 + f(\mathbf{k}),\]
the three cross terms being exactly the three cosines of \(f\). Substituting \(f = |F|^2-3\) and taking \(\varepsilon_{2p}=0\) gives Equation 4.16, whose energies are automatically measured from the Dirac point: at \(K\) one has \(|F|=0\), hence \(f(K)=-3\) and \(E_\pm(K)=0\).
Two remarks. First, the correction is the same for both bands, so it deforms them rigidly and cannot open a gap: \(t'\) multiplies the identity, never \(\sigma_z\). Second, it breaks particle–hole symmetry because it is even in \(|F|\) while the band term is odd: near \(\Gamma\), where \(|F|=3\), the shift is \(+6t'\) for both bands, so the \(\pi^*\) band is pushed up and the \(\pi\) band, already below, is pushed towards \(E_F\) — the antibonding band widens and the bonding one narrows. \(\blacksquare\)
4.3.4 Linearization at the Dirac point
Everything insresting in graphene typically happens within a few hundred meV from the Dirac point, so we now expand \(F(\mathcal{k})\) close to the Dirac point setting \(\mathbf{k} = \mathbf{K} + \mathbf{q}\) with \(|\mathbf{q}|\ll|\mathbf{K}|\). Using \(K_x a/2 = 2\pi/3\), we obtain
\[F(\mathbf{K}+\mathbf{q}) = 1 + 2\cos\!\left(\frac{2\pi}{3}+\frac{q_xa}{2}\right)e^{-iq_ya\sqrt3/2} \simeq -\frac{a\sqrt{3}}{2}\,(q_x - iq_y), \tag{4.17}\]
which is linear in \(\mathbf{q}\). Substituting into Equation 4.15, the \(\pi\) block becomes the celebrated Dirac Hamiltonian
\[\mathcal{H}_K = \hbar v_F\,\boldsymbol{\sigma}\cdot\mathbf{q} = \hbar v_F\begin{bmatrix} 0 & q_x - iq_y \\ q_x + iq_y & 0 \end{bmatrix}, \tag{4.18}\]
where \(\boldsymbol{\sigma} = (\sigma_x,\sigma_y)\) are Pauli matrices acting in the sublattice space, and
\[\hbar v_F = |V(pp\pi)|\,\frac{a\sqrt{3}}{2}\qquad\Longrightarrow\qquad v_F \approx 10^{6}\,\text{m/s} \approx \frac{c}{300} \tag{4.19}\]
defines the Fermi velocity. The eigenvalues for a linear conical dispersion
\[E_\pm(\mathbf{q}) = \pm\hbar v_F|\mathbf{q}|, \tag{4.20}\]
which is identical in form to the relativistic relation \(E = cp\) of a massless particle, with \(v_F\) in place of the speed of light.
CautionThe symmetric \(K'\) valley
Repeating the expansion at the other inequivalent corner, \(\mathbf{K}' = -\mathbf{K}\), a very similar result except for a complex conjugation. The linearized Hamiltonian of the second valley is in fact
\[\mathcal{H}_{K'} = \hbar v_F\,\boldsymbol{\sigma}^{*}\!\cdot\mathbf{q} = \hbar v_F\,(\sigma_x q_x - \sigma_y q_y) = \mathcal{H}_K^{*} . \tag{4.21}\]
with identical eigenvalues \(E_\pm = \pm\hbar v_Fq\), but a different thepseudospin texture. This is can be seen a s time-reversed copy of the state a the \(K\) Dirac point.
CautionFurther corrections: trigonal warping
The cone is only the leading term. The expansion at the second order leads to,
\[E_\pm(\mathbf{q}) \simeq \pm\hbar v_F q\left[1 - \frac{a\sqrt3}{12}\,q\,\cos 3\varphi_{\mathbf{q}}\right]. \tag{4.22}\]
where the \(\cos3\varphi_{\mathbf{q}}\) dependence turns the circular constant-energy contours into rounded triangles. This is trigonal warping, and it could not have been otherwise — the symmetry that survives at \(K\) is the three-fold \(C_{3v}\), not the full rotational invariance that the linearized model accidentally displays.
4.3.5 Pseudospin
The two-component structure of Equation 4.18 is mathematically equivalent to a spin-1/2 degree of freedom.
Definition 4.5: Pseudospin
The eigenvector \((c_A,c_B)\) of Equation 4.10 is also called pseudospin, but it has nothing to do with the real spin of the electron, and its components only quantify the amplitudes on the two sublattices. The eigenstates of \(\mathcal{H}_K = \hbar v_F\,\boldsymbol{\sigma}\cdot\mathbf{q}\) are
\[\psi_\pm(\mathbf{q}) = \frac{\psi_{p_z,A}(\mathbf{K}+\mathbf{q})\pm e^{i\varphi_{\mathbf{q}}}\psi_{p_z,B}(\mathbf{K}+\mathbf{q})}{\sqrt2} = \cos\frac{\theta}{2}\,|\!\uparrow\rangle + e^{i\varphi}\sin\frac{\theta}{2}\,|\!\downarrow\rangle, \tag{4.23}\]
having identified the two sublattice states with the two pseudospin projections, \(|\!\uparrow\rangle \equiv \psi_{p_z,A}\) and \(|\!\downarrow\rangle \equiv \psi_{p_z,B}\), and written the spinor in the usual Bloch-sphere form, with \(\varphi = \varphi_{\mathbf{q}}\) in the conduction band and \(\varphi = \varphi_{\mathbf{q}}+\pi\) in the valence one.
Here the polar angle is always \(\theta = \pi/2\): the two sublattices carry equal weight, \(\cos(\theta/2) = \sin(\theta/2) = 1/\sqrt2\), so the pseudospin never tilts out of the graphene plane and only its azimuth \(\varphi\) matters — it is locked to the direction of propagation, parallel to \(\mathbf{q}\) in the conduction band (\(E>0\)) and antiparallel in the valence one (\(E<0\)). Its projection along the momentum, \(h = \boldsymbol{\sigma}\cdot\hat{\mathbf{q}} = \pm1\), is a conserved quantum number called pseudo-helicity (or chirality).
CautionAnalogy to a Zeeman hamiltonian
Equation 4.18 has exactly the structure of a Zeeman Hamiltonian for a spin-1/2 in a magnetic field,
\[\mathcal{H}_\text{Zeeman} = \frac{g\mu_B}{2}\,\boldsymbol{\sigma}\cdot\mathbf{B} \qquad\longleftrightarrow\qquad \mathcal{H}_K = \hbar v_F\boldsymbol{\sigma}\cdot\mathbf{q},\]
with the crystal momentum playing the role of the magnetic field. In this language, the two bands are then nothing but the two Zeeman branches — pseudospin aligned (\(\pi^*\)) and anti-aligned (\(\pi\)) with the field, split by \(2\hbar v_F q\). The splitting closes only where the effective field vanishes, \(\mathbf{q}=0\). This analogy can be directly connected to a fundamental properties of electrons at the Dirac point which will be later discussed in the course and goes under the name of Berry phase.
CautionPhase texture in the real space
To make the real-space meaning of a Bloch state more vivid, Figure 4.9 shows the phase of the actual \(\pi\) eigenstate over the honeycomb lattice for any chosen \(\mathbf{k}\) in the BZ. Each \(A\) site carries the cell factor \(e^{i\mathbf{k}\cdot\mathbf{R}}\); each \(B\) site carries in addition the pseudospin phase of Equation 4.23, \(c_B/c_A = \mp e^{-i\varphi_{\mathbf{k}}}\) — in phase with \(A\) at \(\Gamma\) for the bonding \(\pi\) band, in antiphase for the antibonding \(\pi^*\). The phase is encoded as a cyclic color wheel (hue \(= \arg\psi\), full saturation). Faint gray stripes across the lattice mark the wavefronts — surfaces of constant cell phase, perpendicular to \(\mathbf{k}\) and spaced by \(2\pi/|\mathbf{k}|\).
Pick a \(\mathbf{k}\) close to \(K\) and then the opposite one, \(\mathbf{q}\to-\mathbf{q}\) at the same energy, and compare the two patterns — switching the reference to \(q\) makes the comparison easier, since only the envelope is left on the screen. In an ordinary crystal the two counter-propagating states are related by a trivial operation: same amplitudes, wavefronts running the other way. Here something more happens. Reversing \(\mathbf{q}\) rotates \(\varphi_{\mathbf{q}}\) by \(\pi\), so by Equation 4.23 the relative phase between the two sublattices flips as well: the pseudospin must be turned upside down, from \(+\hat{\mathbf{q}}\) to \(-\hat{\mathbf{q}}\). Going backwards is not a mere change of sign of the wavevector — it requires a wholesale rearrangement of the internal phases of the wavefunction, sublattice by sublattice.
This is the microscopic reason behind the peculiar conduction properties of graphene. A smooth scattering potential — a remote charged impurity, a long-wavelength strain, anything varying slowly on the scale of \(a\) — acts identically on the two sublattices, i.e. it is proportional to the identity in pseudospin space and cannot flip the pseudospin. Direct backscattering (\(\mathbf{q}\to-\mathbf{q}\)) is thus forbidden, and the scattering probability vanishes as \(\cos^2(\vartheta/2)\) for a deflection by \(\vartheta\). Carriers keep moving forward, which is why mobilities in clean graphene are enormous and why a normally incident electron crosses a smooth potential barrier with probability one (Klein tunnelling). Only sharp, atomic-scale defects — which do distinguish \(A\) from \(B\), and can also scatter between the two valleys — restore ordinary backscattering.
CautionIs the Dirac point an artefact of our approximations?
We reached the cone through a chain of drastic simplifications — orthogonal orbitals, two-center integrals only, nearest neighbors, a single \(p_z\) per atom. It is legitimate to suspect that some neglected term, at some higher order, could open a small gap at \(K\). It cannot, and the reason is symmetry rather than accuracy.
Look at the \(2\times2\) block: any gap at the Dirac point requires a term \(\propto\sigma_z\), i.e. different on-site energies on the two sublattices,
\[\mathcal{H} = \hbar v_F\,\boldsymbol{\sigma}\cdot\mathbf{q} + \Delta\,\sigma_z \quad\Longrightarrow\quad E_\pm = \pm\sqrt{(\hbar v_Fq)^2 + \Delta^2}.\]
But in graphene \(A\) and \(B\) are the same atom in symmetry-equivalent positions: the honeycomb lattice is invariant under an inversion centered at the middle of a C–C bond, which exchanges the two sublattices. As long as that symmetry holds, \(\Delta \equiv 0\) by construction — no matter how many neighbor shells, three-center integrals or overlap corrections we add, since all of them respect the symmetry of the crystal. Higher orders can shift, warp and skew the bands (as in the two boxes above) but they can never produce the one term that would gap the crossing.
The same statement can be made in the language of group theory: at \(K\) the little group is \(C_{3v}\) (three-fold rotation plus the vertical mirrors), which possesses a two-dimensional irreducible representation \(E\); the two \(p_z\) Bloch states at \(K\) transform precisely as that doublet, so their degeneracy is symmetry-enforced, not accidental. Any Hamiltonian with the symmetry of graphene must keep it. And with time-reversal symmetry also present, the two doublets at \(K\) and \(K'\) are pinned to the same energy.
There is even a topological layer of protection: the winding of the pseudospin around the Dirac point (the \(\pi\) Berry phase) is an integer that cannot change continuously. A perturbation may move the Dirac points inside the zone — uniaxial strain, which makes the three bonds inequivalent, does exactly this — but it cannot make one of them disappear on its own: they can only be annihilated in pairs, when two of them merge, which in the honeycomb requires an unrealistically large bond anisotropy.
The gap does open, of course, when the protecting symmetry is genuinely broken: put graphene on hexagonal boron nitride or in a commensurate substrate potential that distinguishes \(A\) from \(B\) (a staggered on-site energy \(\pm\Delta\)), and a mass term appears — that is the massive Dirac fermion discussed in Section A.1. The lesson is the usual one in this course: what protects a degeneracy is symmetry, and approximations that respect the symmetry cannot lift it.
In conclusion…
ImportantTake-home message
- The Bloch sum, turning any set of localized orbitals into Bloch-compliant basis states.
- Empirical tight binding what is the key idea behind the approach? What are its approximation and their justification?
- Graphene extension to sublattice Bloch sums and derivation of the Dirac hamiltonian eigenvalue problem; definition of pseudospin.
You are not expected to remember the contents of the optional sections, but working through them is a good way to test your understanding of the material of this chapter.
Test your understanding
Below we report a few simple band models derived from an extension of our minimalistic TB model to a 2D square lattice, or squarium, and to a 3D cubic lattice, or cubium. These will be a natural playground to test your understanding on density of states and critical points.
CautionExample: 1D tight-binding chain
The simplest non-trivial DOS belongs to the 1D tight-binding chain of Section 4.1: a single \(s\)-wave orbital per site with nearest-neighbor hopping \(\gamma<0\) (on-site energy set to zero). The dispersion is
\[E(k) = 2\gamma\cos(ka), \qquad k\in[-\pi/a, \pi/a],\]
so the band spans \([-2|\gamma|, +2|\gamma|]\). The two stationary points (a minimum at \(\Gamma\) and a maximum at the zone boundary \(X = \pi/a\)) generate the canonical 1D van Hove signature: a \(1/\sqrt{|E - E_c|}\) divergence at each band edge. The DOS admits the closed form
\[D(E) = \frac{1}{\pi\sqrt{4\gamma^2 - E^2}}, \qquad |E| < 2|\gamma|.\]
Drag the iso-energy line in Figure 4.10 to see how the filled (orange) and empty (blue) parts of the band track the cumulative integral of \(D(E)\).
CautionExample: the squarium
The squarium is the natural extension of the previous pedagogical example to a 2D square lattice, we now have four nearest neighbors and the two-center integral yields (site energy set to zero again):
\[E(k_x,k_y) = 2\gamma\cos k_x a + 2\gamma\cos k_y a. \tag{4.24}\]
The band spans \([-4|\gamma|,\,+4|\gamma|]\) with:
- A minimum at \(\Gamma = (0,0)\) — the iso-energetic contour is approximately circular (parabolic dispersion), and \(D(E)\) jumps from zero to a finite value (step, an \(M_0\) singularity).
- Saddle points at \(X = (\pi/a,0)\) and \(Y = (0,\pi/a)\) — the Fermi contour touches the zone boundary and changes topology (Lifshitz transition). The DOS shows a logarithmic van Hove divergence (\(M_1\) singularity).
- A maximum at \(M = (\pi/a,\pi/a)\) — the Fermi contour is again circular (but now an empty contour enclosing unoccupied states), and \(D(E)\) drops back to zero (step, an \(M_2\) singularity).
CautionExample: the cubium
A further extension to 3D leads to the cubium, with an obvious dispersion
\[E(\mathbf{k}) = 2\gamma\cos k_x a + 2\gamma\cos k_y a + 2\gamma\cos k_z a.\]
The density of states can again be expressed in closed form, as a convolution of the 2D result with a one-dimensional density of states. The resulting \(D(E)\) displays two square-root van Hove cusps at \(E = \pm 2|\gamma|\) — in 3D the DOS stays finite at the saddle points, with a divergent one-sided slope — and square-root onsets at the band edges \(E = \pm 6|\gamma|\).