Chapter 3

Topological Insulators

Time-Reversal Symmetry & Kramers Theorem

Topological insulators are protected by time-reversal symmetry (TRS). The anti-unitary time-reversal operator for spin-1/2 particles satisfies:

$$\Theta = i\sigma_y K, \qquad \Theta^2 = -1$$

where $K$ is complex conjugation. The condition $\Theta^2 = -1$ for half-integer spin leads to Kramers theorem: every energy eigenstate at a time-reversal invariant momentum (TRIM) must be at least two-fold degenerate.

Kramers Theorem

For a system with $\Theta^2 = -1$, the states $|\psi\rangle$ and$\Theta|\psi\rangle$ are orthogonal and degenerate. These Kramers pairscannot be split by any time-reversal invariant perturbation.

$$\langle \psi | \Theta\psi \rangle = 0 \quad \text{(orthogonality of Kramers partners)}$$

At the TRIM points $\mathbf{k}_i$ satisfying $\mathbf{k}_i = -\mathbf{k}_i + \mathbf{G}$ (modulo a reciprocal lattice vector), the Bloch Hamiltonian commutes with $\Theta$, enforcing Kramers degeneracy. The way these degeneracies connect across the Brillouin zone distinguishes topologically trivial from nontrivial insulators.

Z₂ Topological Invariant

Unlike the quantum Hall effect with its $\mathbb{Z}$-valued Chern number, time-reversal invariant topological insulators are classified by a $\mathbb{Z}_2$ invariant $\nu \in \{0, 1\}$. The value $\nu = 1$ signals a topological insulator, while $\nu = 0$ is trivial.

$$(-1)^\nu = \prod_{i=1}^{4} \frac{\text{Pf}[w(\Gamma_i)]}{\sqrt{\det[w(\Gamma_i)]}}$$

Here $w_{mn}(\mathbf{k}) = \langle u_m(-\mathbf{k}) | \Theta | u_n(\mathbf{k}) \rangle$ is the sewing matrix. The Pfaffian ratio gives $\pm 1$ at each TRIM, and the product over all four TRIM points in 2D determines $\nu$.

For systems with inversion symmetry (Fu-Kane formula)

$$(-1)^\nu = \prod_{i=1}^{4} \prod_{m=1}^{N} \xi_{2m}(\Gamma_i)$$

where $\xi_{2m}(\Gamma_i)$ is the parity eigenvalue of the $2m$-th occupied band at TRIM point $\Gamma_i$. This enables straightforward computation from band structure calculations.

The $\mathbb{Z}_2$ invariant is robust under adiabatic deformations preserving the bulk gap and TRS. Changing $\nu$ requires closing and reopening the gap -- a topological phase transition.

Kane-Mele Model

The Kane-Mele model extends graphene with intrinsic spin-orbit coupling, providing the first model of a 2D topological insulator:

$$H = -t \sum_{\langle ij \rangle} c_i^\dagger c_j + i\lambda_{\text{SO}} \sum_{\langle\langle ij \rangle\rangle} \nu_{ij}\, c_i^\dagger s_z c_j + i\lambda_R \sum_{\langle ij \rangle} c_i^\dagger (\mathbf{s} \times \hat{\mathbf{d}}_{ij})_z c_j$$

The terms are: (1) NN hopping $t$, (2) intrinsic SOC $\lambda_{\text{SO}}$ between NNN sites with chirality $\nu_{ij} = \pm 1$, and (3) Rashba SOC $\lambda_R$ breaking $s_z$ conservation.

Topological Gap

The spin-orbit term opens a gap at the Dirac points K and K' of graphene:

$$\Delta = 6\sqrt{3}\,\lambda_{\text{SO}}$$

Topological Phase

The system is topological when the SOC-induced gap dominates over Rashba:

$$\nu = 1 \quad \text{for} \quad \lambda_R < 2\sqrt{3}\,\lambda_{\text{SO}}$$

For graphene the intrinsic SOC is extremely weak ($\sim 10\,\mu$eV), but the Kane-Mele framework was soon realized in HgTe quantum wells.

Quantum Spin Hall Effect

A 2D topological insulator with $\nu = 1$ exhibits the quantum spin Hall (QSH) effect. The bulk is insulating, but the edges carry helical edge states: counter-propagating modes with opposite spin polarization.

$$\sigma_{xy}^{\text{spin}} = \frac{e}{2\pi}, \qquad G = \frac{2e^2}{h}$$

The edge states form a single Kramers pair: a right-mover with spin-up and a left-mover with spin-down (or vice versa on the opposite edge). Time-reversal symmetry forbids elastic backscattering between these states, since any scattering event would need to flip the spin -- violating TRS for non-magnetic impurities.

Helical Edge States

The edge dispersion near the crossing point has the form:

