Hofstadter (Chern butterfly) — C++ port of Hofstadter.ipynb

Chern-number butterfly — q_max=180, sep=500, n_mu=8000 (9879×8000, out/Chern_numbers.png).

A self-contained C++17 + OpenMP rewrite of the notebook computation:

1. For every coprime pair (p,q), 1 <= p < q <= q_max (sorted by phi = p/q), builds the q×q complex Hermitian Hofstadter Hamiltonian

   H[m][m]        = 2 cos(ky + 2π·m·p/q)
   H[m][(m+1)%q]  = e^{-i kx},   H[(m+1)%q][m] = e^{+i kx}
   

   on the (sep × sep) Brillouin-zone grid kx ∈ [-π,π], ky ∈ [-π/q, π/q] and extracts its eigenvalues, keeping the band edges bandtop[i] / bandbottom[i] per band index. The default solver is a structured cyclic-tridiagonal secular solver (see Speed below); --dense selects Eigen's SelfAdjointEigenSolver reference path.

2. Evaluates the Chern number of every insulating gap on an n_mu-point chemical-potential grid, with the exact same mask logic and the modular (TKNN Diophantine) solution σ ≡ p⁻¹·r (mod q) of the notebook; 999 marks "not in a gap" (= NaN in the notebook plot).

3. Writes out/pairs.csv, out/chern.csv, out/edges.txt, out/meta.txt and renders the butterfly figure out/Chern_numbers.png (matplotlib's coolwarm_r colormap with log(1+|C|) normalisation, gray for NaN, LUT baked into include/colormap.h).

Layout

hofstadter.cpp                 single translation unit (physics + CLI + PNG)
include/colormap.h             matplotlib coolwarm_r LUT
third_party/stb_image_write.h  vendored single-header PNG writer
reference/hofstadter_ref.py    notebook reference implementation (validation)
Hofstadter.ipynb               original notebook
CMakeLists.txt                 build + Eigen fetch
out/                           generated results (gitignored)

Build & run

CMake fetches Eigen from GitLab on the first configure (tag HOF_EIGEN_TAG, default 3.4.1); a C++17 compiler and OpenMP are the only other requirements.

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j

# notebook parameters: q_max=40 sep=100 n_mu=1000
./build/hofstadter --q_max 40 --sep 100 --n_mu 1000 \
    --outdir out --edges out/edges.txt

./build/hofstadter --help
# useful flags: --q_max N --sep N --n_mu N --threads N --outdir D
#               --edges F  --png P --scale S  --no-png  --dense

Useful CMake options:

-DHOF_NATIVE=OFF        disable -march=native
-DHOF_PROFILE=ON        enable the HOF_PROFILE timers (writes out/prof.csv)
-DHOF_EIGEN_TAG=3.4.0   fetch a different Eigen tag

The build/ directory is generated by CMake and is gitignored; CMake embeds absolute paths in its generated makefiles/CMakeCache there. The source tree itself is path-independent, so delete build/ and reconfigure whenever a clean, relocatable tree is needed.

Speed

Hofstadter.ipynb (Pool(12), numpy/LAPACK, full 100×100 k-grid):

q_max=40  sep=100  n_mu=1000   →  ≈ 28 s   (12 workers, this machine)

This C++ version, same parameters (./build/hofstadter --q_max 40 --sep 100 --n_mu 1000):

dense reference (--dense)   →  ≈ 2.1 s   (20 threads)
structured (default)        →  ≈ 0.78 s  (20 threads)

i.e. ≈ 36× faster than the notebook and ≈ 6.4× faster than the straightforward dense C++ port, with the Chern matrix bit-identical (see validation below). --dense selects the old dense-eigensolver path; the two paths produce the same Chern matrix (band edges agree to ≈ 3e-14).

Why it is faster (all reductions are exact, no change to the discrete k-grid):

• 4× fewer diagonalisations (four-fold k symmetry). The spectrum is invariant under kx → -kx (H(-kx,ky) = conj H(kx,ky)) and under the Bloch time-reversal relation H(-kx,-ky) = H(kx,ky)ᵀ, so spectrum(kx,ky) = spectrum(±kx, ±ky); the notebook solves all 100×100 = 10⁴ k-points per (p,q), this code solves only one representative of every orbit {(i,j), (s-i,j), (i,s-j), (s-i,s-j)} with s = sep-1 (50×50 = 2500 points, the lower-left quadrant) and takes the same per-index max/min. Eigenvalues are sorted ascending at every k, so per-band edges are unchanged to machine precision.
• p ↔ q−p symmetry (2×). H^(q-p)(kx,ky) = H^(p)(kx,-ky) and the ky grid is symmetric, so the two pairs have identical band edges; only p < q/2 pairs are solved and the edges are copied to their mirrors. Chern numbers are still evaluated with each pair's own p.
• Magnetic Brillouin-zone periodicity. The spectrum is periodic in kx with period 2π/q and even in kx, so the 50 lower-left kx representatives collapse to ceil(P/2) distinct spectra, P = (sep-1)/gcd(q,sep-1); the omitted points have identical spectra and are simply not solved.
• Structured cyclic-tridiagonal eigensolver (the big one). The q×q matrix is bordered by its leading (q−1)×(q−1) open chain B. With B = Q diag(μ) Qᵀ and w_i = γ Q_{i,0} + Q_{i,q-2} (γ = e^{-i q kx}), the eigenvalues are the roots of the scalar secular equation (a_{q-1}−λ) − Σ |w_i|²/(μ_i−λ) = 0, one per interval between consecutive μ_i. B is diagonalised once per ky and reused for every kx; each kx then costs O(q²) (bracketed Newton on a monotone function) instead of the dense O(q³) complex Householder + QR. This is ~2–3.5× on the eigensolver and, combined with the symmetries above, gives the overall ~6.4×.
• Parallelised over the (p,q) pairs with OpenMP schedule(dynamic); heavy pairs (large q) are spread evenly.

A profiling pass (-DHOF_PROFILE=ON) shows the time is now split between the per-ky open-chain eigendecomposition and the per-kx secular root-finding; matrix setup, Chern counting and I/O are negligible.

Scaling

Measured on this machine (sep=100, n_mu=1000, 20 threads, --no-png):

q_max   pairs   structured
 40      489      0.79 s
 60     1101      3.91 s
 80     1965     13.16 s
100     3043     31.60 s
120     4385     70.22 s

The runtime grows roughly as q_max^4.1–4.4: the number of pairs grows as q_max², and the per-pair work is dominated by the O(q³) per-ky open-chain eigendecomposition plus the O(q²) per-kx secular root-finding. Extrapolating, q_max=200 (12,231 pairs) is expected to take ≈ 8–12 minutes for the computation (--no-png), possibly a little more if the O(q³) term dominates at large q. Without --no-png the default --scale 2 would render a 24,462 × 2,000 PNG (~147 MB raw), so use --no-png or --scale 1 for large q_max.

Validation

reference/hofstadter_ref.py is the notebook code, run once with 12 workers (q_max=40, sep=100, n_mu=1000) to produce a reference pickle. Against it the original dense C++ path gave:

band edges (all 489 pairs × q bands):  max |Δ| ≤ 2.5e-14
Chern matrix  (489 × 1000):             0 / 489000 cells differ
NaN (999) pattern:                      identical

The structured path (default) is checked against the dense path (--dense) on the same parameters:

Chern matrix:   byte-identical (0 / 489000 cells differ)
band edges:     max |Δ| ≤ 3.1e-14
rendered PNG:   identical md5

It is also checked for odd/even sep and several q_max values (Chern matrix identical, band edges ≤ 1.5e-14 between the two paths).

To rerun the reference check:

python3 -m venv .venv && .venv/bin/pip install numpy scipy
.venv/bin/python reference/hofstadter_ref.py \
    --q_max 40 --sep 100 --n_mu 1000 --workers 12 --out reference/ref_full.pkl
# then compare out/pairs.csv, out/chern.csv and out/edges.txt
# against the pickle (or configure with -DHOF_PROFILE=ON and inspect
# out/prof.csv for the per-pair solver breakdown)

Notes

• Python 2-arg % and the notebook's mod_inv are reproduced literally; values involved are non-negative so C++ and Python semantics coincide.
• q_max=40 gives 489 coprime pairs (not 499); q_max=200 gives 12,231.
• The notebook docstring for coprime_pairs ("q is a factor of n") does not match its code (all coprime pairs are enumerated); the code is followed.