$$H_{\text{edge}} = v_F k \sigma_z$$

The right-mover carries spin-up and the left-mover carries spin-down, forming a 1D helical liquid. Unlike the chiral edge states of the quantum Hall effect, helical states come in time-reversed pairs and cannot be fully gapped without breaking TRS.

The quantized conductance $G = 2e^2/h$ arises from one Kramers pair of edge modes, each contributing $e^2/h$. This was dramatically confirmed in HgTe/CdTe quantum wells by Konig et al. (2007), validating the Bernevig-Hughes-Zhang (BHZ) model prediction.

3D Topological Insulators

In three dimensions, the $\mathbb{Z}_2$ classification involves four invariants$(\nu_0;\, \nu_1\nu_2\nu_3)$. The strong topological insulator has $\nu_0 = 1$, featuring robust metallic surface states on all surfaces.

$$H_{\text{surface}} = v_F(\hat{z} \times \boldsymbol{\sigma}) \cdot \mathbf{k} = v_F(k_x \sigma_y - k_y \sigma_x)$$

The surface states form a single Dirac cone -- an odd number of Dirac cones that cannot exist in a purely 2D system (the fermion doubling theorem would require an even number). This is the hallmark of a 3D topological insulator.

Spin-Momentum Locking

The spin is locked perpendicular to the momentum in the surface plane:$\langle \mathbf{s} \rangle = \frac{\hbar}{2}\frac{\hat{z} \times \mathbf{k}}{|\mathbf{k}|}$. This creates a helical spin texture around the Dirac cone, directly measurable by spin-resolved ARPES.

Topological Protection

The Dirac point is pinned to a TRIM. Backscattering ($\mathbf{k} \to -\mathbf{k}$) is forbidden by TRS since it requires a $\pi$ spin rotation. The Berry phase around the Fermi surface is $\pi$, suppressing weak localization and producing weak anti-localization instead.

Topological Insulator Materials

Several material families with strong spin-orbit coupling host topological insulator phases:

Material	Dimension	Bulk Gap	Key Feature
HgTe/CdTe QW	2D	~10 meV	First experimental QSH system (2007)
InAs/GaSb QW	2D	~4 meV	Type-II quantum well, tunable gap
Bi$_2$Se$_3$	3D	~0.3 eV	Single Dirac cone, large gap
Bi$_2$Te$_3$	3D	~0.17 eV	Warped Dirac cone, thermoelectric
Sb$_2$Te$_3$	3D	~0.28 eV	p-type bulk, single Dirac cone
Bi$_{1-x}$Sb$_x$	3D	~30 meV	First 3D TI (2008), complex surface

Bi$_2$Se$_3$ is considered the "hydrogen atom" of 3D topological insulators: it has a single, well-isolated Dirac cone at the $\Gamma$ point, a relatively large bulk gap (~0.3 eV), and a simple rhombohedral crystal structure. The band inversion occurs between Se $p$ and Bi $p$ orbitals near the $\Gamma$ point.

Experimental Signatures

ARPES (Angle-Resolved Photoemission)

ARPES directly images the surface electronic structure, revealing the Dirac cone dispersion $E = \pm v_F |\mathbf{k}|$ and the spin texture via spin-resolved measurements. The linear crossing at the $\Gamma$ point with a single Dirac cone is the definitive signature of a strong 3D TI.

Transport Measurements

In 2D TIs, the quantized edge conductance $G = 2e^2/h$ is measured in short channels. In 3D TIs, signatures include weak anti-localization (WAL), Shubnikov-de Haas oscillations with a $\pi$ Berry phase shift, and Aharonov-Bohm oscillations in nanowires confirming the surface state nature.

STM & Quasiparticle Interference

Scanning tunneling microscopy reveals the absence of backscattering from non-magnetic impurities. The Fourier transform of quasiparticle interference patterns shows the suppression of $\mathbf{q} = 2\mathbf{k}_F$ scattering vectors, confirming the topological protection of surface states.

Majorana Fermions & Topological Superconductors

The interface between a topological insulator and an s-wave superconductor realizes an effective topological superconductor. The proximity-induced pairing on the TI surface state takes the form:

$$H = v_F(\mathbf{k} \times \boldsymbol{\sigma})\cdot\hat{z} - \mu + \Delta(k)\tau_x$$

A vortex in the superconducting order parameter $\Delta$ on the TI surface binds a Majorana zero mode -- a zero-energy state that is its own antiparticle: $\gamma = \gamma^\dagger$. These non-Abelian anyons obey exchange statistics that could enable fault-tolerant topological quantum computation.

The Fu-Kane proposal (2008) showed that the surface of a 3D TI proximitized by an s-wave superconductor mimics a spinless $p_x + ip_y$ superconductor, but without requiring spin-triplet pairing. Each vortex core hosts a single Majorana bound state, and braiding these vortices performs topologically protected quantum gates.

Periodic Table of Topological Insulators

The complete classification of free-fermion topological phases depends on spatial dimension$d$ and the discrete symmetries: time-reversal $\Theta$, particle-hole $\Xi$, and chiral $\Pi = \Theta\Xi$. The ten Altland-Zirnbauer symmetry classes yield the following topological invariants:

Class	$\Theta$	$\Xi$	$\Pi$	d=1	d=2	d=3
A (Unitary)	0	0	0	--	$\mathbb{Z}$	--
AIII (Chiral)	0	0	1	$\mathbb{Z}$	--	$\mathbb{Z}$
AII (TI)	-1	0	0	--	$\mathbb{Z}_2$	$\mathbb{Z}_2$
AI	+1	0	0	--	--	--
BDI	+1	+1	1	$\mathbb{Z}$	--	--
D (SC)	0	+1	0	$\mathbb{Z}_2$	$\mathbb{Z}$	--
DIII	-1	+1	1	$\mathbb{Z}_2$	$\mathbb{Z}_2$	$\mathbb{Z}$
C	0	-1	0	--	$\mathbb{Z}$	--
CI	+1	-1	1	--	--	$\mathbb{Z}$
CII	-1	-1	1	$\mathbb{Z}$	--	$\mathbb{Z}_2$

The highlighted AII row contains the time-reversal invariant topological insulators ($\mathbb{Z}_2$ in 2D and 3D). The D row describes topological superconductors hosting Majorana modes. This "periodic table" exhibits Bott periodicity with period 8 in dimension, unifying all topological phases of non-interacting fermions.

Detailed Derivations

Derivation 1: Kramers Theorem from $\Theta^2 = -1$

We prove that every eigenstate of a time-reversal invariant Hamiltonian with half-integer spin is at least two-fold degenerate.

Step 1: For spin-1/2 particles, the time-reversal operator is $\Theta = i\sigma_y K$ where $K$ is complex conjugation. Compute $\Theta^2$:

$$\Theta^2 = (i\sigma_y K)(i\sigma_y K) = i\sigma_y \cdot (i\sigma_y K) \cdot K = i\sigma_y \cdot (-i\sigma_y^*) = i\sigma_y \cdot (-i\sigma_y) = -\sigma_y^2 = -\mathbf{1}$$

where we used $K(i\sigma_y) = -i\sigma_y^* = -i\sigma_y$ (since $\sigma_y$ is purely imaginary: $\sigma_y^* = -\sigma_y$).

Step 2: Let $|\psi\rangle$ be an eigenstate of $H$ with energy $E$. Since $[H, \Theta] = 0$, the state $\Theta|\psi\rangle$ is also an eigenstate with the same energy $E$.

Step 3: We show $|\psi\rangle$ and $\Theta|\psi\rangle$ are orthogonal. Compute the overlap using the anti-unitary property $\langle \Theta\alpha|\Theta\beta\rangle = \langle\beta|\alpha\rangle$:

$$\langle\psi|\Theta\psi\rangle = \langle\Theta(\Theta\psi)|\Theta\psi\rangle^* = \langle\Theta^2\psi|\Theta\psi\rangle^*$$

Step 4: Since $\Theta^2 = -1$, we have $\Theta^2|\psi\rangle = -|\psi\rangle$, so:

$$\langle\psi|\Theta\psi\rangle = \langle -\psi|\Theta\psi\rangle^* = -\langle\psi|\Theta\psi\rangle^*$$

$$\langle\psi|\Theta\psi\rangle = 0$$

Therefore $|\psi\rangle$ and $\Theta|\psi\rangle$ are orthogonal and degenerate -- this is the Kramers doublet. No time-reversal invariant perturbation can split this degeneracy.

Derivation 2: Surface Dirac Cone Hamiltonian from $\mathbf{k}\cdot\mathbf{p}$ Theory

We derive the effective surface Hamiltonian for a 3D topological insulator by imposing time-reversal symmetry and the surface boundary condition.

Step 1: Consider the low-energy bulk Hamiltonian of a 3D TI such as Bi$_2$Se$_3$, expanded around the $\Gamma$ point. The four-band $\mathbf{k}\cdot\mathbf{p}$ model (Zhang et al., 2009) is:

$$H(\mathbf{k}) = \epsilon_0(\mathbf{k})\mathbf{1}_{4\times 4} + \begin{pmatrix} M(\mathbf{k}) & A_1 k_z & 0 & A_2 k_- \\ A_1 k_z & -M(\mathbf{k}) & A_2 k_- & 0 \\ 0 & A_2 k_+ & M(\mathbf{k}) & -A_1 k_z \\ A_2 k_+ & 0 & -A_1 k_z & -M(\mathbf{k}) \end{pmatrix}$$

where $k_\pm = k_x \pm ik_y$ and $M(\mathbf{k}) = M_0 - B_1 k_z^2 - B_2(k_x^2 + k_y^2)$. The band inversion occurs when $M_0 B_1 > 0$ or $M_0 B_2 > 0$.

Step 2: At the surface (take $z = 0$ as the surface, $z > 0$ as vacuum), replace $k_z \to -i\partial_z$. The surface state wavefunction decays exponentially into the bulk: $\psi(z) \propto e^{-\kappa z}$.

Step 3: For $k_x = k_y = 0$, the Hamiltonian decouples into two $2\times 2$ blocks. Each block is:

$$h(k_z) = A_1 k_z \sigma_x + M(k_z)\sigma_z$$

Replacing $k_z \to -i\partial_z$ and solving $h(-i\partial_z)\psi = 0$ for the zero-energy bound state with the ansatz $\psi = (e^{\lambda_1 z} - e^{\lambda_2 z})\chi$, we find:

$$\lambda_{1,2} = \frac{A_1 \pm \sqrt{A_1^2 - 4B_1 M_0}}{2B_1}$$

The surface state exists (normalizable for $z < 0$, i.e., into the bulk) when $M_0/B_1 > 0$ -- precisely the band-inversion condition.

Step 4: Project the full Hamiltonian onto the surface-state subspace at finite $k_x, k_y$. The $A_2 k_\pm$ terms couple the two blocks, yielding the $2\times 2$ effective surface Hamiltonian:

$$H_\text{surf} = A_2(k_x \sigma_y - k_y \sigma_x) = v_F(\hat{z}\times\boldsymbol{\sigma})\cdot\mathbf{k}$$

with $v_F = A_2/\hbar$. The spectrum is $E = \pm \hbar v_F |\mathbf{k}|$, a single massless Dirac cone.

Step 5: The spin-momentum locking follows directly: $\langle\boldsymbol{\sigma}\rangle = \hat{z}\times\hat{\mathbf{k}}$ for the upper cone, meaning the spin is always perpendicular to the momentum and tangent to the constant-energy contour, forming a left-handed helical texture.

Derivation 3: Z$_2$ Invariant via the Fu-Kane Parity Criterion

For centrosymmetric crystals, Fu and Kane (2007) showed that the Z$_2$ invariant can be computed simply from the parity eigenvalues of the occupied bands at TRIM points.

Step 1: At a TRIM point $\Gamma_i$ where $\mathbf{k}_i = -\mathbf{k}_i + \mathbf{G}$, the Bloch Hamiltonian commutes with both $\Theta$ (time reversal) and $P$ (spatial inversion). Since $[P, \Theta] = 0$, we can simultaneously label states by their parity eigenvalue $\xi = \pm 1$.

Step 2: At each TRIM, the occupied bands come in Kramers pairs $|u_{2m-1}\rangle$ and $|u_{2m}\rangle = \Theta|u_{2m-1}\rangle$. Since parity commutes with time reversal, both members of a Kramers pair share the same parity eigenvalue $\xi_{2m}(\Gamma_i)$.

Step 3: Fu and Kane showed that the sewing matrix Pfaffian ratio $\delta_i = \text{Pf}[w(\Gamma_i)]/\sqrt{\det[w(\Gamma_i)]}$ simplifies in the presence of inversion symmetry to:

$$\delta_i = \prod_{m=1}^{N} \xi_{2m}(\Gamma_i)$$

where the product is over all occupied Kramers pairs (one parity eigenvalue per pair).

Step 4: The Z$_2$ invariant is then:

$$(-1)^\nu = \prod_{i=1}^{4}\delta_i = \prod_{i=1}^{4}\prod_{m=1}^{N}\xi_{2m}(\Gamma_i)$$

In 2D there are 4 TRIM points; in 3D there are 8 TRIM points and the strong invariant $\nu_0$ uses the product over all 8.

Step 5 (Example: Bi$_2$Se$_3$): At the $\Gamma$ point, band inversion swaps the parity of the Se $p_z$ and Bi $p_z$ orbitals. The product of parities at $\Gamma$ gives $\delta_\Gamma = -1$, while the three equivalent L points each give $\delta_L = +1$. Therefore:

$$(-1)^\nu = (-1)\cdot(+1)^3 = -1 \quad \Longrightarrow \quad \nu = 1 \;\;\text{(topological)}$$

This is the simplest practical route to identifying topological insulators from ab initio band structure calculations: one only needs to compute parity eigenvalues at the high-symmetry points, a straightforward post-processing step in any DFT code.

Real-World Applications of Topological Insulators

1. Spintronics and Spin-Charge Conversion

The spin-momentum locking on TI surfaces generates a highly efficient spin-charge conversion mechanism. A charge current on the TI surface automatically produces a net spin polarization, and conversely, injecting a spin current generates a charge voltage. Bi$_2$Se$_3$ and related materials exhibit spin-charge conversion efficiencies (characterized by the inverse Edelstein length) up to 100 times larger than conventional heavy metals like Pt, making them promising for spin-orbit torque magnetic random access memory (SOT-MRAM).

2. Topological Quantum Computing

The interface between a 3D TI and an s-wave superconductor hosts Majorana zero modes in vortex cores (Fu-Kane proposal). These non-Abelian anyons could serve as the building blocks of topological qubits that are inherently protected against local sources of decoherence. Experimental signatures of proximity-induced superconductivity on TI surfaces have been observed in Bi$_2$Se$_3$/NbSe$_2$ and Bi$_2$Te$_3$/FeTe$_{0.55}$Se$_{0.45}$ heterostructures.

3. Thermoelectric Energy Conversion

Many topological insulator materials (Bi$_2$Te$_3$, Bi$_2$Se$_3$, Sb$_2$Te$_3$) are also the best-known thermoelectric materials, used in Peltier coolers and thermoelectric generators. The surface states contribute to the Seebeck coefficient and thermal conductivity, and the topological band structure favors the combination of high electrical conductivity and low thermal conductivity that maximizes the thermoelectric figure of merit ZT.

4. Infrared and Terahertz Detection

The gapless Dirac surface states of 3D TIs absorb light across a broad frequency range, including the technologically important terahertz gap. TI-based photodetectors have demonstrated room-temperature response from THz to mid-infrared frequencies, with applications in imaging, spectroscopy, and telecommunications. The high surface-to-volume ratio in TI nanostructures enhances the surface-state contribution to the photoresponse.

5. Catalysis and Topological Surface Chemistry

The metallic surface states of TIs provide a high density of states at the Fermi level on an otherwise insulating bulk, creating catalytically active surfaces. Bi$_2$Se$_3$ and Bi$_2$Te$_3$ have shown enhanced activity for hydrogen evolution and CO$_2$ reduction reactions. The spin texture of the surface states may also influence adsorption geometries and reaction pathways through spin-selective chemistry.

Historical Context

The theoretical prediction of topological insulators emerged from two parallel developments in 2005-2006. Charles Kane and Eugene Mele at the University of Pennsylvania proposed the quantum spin Hall effect in graphene with spin-orbit coupling (2005), introducing the Z$_2$ topological invariant that distinguishes the new phase from ordinary band insulators. Independently, B. Andrei Bernevig and Shou-Cheng Zhang at Stanford predicted (2006) that the quantum spin Hall effect would be observable in HgTe/CdTe quantum wells, where the band inversion between the $\Gamma_6$and $\Gamma_8$ bands provides the topological gap. This prediction was confirmed experimentally by Markus Konig et al. in Laurens Molenkamp's group at the University of Wurzburg in 2007, providing the first observation of the quantum spin Hall effect through quantized edge conductance$G = 2e^2/h$.

The extension to three dimensions was made by Liang Fu, Charles Kane, and Eugene Mele (2007), and independently by Joel Moore and Leon Balents (2007) and Rahul Roy (2009). Fu and Kane's parity criterion enabled rapid identification of candidate materials. The first 3D TI, the alloy Bi$_{1-x}$Sb$_x$, was confirmed by ARPES experiments by D. Hsieh et al. in M. Z. Hasan's group at Princeton in 2008. The "second generation" TIs -- Bi$_2$Se$_3$, Bi$_2$Te$_3$, and Sb$_2$Te$_3$ -- were predicted by H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang (2009) and quickly confirmed by ARPES. The broader classification of topological phases via the ten-fold way (Altland-Zirnbauer symmetry classes) was established by Andreas Schnyder, Shinsei Ryu, Akira Furusaki, and Andreas Ludwig (2008), and independently by Alexei Kitaev (2009). Duncan Haldane, David Thouless, and Michael Kosterlitz shared the 2016 Nobel Prize in Physics for their foundational work on topological phases of matter.

Python Simulation: Kane-Mele Band Structure

Bulk band structure of the Kane-Mele model showing the topological gap at the Dirac points, helical edge states in a zigzag nanoribbon, and the surface Dirac cone of a 3D topological insulator.

Kane-Mele Model: Bands, Edge States & Surface Dirac Cone

Python

Band structure with spin-orbit gap, ribbon edge states, and 3D TI surface dispersion

script.py125 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

# ──────────────────────────────────────────────────────────
# Kane-Mele Model: Band Structure & Edge States
# Graphene with intrinsic spin-orbit coupling
# ──────────────────────────────────────────────────────────

# Lattice and reciprocal vectors for honeycomb
a1, a2 = np.array([1.0, 0.0]), np.array([0.5, np.sqrt(3)/2])
b1 = 2*np.pi * np.array([1.0, -1.0/np.sqrt(3)])
b2 = 2*np.pi * np.array([0.0, 2.0/np.sqrt(3)])
Gamma, K, M = np.array([0.,0.]), (2*b1+b2)/3, (b1+b2)/2

# Nearest-neighbor and NNN vectors
d1, d2, d3 = np.array([0., 1/np.sqrt(3)]), np.array([.5, -.5/np.sqrt(3)]), np.array([-.5, -.5/np.sqrt(3)])
nnn = [a1, a2, a2-a1, -a1, -a2, a1-a2]
nu_signs = [1, 1, 1, -1, -1, -1]

t, t_so, t_R = 1.0, 0.06, 0.05

def kane_mele_H(kx, ky, sz=1):
    k = np.array([kx, ky])
    f_k = sum(np.exp(1j*np.dot(k, d)) for d in [d1, d2, d3])
    g_k = sum(nu*np.sin(np.dot(k, v)) for nu, v in zip(nu_signs, nnn))
    return np.array([[sz*t_so*g_k, -t*f_k], [-t*np.conj(f_k), -sz*t_so*g_k]], dtype=complex)

def k_path(points, npts=100):
    kpts, dists = [], [0.0]
    for i in range(len(points)-1):
        for j in range(npts):
            kpts.append(points[i] + j/npts*(points[i+1]-points[i]))
            if len(kpts) > 1: dists.append(dists[-1] + np.linalg.norm(kpts[-1]-kpts[-2]))
    kpts.append(points[-1])
    dists.append(dists[-1] + np.linalg.norm(kpts[-1]-kpts[-2]))
    return np.array(kpts), np.array(dists)

path_pts = [Gamma, K, M, Gamma]
kpts, kdist = k_path(path_pts, 120)

# Compute bulk bands for both spin sectors
bands_up, bands_dn = np.zeros((len(kpts), 2)), np.zeros((len(kpts), 2))
for i, kp in enumerate(kpts):
    bands_up[i] = np.sort(np.linalg.eigvalsh(kane_mele_H(kp[0], kp[1], +1)).real)
    bands_dn[i] = np.sort(np.linalg.eigvalsh(kane_mele_H(kp[0], kp[1], -1)).real)

# ── Zigzag ribbon for edge states ──
Ny, nk_ribbon = 40, 200
kx_vals = np.linspace(-np.pi, np.pi, nk_ribbon)
ribbon_bands = np.zeros((nk_ribbon, 4*Ny))

for ik, kx in enumerate(kx_vals):
    H = np.zeros((4*Ny, 4*Ny), dtype=complex)
    for iy in range(Ny):
        for sz_idx, sz in enumerate([+1, -1]):
            iA, iB = 4*iy + 2*sz_idx, 4*iy + 2*sz_idx + 1
            H[iA, iB] += -t*(1 + np.exp(-1j*kx))
            H[iB, iA] += -t*(1 + np.exp(1j*kx))
            H[iA, iA] += t_so*sz*2*np.sin(kx)
            H[iB, iB] -= t_so*sz*2*np.sin(kx)
            if iy < Ny - 1:
                jA, jB = 4*(iy+1)+2*sz_idx, 4*(iy+1)+2*sz_idx+1
                H[iB, jA] += -t; H[jA, iB] += -t
                H[iA, jA] += 1j*t_so*sz; H[jA, iA] -= 1j*t_so*sz
                H[iB, jB] -= 1j*t_so*sz; H[jB, iB] += 1j*t_so*sz
    ribbon_bands[ik] = np.sort(np.linalg.eigvalsh(H).real)

# ── Surface Dirac cone (3D TI model) ──
kx_s, ky_s = np.linspace(-0.5, 0.5, 100), np.linspace(-0.5, 0.5, 100)
KX, KY = np.meshgrid(kx_s, ky_s)
E_plus = np.sqrt(KX**2 + KY**2)
E_minus = -E_plus

# ── Plotting ──
fig = plt.figure(figsize=(18, 6))
fig.patch.set_facecolor('#0a0a1a')
gs = GridSpec(1, 3, figure=fig, width_ratios=[1, 1, 1.2])

for idx, ax_gs in enumerate(gs):
    if idx < 2:
        ax = fig.add_subplot(ax_gs)
        ax.set_facecolor('#0a0a1a')
        for sp in ax.spines.values(): sp.set_edgecolor('#334155')
        ax.tick_params(colors='#94a3b8')
        if idx == 0:
            ax.plot(kdist, bands_up[:,0], color='#a855f7', lw=2, label='Spin up')
            ax.plot(kdist, bands_up[:,1], color='#a855f7', lw=2)
            ax.plot(kdist, bands_dn[:,0], color='#06b6d4', lw=2, ls='--', label='Spin down')
            ax.plot(kdist, bands_dn[:,1], color='#06b6d4', lw=2, ls='--')
            tp = [kdist[0], kdist[120], kdist[240], kdist[-1]]
            ax.set_xticks(tp)
            ax.set_xticklabels(['$\\Gamma$','K','M','$\\Gamma$'], color='#94a3b8')
            for p in tp: ax.axvline(p, color='#334155', lw=0.5)
            ax.set_ylabel('E / t', color='#94a3b8'); ax.set_ylim(-3.5, 3.5)
            ax.set_title('Kane-Mele Bulk Bands', color='white', fontsize=14)
            ax.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
        else:
            n_mid = 2*Ny
            for n in range(4*Ny):
                c = '#a855f7' if abs(n-n_mid)<8 else '#334155'
                a = 1.0 if abs(n-n_mid)<8 else 0.3
                ax.plot(kx_vals, ribbon_bands[:,n], color=c, lw=1.5 if a==1 else 0.3, alpha=a)
            ax.set_xlabel('$k_x$', color='#94a3b8'); ax.set_ylabel('E / t', color='#94a3b8')
            ax.set_title('Zigzag Ribbon: Edge States', color='white', fontsize=14)
            ax.set_ylim(-0.5, 0.5); ax.axhline(0, color='#64748b', lw=0.5, ls='--')

ax3 = fig.add_subplot(gs[2], projection='3d')
ax3.set_facecolor('#0a0a1a')
ax3.plot_surface(KX, KY, E_plus, alpha=0.7, color='#a855f7', edgecolor='none')
ax3.plot_surface(KX, KY, E_minus, alpha=0.7, color='#06b6d4', edgecolor='none')
for lbl, axis in [('$k_x$', ax3.set_xlabel), ('$k_y$', ax3.set_ylabel), ('E', ax3.set_zlabel)]:
    axis(lbl, color='#94a3b8', fontsize=10, labelpad=5)
ax3.set_title('3D TI Surface Dirac Cone', color='white', fontsize=14, pad=10)
ax3.tick_params(colors='#94a3b8')
for pane in [ax3.xaxis.pane, ax3.yaxis.pane, ax3.zaxis.pane]:
    pane.fill = False; pane.set_edgecolor('#1a1a2e')

plt.tight_layout()
plt.savefig('kane_mele_bands.png', dpi=130, bbox_inches='tight', facecolor=fig.get_facecolor())
print(f"Kane-Mele: t={t}, t_SO={t_so}, t_R={t_R}, gap~{4*np.sqrt(3)*t_so:.4f}t")
print(f"Ribbon width: {Ny} cells. Edge states cross gap at TRIM points.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Simulation: Z₂ Invariant Computation

Computes the Z$_2$ topological invariant as a function of spin-orbit coupling strength, demonstrating the transition from trivial insulator to quantum spin Hall phase.

Z₂ Topological Invariant: Phase Transition

Fortran

Berry phase computation of Z₂ invariant across trivial-to-topological transition

program.f90105 lines

program z2_invariant
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: Nk = 60, Nlam = 40
  real(dp), parameter :: pi = 3.14159265358979323846_dp
  real(dp) :: t, t_so, lam_min, lam_max, dlam, kx, ky, dk
  real(dp) :: gap, min_gap, tr, det_val, E_plus, E_minus
  real(dp) :: f_re, f_im, g_val, norm_val, theta, total_berry
  integer :: i, j, il, il2, z2
  real(dp) :: Hre(2,2), Him(2,2)
  real(dp) :: d1x, d1y, d2x, d2y, d3x, d3y
  real(dp) :: nnn_x(6), nnn_y(6)
  integer  :: nu(6)
  real(dp) :: overlap_re, overlap_im
  real(dp) :: u_re(2), u_im(2), u_prev_re(2), u_prev_im(2)

d1x = 0.0_dp;  d1y = 1.0_dp / sqrt(3.0_dp)
  d2x = 0.5_dp;  d2y = -0.5_dp / sqrt(3.0_dp)
  d3x = -0.5_dp; d3y = -0.5_dp / sqrt(3.0_dp)
  nnn_x = (/1.0_dp, 0.5_dp, -0.5_dp, -1.0_dp, -0.5_dp, 0.5_dp/)
  nnn_y = (/0.0_dp, sqrt(3.0_dp)/2, sqrt(3.0_dp)/2, &
            0.0_dp, -sqrt(3.0_dp)/2, -sqrt(3.0_dp)/2/)
  nu = (/1, 1, 1, -1, -1, -1/)

t = 1.0_dp; lam_min = 0.0_dp; lam_max = 0.30_dp
  dlam = (lam_max - lam_min) / real(Nlam - 1, dp)
  dk = 2.0_dp * pi / real(Nk, dp)

write(*,'(A)') "======================================================"
  write(*,'(A)') "   Z2 Topological Invariant: Kane-Mele Model"
  write(*,'(A)') "   Transition between trivial and topological phases"
  write(*,'(A)') "======================================================"
  write(*,'(A)') ""
  write(*,'(A)') "  lambda_SO     min_gap     Z2    Phase"
  write(*,'(A)') "  ---------    ---------   ----   -----"

do il = 1, Nlam
    t_so = lam_min + real(il - 1, dp) * dlam
    min_gap = 1.0d10; total_berry = 0.0_dp

do j = 0, Nk
      ky = real(j, dp) * dk
      do i = 0, Nk
        kx = real(i, dp) * dk
        f_re = cos(kx*d1x+ky*d1y) + cos(kx*d2x+ky*d2y) + cos(kx*d3x+ky*d3y)
        f_im = sin(kx*d1x+ky*d1y) + sin(kx*d2x+ky*d2y) + sin(kx*d3x+ky*d3y)
        g_val = 0.0_dp
        do il2 = 1, 6
          g_val = g_val + real(nu(il2),dp)*sin(kx*nnn_x(il2)+ky*nnn_y(il2))
        end do

Hre(1,1) = t_so*g_val;  Him(1,1) = 0.0_dp
        Hre(1,2) = -t*f_re;     Him(1,2) = -t*f_im
        Hre(2,1) = -t*f_re;     Him(2,1) = t*f_im
        Hre(2,2) = -t_so*g_val; Him(2,2) = 0.0_dp

tr = Hre(1,1) + Hre(2,2)
        det_val = Hre(1,1)*Hre(2,2) - Hre(1,2)*Hre(2,1) - Him(1,2)*Him(2,1)
        E_plus = 0.5_dp*(tr + sqrt(max(tr*tr - 4.0_dp*det_val, 0.0_dp)))
        E_minus = 0.5_dp*(tr - sqrt(max(tr*tr - 4.0_dp*det_val, 0.0_dp)))
        gap = E_plus - E_minus
        if (gap < min_gap) min_gap = gap

norm_val = sqrt((Hre(1,1)-E_minus)**2 + Hre(1,2)**2 + Him(1,2)**2)
        if (norm_val > 1.0d-10) then
          u_re(1) = -Hre(1,2)/norm_val; u_im(1) = -Him(1,2)/norm_val
          u_re(2) = (Hre(1,1)-E_minus)/norm_val; u_im(2) = 0.0_dp
        else
          u_re = (/1.0_dp, 0.0_dp/); u_im = (/0.0_dp, 0.0_dp/)
        end if

if (i > 0) then
          overlap_re = u_prev_re(1)*u_re(1) + u_prev_im(1)*u_im(1) &
                     + u_prev_re(2)*u_re(2) + u_prev_im(2)*u_im(2)
          overlap_im = u_prev_re(1)*u_im(1) - u_prev_im(1)*u_re(1) &
                     + u_prev_re(2)*u_im(2) - u_prev_im(2)*u_re(2)
          total_berry = total_berry + atan2(overlap_im, overlap_re)
        end if
        u_prev_re = u_re; u_prev_im = u_im
      end do
    end do

z2 = mod(nint(total_berry / pi), 2)
    if (z2 < 0) z2 = z2 + 2

if (t_so < 1.0d-6) then
      write(*,'(F10.4, F12.6, I7, A)') t_so, min_gap, 0, "   Trivial (graphene)"
    else if (min_gap < 0.01_dp) then
      write(*,'(F10.4, F12.6, I7, A)') t_so, min_gap, z2, "   Gap closing"
    else if (z2 == 1) then
      write(*,'(F10.4, F12.6, I7, A)') t_so, min_gap, z2, "   Topological (QSH)"
    else
      write(*,'(F10.4, F12.6, I7, A)') t_so, min_gap, z2, "   Trivial"
    end if
  end do

write(*,'(A)') ""
  write(*,'(A)') "======================================================"
  write(*,'(A)') "  Summary:"
  write(*,'(A)') "  - lambda_SO = 0: gapless graphene (trivial)"
  write(*,'(A)') "  - lambda_SO > 0: topological gap opens"
  write(*,'(A)') "  - Z2 = 1 indicates quantum spin Hall insulator"
  write(*,'(A)') "  - Protected helical edge states at boundaries"
  write(*,'(A)') "======================================================"
end program z2_invariant

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Chapter 2: Quantum Hall Effect Course Overview →