Notes on the interactive figures

Quarto interactive figures are powered by Observable JavaScript (OJS) code cells embedded directly in each chapter source. This appendix opens with a short tour of the technologies used across the book, then drills into the algorithms behind a handful of figures worth a deeper look.

Technology stack for OJS figures

Every interactive figure follows the same scaffolding:

//| label: fig-something
import { renderSomething } from "./_interactive/chNN-something.js"
const W = Math.min(width, 860);
const container = document.createElement("div");
container.style.cssText =
  "width:100%; max-width:" + W + "px; margin:0 auto;";
yield container;
renderSomething(container);

The OJS cell only sets up the container and yields it to the page — all the real work happens inside an ES-module under _interactive/chNN-…js. The libraries used across the book:

Three.js (https://cdn.jsdelivr.net/npm/three@…) — the WebGL workhorse. Used for every 3D figure that is not an atomic structure: Brillouin-zone polyhedra, Fermi surfaces, marching-cubes iso-surfaces. We rely on OrbitControls for camera, on localClippingEnabled + THREE.Plane for per-fragment clipping, and on MeshPhysicalMaterial for nice physically-based shading with metalness, roughness and clearcoat. Resize is handled by a ResizeObserver on the container that updates the renderer size and camera aspect, coalesced with requestAnimationFrame so resize bursts don’t thrash the GPU.
3Dmol.js (loaded once globally by _header.html as the $3Dmol global) — used whenever the figure shows actual atoms and bonds (Bravais lattices, crystal structures, defects, COD entries). Custom geometries (transparent unit-cell boxes, symmetry-axis arrows, plane mirrors) are built with viewer.addCustom({vertexArr, normalArr, faceArr, color, opacity}). Two project-specific conventions are documented in the main CLAUDE.md: the atoms- visible checkbox uses sphere.scale = 0.01 (not 0, which would also hide the bonds), and double-sided transparent surfaces duplicate the vertex list with reversed winding so that both faces shade correctly.
Canvas 2D — the default for plots, dispersion curves, DOS, BZ cross-sections, schematics, etc. Drawn with the standard getContext("2d") API; sizing follows the responsive pattern (canvas.width = container.clientWidth · DPR, style.width = container.clientWidth + "px") and re-paint is triggered by the same ResizeObserver + RAF combination used for Three.js.
SVG — used sparingly, only where the figure benefits from scalable vector text or hyperlinks (notably some of the 2D Bragg / Miller / wallpaper diagrams).
KaTeX_Math font — loaded by Quarto’s HTML format. Every label containing a mathematical symbol (variables in italic, Greek, subscripts, units like Å or eV) is drawn in '"KaTeX_Math","Times New Roman",serif', while UI labels (button text, dropdown options) stay in sans-serif. A small in-canvas helper, drawLabel(ctx, segments, x, y, size, color), takes a list of {t, i?, sub?} segments and assembles “italic variable + roman subscript + roman units” strings the way a real LaTeX paper would.
Reactive widgets — when a figure exposes parameters to the reader they live either inside the figure’s own DOM (HTML <input type="range">, <select>, custom segmented buttons) or as Inputs.range / Inputs.select cells in the OJS block above. The in-figure DOM widgets call back into the figure’s draw() directly (and use requestAnimationFrame coalescing for high-frequency drag events such as slider input).
MathJax — handles all body-text equations. Quarto switches it on by default; the figure layer never touches it.
Animation primitives — most “the curves slide between two configurations” effects (e.g. extended ↔︎ reduced zone toggle, progressive filling of the n-th BZ, layer-by-layer reveal in close-packing, camera fly-throughs) are implemented with a requestAnimationFrame loop that drives a single scalar parameter $t \in [0, 1]$ from a Start value to a Target value. Each frame, $t$ is updated from the elapsed time, smoothed with a cubic ease-in/out
```
const easeInOut = x => x < 0.5
  ? 2*x*x
  : 1 - Math.pow(-2*x + 2, 2) / 2;
```
and the figure’s geometry is rebuilt by linear interpolation between the two endpoint states. For 1D / 2D plot states the interpolation is a plain lerp(a, b, t); for camera quaternions in 3D figures it is a normalised LERP (NLERP) — linear blend followed by renormalisation, cheaper than spherical SLERP and visually indistinguishable for the small angle gaps used here. The same RAF gates re-paints to the screen refresh rate (no unnecessary compute when the animation has settled or the figure is off-screen — an IntersectionObserver pauses the loop in that case).
StatiCrypt — only chapters 4–10 are protected behind a password (the unfinished portion of the book). Chapters 1–3 and this appendix are served plain so the open textbook can be linked to freely.

Fig. 1.14 — Periodic table of elemental crystals

Source: _interactive/ch01-periodic-table.js.

A 2D HTML grid of the periodic table; clicking on an element opens a 3Dmol viewer of its room-temperature crystal structure. The element data is a static table baked into the file; no on-the-fly fetching. The structure is one of BCC, FCC, HCP, DIA(diamond) (or null for elements that crystallise in a more complex / molecular form not covered by the simple primitive-cell viewer).

Lattice parameters are in Å; for HCP, both a and c are given. The full table currently shipped:

Element	Z	Structure	a (Å)	c (Å)	Notes
H	1	HCP	3.75	6.12	T < 14 K
He	2	HCP	3.57	5.83	T < 4 K, P > 25 bar
Li	3	BCC	3.51	—
Be	4	HCP	2.29	3.58
C	6	DIA	3.57	—	metastable (graphite is stable)
Ne	10	FCC	4.43	—	T < 25 K
Na	11	BCC	4.29	—
Mg	12	HCP	3.21	5.21
Al	13	FCC	4.05	—
Si	14	DIA	5.43	—
Ar	18	FCC	5.26	—	T < 84 K
K	19	BCC	5.23	—
Ca	20	FCC	5.58	—
Sc	21	HCP	3.31	5.27
Ti	22	HCP	2.95	4.69
V	23	BCC	3.02	—
Cr	24	BCC	2.88	—
Fe	26	BCC	2.87	—	(α-Fe, room T)
Co	27	HCP	2.51	4.07
Ni	28	FCC	3.52	—
Cu	29	FCC	3.61	—
Zn	30	HCP	2.66	4.95
Ge	32	DIA	5.66	—
Kr	36	FCC	5.72	—	T < 116 K
Rb	37	BCC	5.59	—
Sr	38	FCC	6.08	—
Mo	42	BCC	3.15	—
Rh	45	FCC	3.80	—
Pd	46	FCC	3.89	—
Ag	47	FCC	4.09	—
Sn	50	DIA	6.49	—	T < 13 °C (α-Sn)
Xe	54	FCC	6.20	—	T < 161 K
Cs	55	BCC	6.05	—
Ba	56	BCC	5.02	—
W	74	BCC	3.16	—
Pt	78	FCC	3.92	—
Au	79	FCC	4.08	—
Pb	82	FCC	4.95	—

(Lanthanides and actinides — many of which are DHCP, double-HCP — are also in the file; the table above is a representative subset. For the full list see the source file.)

When the user clicks on an element with a non-null struct, the figure spawns a fresh $3Dmol.createViewer and adds atoms at the positions implied by the structure and the element’s a, c parameters. Only structural geometry is shown — no bonds beyond nearest-neighbour visual sticks.

Fig. 1.10 — Close-packing (HCP vs FCC)

Source: _interactive/ch01-close-packing.js.

A 3D-rendered demonstration of how the two close-packed structures emerge from the same triangular layer: starting from a single hexagonal layer of touching spheres, the figure animates the addition of a second layer in the hollows, then offers two choices for the third layer (above the first → ABA stacking → HCP, or above the other hollows → ABC → FCC).

The interesting bit on the implementation side is the camera choreography. As the user toggles between layer 1, layer 2, HCP, FCC, the camera doesn’t snap — it interpolates between two viewpoints (each viewpoint = a position vector + an up vector + a target vector, internally encoded as a quaternion + distance). Each viewpoint slot has 8 numbers (3 position + 3 target + 2 quaternion components in our reduced parametrisation, with sign- disambiguation on the second half so the quaternion takes the short path):

function lerp(a, b, t) { return a + (b - a) * t; }
function step(s, e, raw) {
  const t = easeInOut(raw);
  const v = s.map((si, i) => lerp(si, e[i], t));
  /* sign-disambiguate quaternion components 4..7 for the short arc */
  if (s[4]*e[4] + s[5]*e[5] + s[6]*e[6] + s[7]*e[7] < 0) {
    for (let k = 4; k <= 7; k++) v[k] = lerp(s[k], -e[k], t);
  }
  applyView(v);
  if (raw < 1) animId = requestAnimationFrame(step);
}

The “appearance” of the new layer of spheres is itself animated: each sphere’s vertical coordinate is interpolated from a “high” parking position down to its destination over the same $t$ window, so the two effects (camera glide + layer drop) read as a single fluid motion. Spheres / sticks themselves are 3Dmol.js, but the interpolation logic and the RAF loop live entirely in the figure’s JS file.

Fig. 3.1 — Reduced/extended/repeated zone schemes

Source: _interactive/ch03-zone-schemes.js.

What is “calculated” here is in fact the simplest possible band: the 1D free-electron parabola

\[ E(k) = \frac{\hbar^2 k^2}{2m}. \]

The figure draws three equivalent representations of this same single-particle dispersion:

Extended zone: the bare parabola, plotted over a wide $k$ range that crosses several Brillouin zones.
Reduced zone: the parabola is folded back into the first BZ $[-\pi/a, \pi/a]$ by repeatedly applying $k \to k \pm 2\pi/a$ and re-labelling each piece as a separate band $E_n(k)$.
Repeated zone: the reduced-zone bands are then translated by every reciprocal lattice vector to fill the entire $k$-axis with periodic copies — emphasising that band energies are periodic in $k$ with period $2\pi/a$.

There is no eigensolver involved. The “model” is the trivial empty lattice (zero potential), and the only computation is the integer shift $k_{\text{red}} = k_{\text{ext}} - n \cdot 2\pi/a$ for the folding.

The animation. Toggling between Extended and Reduced drives a parameter $t_\text{anim} \in [0, 1]$ over ANIM_MS = 1000 ms. Each band segment $n$ is plotted as $E_n(k - n \cdot 2\pi/a \cdot t_\text{anim})$: at $t = 1$ the segments sit on the bare parabola (extended), at $t = 0$ they have all collapsed back into the first BZ (reduced). The cubic ease-in/ out described in the preamble shapes the motion. The visual effect is that the high-energy parabola “pieces” fly inward and stack above the first BZ — exactly the geometric construction of band-folding worked out in the chapter text.

References: Ashcroft & Mermin, ch. 8; Kittel, ch. 7.

Fig. 3.2 — Successive Brillouin zones (2D square lattice)

Source: _interactive/ch03-brillouin-zones.js.

A 2D square reciprocal lattice with a discrete control that progressively highlights the first $n$ Brillouin zones, $n = 1 \ldots 7$. The $n$-th BZ is the locus of $\mathbf{k}$-points whose Voronoi rank in the reciprocal lattice is exactly $n$: the point has $n - 1$ reciprocal vectors strictly closer to it than the origin. The classification is done on a fine pixel raster.

The animation. Stepping from $n$ to $n + 1$ doesn’t simply add the new region — it animates the folding of the higher zones into the first by a sequence of integer reciprocal-lattice translations, one group at a time. Each group fold takes ANIM_MS_PER_GROUP = 1000 ms (with a minimum total animation of ANIM_MS_MIN = 900 ms), with the same easeInOut shape used elsewhere in the book. The result is that the user sees each fragment of the $n$-th zone slide across the lattice and snap into its reduced-zone position inside the first BZ — making explicit the fact that every $n$-th BZ has the same total area as the first, just chopped up and dispersed.

The same primitive serves a second purpose: it is the visual analogue of the band-folding animation in Fig. 3.1, generalised to 2D. Watching them both makes the otherwise abstract notion of “what the higher Brillouin zones really look like” concrete.

Fig. 3.4 — Split-step wavepacket propagation

Source: _interactive/ch03-splitstep.js.

Solves the time-dependent Schrödinger equation in 1D,

\[ i\hbar \,\partial_t \psi(x,t) = \left[ T + V(x) \right] \psi(x,t), \qquad T = -\frac{\hbar^2}{2m}\,\partial_x^2, \]

for a Gaussian wavepacket scattering off a user-chosen potential landscape (single barrier, double-barrier resonance, periodic chain of barriers).

The integrator is the split-step Fourier method (also known as the split-operator or Feit–Fleck–Steiger method, originally introduced for nonlinear optics and quickly adopted for QM). The exact propagator over a small time step $\Delta t$ is approximated as

\[ e^{-i (T+V) \Delta t / \hbar} \approx e^{-i V \Delta t / 2\hbar}\, e^{-i T \Delta t / \hbar}\, e^{-i V \Delta t / 2\hbar} \;+\; \mathcal{O}(\Delta t^3), \]

(symmetric Strang splitting). Each factor is diagonal in either real or momentum space:

$e^{-i V \Delta t / 2\hbar}$ is just pointwise multiplication on the spatial grid;
$e^{-i T \Delta t / \hbar}$ is pointwise multiplication on the $k$-grid, where $T = \hbar^2 k^2 / 2m$ is diagonal.

A Fast Fourier Transform moves between the two representations. One time step is therefore: half-kick in real space, FFT, full drift in momentum space, inverse FFT, second half-kick:

function step(psi, V, dt) {
  for (let i = 0; i < N; i++) psi[i] = mul(psi[i], expV[i]);   /* V/2 */
  fft(psi, +1);
  for (let i = 0; i < N; i++) psi[i] = mul(psi[i], expT[i]);   /* T   */
  fft(psi, -1);
  for (let i = 0; i < N; i++) psi[i] = mul(psi[i], expV[i]);   /* V/2 */
}

The two phase tables expV[i] = exp(-i V[i] dt/2 ℏ) and expT[i] = exp(-i ℏ k[i]² dt / 2m) are pre-computed once. The FFT is a hand-rolled radix-2 Cooley–Tukey on Float64Array buffers (N = 2048 for the figure), with bit-reversal index permutation done in place. All working arrays — psi, Vout, the two phase tables — are allocated once at start-up and overwritten in place: this kills all per-step GC, which had been responsible for visible RAF stalls on tablet hardware.

References: Numerical Recipes (3rd ed.), §13.1 (FFT) and §20.3.4 (“Split-operator method for the Schrödinger equation”); Feit, Fleck, Steiger, Solution of the Schrödinger equation by a spectral method, J. Comput. Phys. 47, 412 (1982).

Fig. 3.5 — N-barrier transmission via transfer matrices

Source: _interactive/ch03-n-barriers.js.

Computes the transmission probability $T(E)$ for a 1D chain of $N$ identical rectangular barriers (height $V_0$, width $b$, gap $g$) using the transfer matrix method.

For a single rectangular barrier sitting in $[0, b]$, the wavefunction on the left and right is a superposition of two plane waves:

\[ \psi_{L}(x) = A_L e^{ikx} + B_L e^{-ikx}, \qquad \psi_{R}(x) = A_R e^{ikx} + B_R e^{-ikx}, \]

with $k = \sqrt{2 m E}/\hbar$ outside and $q = \sqrt{2 m (V_0 - E)}/\hbar$ (real or imaginary depending on $E \lessgtr V_0$) inside. Continuity of $\psi$ and $\psi'$ at the two interfaces $x=0$ and $x=b$ closes the linear system, yielding a $2\times2$ matrix $M_b(E)$ with

\[ \begin{pmatrix} A_R \\ B_R \end{pmatrix} = M_b(E) \begin{pmatrix} A_L \\ B_L \end{pmatrix}. \]

The closed-form entries of $M_b$ involve $\cosh / \sinh$ for $E < V_0$ (tunnelling regime) and $\cos / \sin$ for $E > V_0$ (above-barrier). A free region of width $g$ is also a $2\times2$ matrix $M_g(E) = \mathrm{diag}(e^{ikg}, e^{-ikg})$.

For $N$ barriers the total transfer matrix is just the product:

\[ M_{\mathrm{tot}}(E) = M_b \, M_g \, M_b \, M_g \, \cdots \, M_b \]

with $N$ copies of $M_b$ alternated with $N-1$ copies of $M_g$. Once $M_{\mathrm{tot}} = \begin{pmatrix} m_{11} & m_{12} \\ m_{21} & m_{22}\end{pmatrix}$ is known, the transmission probability is

\[ T(E) = \frac{1}{|m_{22}|^2}. \]

(With our boundary conditions: incident (A_L, 0), transmitted (A_R, 0), reflected only on the left.)

The figure scans $E$ over a grid of $\sim 1000$ values and plots $T(E)$. Two qualitative regimes are visible:

$N = 1$: smooth transmission, with tunnelling for $E < V_0$.
$N = 2$: sharp resonances at energies that match the quasi-bound levels of the inter-barrier well (perfect transmission $T = 1$ at the resonance peaks).
$N \to \infty$: the discrete resonance peaks broaden into continuous transmission bands separated by gaps — the band-formation argument the chapter is about.

The energy grid is a pre-allocated Float64Array (re-used across re-evaluations to avoid GC pressure on tablet).

References: Numerical Recipes §11.2 (matrix products); Cohen- Tannoudji, Diu, Laloë, Quantum Mechanics vol. 1, complement III; Ashcroft & Mermin §8 for the band-formation interpretation.

Figs. 4.6 / 4.7 / 4.8 — Density of states (1D, 2D, 3D)

Sources: _interactive/ch03-1d-dos.js, _interactive/ch03-squarium-dos.js, _interactive/ch03-cubium-dos.js.

These three figures share a common goal — computing the density of states (DOS) $D(E)$ for a tight-binding band — but use three different algorithms, scaling with the dimensionality.

1D — closed form (Fig. 4.6)

The dispersion is $E(k) = -2 t \cos(k a)$ on the interval $[-\pi/a, \pi/a]$. The DOS is

\[ D(E) = \frac{1}{\pi}\,\left|\frac{dk}{dE}\right| = \frac{1}{\pi\sqrt{4t^2 - E^2}}, \]

an explicit closed form with a square-root van Hove singularity at each band edge. The figure simply evaluates this formula on a grid of $E$. No numerics beyond Math.sqrt.

2D square — elliptic integral (Fig. 4.7)

The dispersion is $E(\mathbf{k}) = -2 t (\cos k_x a + \cos k_y a)$ over $[-\pi/a, \pi/a]^2$. The DOS, integrated over the BZ, gives

\[ D(E) = \frac{1}{2\pi^2 t}\, K\!\left(\sqrt{1 - (E/4t)^2}\right), \]

where $K(m)$ is the complete elliptic integral of the first kind:

\[ K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - m \sin^2 \theta}}. \]

The signature logarithmic van Hove singularity at $E = 0$ (the saddle point of the dispersion) is the divergence of $K(m)$ as $m \to 1$. We evaluate $K(m)$ with the arithmetic-geometric mean (AGM) recurrence, which converges quadratically — six iterations are enough to reach machine precision:

function ellipticK(m) {
  let a = 1, b = Math.sqrt(1 - m);
  for (let i = 0; i < 8; i++) {
    const t = (a + b) / 2; b = Math.sqrt(a * b); a = t;
  }
  return Math.PI / (2 * a);
}

Reference: Numerical Recipes §6.11 (“Elliptic integrals and Jacobi elliptic functions”), eq. (6.11.3) for AGM.

3D cubic — marching cubes (Fig. 4.8)

The dispersion is $E(\mathbf{k}) = -2 t (\cos k_x + \cos k_y + \cos k_z) a$ over the cube. There is no clean closed form for the 3D DOS of this dispersion, so we evaluate it by histogramming: sample $E(\mathbf{k})$ on a regular grid of $N^3$ k-points and bin the results into a fine $E$-grid; the histogram approximates $D(E)$ up to the bin width.

The Fermi iso-surface (left panel, 3D) is then extracted with Marching Cubes (Lorensen & Cline 1987): for each voxel of the sampling grid, look up the topology of the implicit surface $E(\mathbf{k}) = E_F$ from a 256-entry table and emit triangles that interpolate the surface across voxel edges.

We exploit the cubium’s $O_h$ point-group symmetry (48 operations) to evaluate the dispersion only on the canonical wedge $k_x \ge k_y \ge k_z \ge 0$ (factor $\sim 48$ saving), then splat the results onto the up-to-6 wedge permutations and finally mirror the resulting mesh through the three coordinate planes to fill the full BZ. The truncated-cube clipping is done per-fragment by THREE.Plane clipping planes attached to the material.

References: Lorensen & Cline, Marching Cubes, Computer Graphics 21:163, 1987; Chernyaev, Marching Cubes 33: Construction of Topologically Correct Isosurfaces (1995, the variant implemented by the npm isosurface package we use).

Fig. 8.2 — Quantum Bloch oscillator (real-space split-step Fourier)

Source: _interactive/ch04-bloch-oscillator.js.

The figure is a genuine real-space solution of the time-dependent Schrödinger equation for a wave packet in a tilted periodic potential, integrated with the spectrally-accurate split-step Fourier method (the same propagator as the resonant-tunnelling figure, Fig. 3.4 — Split-step wavepacket propagation) and seeded with a Wannier–Stark coherent state — the band analogue of a harmonic-oscillator coherent state. The top panel shows the packet with the same conventions as the resonant-tunnelling figure ($|\psi|$ yellow envelope band, Re $\psi$ blue, Im $\psi$ red) inside the washboard potential; the bottom panel shows the folded quasimomentum density $\sum_G|\psi(k+G)|^2$ in the repeated-zone scheme (solid in the FBZ, fainter $\pm 2\pi/a$ replicas outside), drifting and Bragg-reflecting at the zone boundaries over the three lowest bands $E(k)$.

Physical model

In natural units $\hbar = m = a = 1$, the Hamiltonian is

\[ \mathcal{H} = -\frac{1}{2}\,\partial_x^2 + V(x), \qquad V(x) = -V_0\cos\!\left(\frac{2\pi x}{a}\right) - F x, \]

a cosine lattice of depth $V_0$ tilted by the uniform field $F = eEa$. Three parameters are exposed: $V_0$, $F$ and the integer Wannier–Stark rung $n$ that sets the launch site.

Algorithm

The figure runs in three stages.

Wannier–Stark coherent packet. The Wannier–Stark ladder is exactly equispaced ($E_n = \bar E_1 - nFa$, $\hbar\omega_B = Fa$), like the harmonic-oscillator spectrum — so a Gaussian superposition of ladder rungs plays the role of a coherent state: its rung weights are stationary and it revives every Bloch period instead of spreading. The state is synthesised directly in the band-1 Bloch basis on the exact Born–von-Kármán $k$-grid of the box ($N_\text{cells} = L/a = 90$ bins),

\[ \psi(x,0) \;=\; \sum_{k}\, g(k-k_0)\; e^{i\theta(k)}\, e^{-ikx_0}\;\psi_{1,k}(x), \qquad \theta(k) = -\frac{1}{F}\int_0^{k}\!\bigl[E_1(k') - \bar E_1\bigr]dk', \]

with $g$ a Gaussian envelope centred at $k_0 = -\pi/2a$ and the orbit centred on an automatically chosen Wannier–Stark rung (see below). The matched width $\sigma_k^* = \sqrt{F/(2v_\text{max})}$ ($v_\text{max}$ = peak band-1 group velocity) balances the real-space envelope $1/(2\sigma_k)$ against the orbit shear $v_\text{max}\sigma_k/F$: it is the width of the true coherent state, with minimum breathing (max/min width $\approx\sqrt2$ over one $T_B$, $\Delta x\cdot\Delta k$ within $\sqrt2$ of $1/2$). The figure ships exactly this matched width; a boosted $\sigma_k = b\,\sigma_k^*$ would squeeze the state and make it breathe by $\sim b^2$ (with $b=2$: $\Delta x$ from $1.4\,a$ to $4.8\,a$ over the period). The Bloch functions $\psi_{1,k}$ are built from the eigenvectors of the same plane-wave Jacobi solve used for the display bands, gauge-fixed to $u_{1,k}(0) > 0$, so the packet is band-1 pure. The orbit phase $\theta(k)$ places each $k$-component on the classical trajectory $x(k) = [E_1(k)-\bar E_1]/F$ and is what keeps the shear bounded — with the opposite sign or no phase the maximum $\Delta x$ grows to $8$–$12\,a$. Revival fidelity after one $T_B$: $0.999$.
Split-step Fourier propagation. Each substep applies the potential and kinetic operators in their diagonal bases,

\[ \psi \;\leftarrow\; e^{-iV\Delta t}\,\psi, \qquad \tilde\psi \;\leftarrow\; e^{-i k^2\Delta t/2}\,\tilde\psi, \]

with the FFT carrying $\psi$ between real and reciprocal space (a radix-2 Cooley–Tukey transform). Because the kinetic step is exact in $k$-space there is no finite-difference dispersion, so the wavefunction stays clean (this is what fixes the “noisy” look of a stencil propagator). A smooth imaginary absorbing layer in the outer $8\,a$ of each edge stops the periodic FFT from wrapping. Because this is the full Hamiltonian (not a single band), lowering $V_0$ makes the gap small enough for the packet to leak across it — Zener tunnelling appears naturally.
Bands and folded momentum density. The three lowest bands $E(k)$ of the field-free lattice are found at each $k$ by diagonalising the plane-wave Hamiltonian $H_{ml}(k) = \tfrac12(k+G_m)^2\delta_{ml} - \tfrac{V_0}{2}\delta_{|m-l|,1}$ ($G_m = 2\pi m/a$, $N_\text{PW} = 7$) with the same symmetric Jacobi eigensolver used by the NFE figures (Pipeline 2 — Ca / Al (15-PW Hamiltonian + marching cubes)), drawn over three Brillouin zones $k \in [-3\pi/a, 3\pi/a]$ in the Fig 3.7 conventions: all bands in the NFE band blue, over the dashed grey folded free-electron parabola $E = (k+G)^2/2$, with the energy axis in the dimensionless units $2mEa^2/\pi^2\hbar^2$; the energy window includes the bottom of band 3, whose steep free-electron arms are clipped at the top of the frame. The same band edges feed the allowed-band stripes of the real-space panel (bands 1, 2 and 3 shaded in light grey, each tilted by $-Fx$; the gaps stay unshaded), with the packet’s oscillation confined inside the tilted band-1 stripe. Each frame one clean FFT of the state gives the folded quasimomentum density $\sum_G |\psi(k+G)|^2$ on the $N_\text{cells}$ FBZ bins (every FFT bin is accumulated onto its reduced-zone residue); it is drawn in the repeated-zone scheme — solid inside the FBZ, and as fainter replicas translated by $\pm 2\pi/a$ outside, to emphasise that the quasimomentum is defined modulo a reciprocal-lattice vector.

The render loop is gated by an IntersectionObserver and the Page Visibility API, so it stops stepping when scrolled off-screen or the tab is hidden.

Parameters

Symbol	Meaning	Default	Range / value
$V_0$	cosine lattice depth	3.0	fixed (no slider)
$F = eEa$	uniform field (tilt), in $\hbar^2/ma^2$	0.25	slider 0.15 – 0.45 (re-seeds the packet)
$x_c$	orbit centre	$0$	the oscillation is centred in the frame (exactness pass below)
$\sigma_k$	seed envelope width	$\sqrt{F/2v_\text{max}}$ (matched)	clamped to $[2.5\,\delta k,\ 0.8/a]$
speed	split-step substeps per frame	10	slider 1 – 20
$x$ view	real-space window (FIXED)	$\pm 25\,a$	independent of $F$
$E$ view	energy window (FIXED vs $F$)	$\pm(V_0 + 6)$	set by $V_0$ alone
$N$	real-space grid points (FFT)	2048	fixed ($L = 90a$)
$M$	$k$-grid points (reference bands)	300	fixed
$N_\text{PW}$	plane waves for the band solve	7	fixed

The real-space axes are fixed: neither window rescales with $F$, so changing the field leaves the lattice period on screen unchanged and only steepens the visible washboard slope, while the dashed $\langle E\rangle$ line moves with Zener decay. The orbit is centred in the frame ($x_c = 0$): the packet’s swing $\pm W/(2F)$ — the band-1 allowed window along its energy line — is symmetric about the middle of the $x$ axis at every allowed $F$. Two corrections make this exact. First, the launch carries the offset $x_0 = x_c - [\tfrac12(E_\text{min}{+}E_\text{max}) - \bar E_1]/F$: the orbit-phase reference $\bar E_1$ (the $k$-average, required for the periodicity of $\theta$) is not the band midpoint, and without this offset the swing would be visibly asymmetric. Second, an exactness pass: the synthesised packet’s $\langle E\rangle$ carries small envelope-curvature corrections ($\sim\sigma_k^2$), so it is measured once and the launch rigidly shifted by $\delta = (\langle E\rangle - E_\text{mid})/F$, pinning the classical orbit centre at $x = 0$. Over the $F$ range $0.15$–$0.45$ the swing $\pm W/(2F)$ grows from $\pm 4\,a$ up to $\pm 12\,a$, always well inside the $\pm 25\,a$ window and clear of the absorbers.

References: Wannier, Phys. Rev. 117, 432 (1960) (Stark ladder); Hartmann et al., New J. Phys. 6, 2 (2004) for a modern review of Bloch oscillations; Feit, Fleck & Steiger, J. Comput. Phys. 47, 412 (1982) for the split-step Fourier integrator.

Fig. 8.3 — Zener tunnelling and the complex band structure

Source: _interactive/ch04-zener.js.

Two side-by-side panels of widths 1 : 3 and equal height give the two complementary pictures of interband tunnelling: the complex band structure in $k$-space on the left (a narrow, tall rotatable 3-D view), the field-tilted bands in real space on the right (Canvas 2D, fixed field — no controls).

Left panel — Kronig–Penney complex band structure [Three.js]

This is the algorithmically new part. The Kronig–Penney dispersion relation, in units $a = 1$ and $\tilde E = 2mEa^2/\pi^2\hbar^2$ (the Fig 3.7 convention), reads

\[ \cos(ka) = f(E) = \cos(qa) + P\,\frac{\sin(qa)}{qa}, \qquad qa = \pi\sqrt{\tilde E}, \]

with comb strength $P = 3\pi/2$. Where $|f| \le 1$ the quasimomentum is real, $k = \pm\arccos f$ (both branches drawn in blue): these are the first four allowed bands, with band edges pinned at $\tilde E = n^2$ on one side of each gap (there $qa = n\pi$ and $f = (-1)^n$ exactly). Where $|f| > 1$ the dispersion relation is still satisfied, but only by a complex quasimomentum

\[ k = k_0 \pm i\kappa, \qquad \kappa = \operatorname{arccosh}|f|, \]

with $k_0 = 0$ where $f > +1$ (zone-centre gaps and the region below band 1) and $k_0 = \pm\pi$ where $f < -1$ (zone-boundary gaps). The panel renders this as a rotatable 3-D scene (Three.js, loaded on demand from a CDN via an import-map, with OrbitControls): energy $E$ is the vertical axis — passing through the origin of the $k$-plane — and the horizontal plane carries $\mathrm{Re}\,k$ and $\mathrm{Im}\,k$. Everything is enclosed in a tall parallelepiped of square footprint $2\pi/a \times 2\pi/a$, so $\kappa$ is drawn on the same scale as $k$ and no clamping is needed (the deepest tail, $\kappa = \operatorname{arccosh}(1{+}P) \approx 2.43$, stays inside). The real bands are blue curves in the $\mathrm{Im}\,k = 0$ plane; both complex-conjugate branches $\pm i\kappa$ are always drawn (in red) lifting into the two $\pm\mathrm{Im}\,k$ directions from the pinning line $k_0$, forming the characteristic loops that bridge each gap from band top to band bottom. Six bands are computed but the energy window stops just above the top of band 4: higher bands are clipped at the box top (band 5 enters as a short stub, the gap-4 loop is complete) rather than rescaling the view. Band edges are located by bracketing and bisecting the crossings $f = \pm 1$ (50 bisection steps). All curves are evaluated in closed form and built once as fat-line geometries (Line2 + LineMaterial, which honour a pixel linewidth where the classic THREE.Line does not); the render loop only runs OrbitControls damping and is gated by an IntersectionObserver + the Page Visibility API (§D-bis) so it idles when off-screen.

/* ╔══════════════════════════════════════════════════════════╗
   ║  ALGORITHM: Kronig–Penney complex band structure          ║
   ║  f(E) = cos(qa) + P·sin(qa)/(qa), qa = π√Ẽ (Ẽ=2mEa²/π²ℏ²).║
   ║  |f|≤1 → real k = ±acos(f)/a (bands, blue). |f|>1 → gap:  ║
   ║  k = k0 ± iκ with κ = arccosh(|f|)/a, k0 = 0 for f>+1 and ║
   ║  k0 = ±π/a for f<−1 (evanescent branches, red).           ║
   ╚══════════════════════════════════════════════════════════╝ */

Right panel — bent (field-tilted) bands in real space [Canvas 2D]

The two-band spectrum is tilted by the uniform field, so the two gap edges become straight lines in real space, bending downhill with the same sign convention as the Bloch-oscillator figure:

\[ E_v^\text{top}(x) = -eEx - E_0, \qquad E_c^\text{bot}(x) = -eEx + E_0, \]

with half-gap $E_0$ and gap $E_g = 2E_0$. The allowed regions carry the same light-grey shading as the band stripes of the Bloch-oscillator figure, extending indefinitely above the conduction-band bottom and below the valence-band top (there is no outer band edge in this picture); the two gap edges are thin black lines, again matching Fig 4.2. At a fixed energy the valence-band top and conduction-band bottom are separated by the classically forbidden window of width

\[ d = \frac{E_g}{eE}, \]

marked by two ticks at the turning points $x = \pm d/2$ on the zero-energy axis. On the valence side of that axis an incoming electron is drawn as a filled arrow (yellow outline and shading, the Fig 3.3 $|\psi|$ yellow) pointing at the forbidden window; on the conduction side a yellow question mark asks whether the electron makes it across the gap — the quantitative answer being carried by the evanescent branches of the left panel and by the WKB estimate in the text. The vertical energy scale is fixed (decoupled from the field tilt) so that the gap $E_g$ occupies about half the panel height. The field is fixed at $eEa/E_g = 0.20$ (no slider); the complex band structure on the left is field-independent anyway. The right panel is a static Canvas 2D drawing, repainted on resize via the shared ResizeObserver + RAF pattern.

Parameters

Symbol	Value / default	Meaning
$P$	$3\pi/2$	Kronig–Penney comb strength
$N_\text{bands}$	6 computed	bands 1–4 in window, higher clipped
$E_\text{TOP}$	band-4 top $+\,2$ ($= 18$)	energy window (clipping plane)
$S_X / S_Y / S_Z$	0.5 / 1.6 / 0.5	3-D half-extents (Re k, E, Im k); footprint square $2\pi/a$, drawn at half the $E$ scale
line widths	2.5 px curves / 1.2 px frame	fat lines via `Line2`/`LineMaterial`
$E_0$	0.5	half-gap; $E_g = 2E_0 = 1.0$ (right panel)
$E_\text{VIS}$	$\pm 1.10$	right-panel $E$ window (gap ≈ 45 %)
widths / height	1 : 3 / half the row width	panel layout
$X$ window	$\pm 4.8$	right-panel real-space ($x/a$)
$eEa/E_g$	0.20 (fixed), slope $-eE$	field strength (right panel)

References: Zener, Proc. R. Soc. Lond. A 145, 523 (1934); Kohn, Phys. Rev. 115, 809 (1959) for the analytic continuation to complex band structure; Ashcroft & Mermin, ch. 12.

Fig. 5.5 — NFE Fermi surfaces of FCC

Source: _interactive/ch04-nfe-fcc.js.

The most algorithmically dense figure of the book. Three pseudo- potential presets are exposed via segmented selector (Free electron / Ca / Al), and the figure swaps between two parallel pipelines under the hood — analytic for the empty lattice, full diagonalisation for the others.

The 15-plane-wave model

The Hamiltonian is built in a finite plane-wave basis $|\,\mathbf{k} + \mathbf{G}_i\,\rangle$ truncated to the first three reciprocal-lattice shells: $\mathbf{G}_0 = 0$, the eight $\langle 111 \rangle \cdot 2\pi/a$ vectors, and the six $\langle 200 \rangle \cdot 2\pi/a$ vectors — 15 plane waves total. On this basis,

\[ H_{ij}(\mathbf{k}) = \frac{\hbar^2}{2m}|\mathbf{k}+\mathbf{G}_i|^2 \,\delta_{ij} + V_{|\mathbf{G}_i - \mathbf{G}_j|}\,(1-\delta_{ij}), \]

where the off-diagonal Fourier components $V_G$ are non-zero only when $\mathbf{G}_i - \mathbf{G}_j$ falls in the $\{111\}$ or $\{200\}$ shell — giving exactly two material parameters, $V_{111}$ and $V_{200}$, the Ashcroft pseudopotential coefficients. The preset table:

Preset	$V_{111}$ (eV)	$V_{200}$ (eV)	Source
Free electron	0	0	empty-lattice limit
Ca	0.34	0.68	Animalu–Heine 1965
Al	0.244	0.765	Ashcroft 1963

The lattice parameter is the FCC a = 0.405 nm of aluminium for all presets (kept fixed so the BZ geometry doesn’t change between presets and the user can compare iso-surfaces directly).

Pipeline 1 — Free electron (analytic, Harrison construction)

When $V_{111} = V_{200} = 0$, the Hamiltonian is diagonal and the band-$n$ iso-surfaces inside the FBZ are pieces of spheres centred on $-\mathbf{G}$ (the standard Harrison construction):

Band 1: a single sphere centred on $\Gamma$, radius $k_F = \sqrt{2 m E_F}/\hbar$.
Band 2: 14 caps from $\mathbf{G} \in \{\langle 111 \rangle, \langle 200 \rangle\} \cdot 2\pi/a$, each clipped to its Voronoi sub-region inside the FBZ — so caps stop exactly where two neighbouring $\mathbf{G}$’s are equidistant from $\mathbf{k}$.
Band 3: 14 × 13 ordered pairs $(\mathbf{G}_b, \mathbf{G}_a)$, one cap per pair from the sphere around $-\mathbf{G}_b$ clipped to the rank-2 wedge “$\mathbf{G}_a$ is closer than $\mathbf{G}_b$, every other $\mathbf{G}_c$ is farther than $\mathbf{G}_b$”.

Geometrically each cap starts as a THREE.SphereGeometry (UV- tessellated) and is then clipped at the mesh level by Sutherland–Hodgman 3D triangle clipping against:

the 14 FBZ planes (6 X-faces + 8 L-faces of the truncated octahedron),
the 13 Voronoi sub-region planes (band 2) or the 1 + 12 wedge planes (band 3).

/* Sutherland–Hodgman clip of a single triangle by one half-space */
for (let k = 0; k < 3; k++) {
  if (sCur >= 0) poly.push(vCur);                 /* keep visible vertex */
  if ((sCur >= 0) !== (sNxt >= 0)) {              /* edge crosses plane  */
    const t = sCur / (sCur - sNxt);
    poly.push(lerp(vCur, vNxt, t));               /* exact intersection  */
  }
}

Vertex normals are then assigned analytically as $\hat{\mathbf{n}} = (\mathbf{v} - \mathbf{c}_{\text{sphere}})/k_F$, which is the exact outward normal of a sphere at any vertex — even for the new vertices generated on the cut edges. The result is shading nearly indistinguishable from a true analytic surface even with very modest UV-tessellation density.

References: Harrison, Pseudopotentials in the Theory of Metals (Benjamin, 1966) — Harrison construction; Sutherland & Hodgman, Reentrant Polygon Clipping, Comm. ACM 17, 32 (1974).

Pipeline 2 — Ca / Al (15-PW Hamiltonian + marching cubes)

When $V_{111}$ or $V_{200}$ is non-zero the Hamiltonian is no longer diagonal and analytic Harrison spheres no longer apply. We then diagonalise the $15\times15$ Hermitian Hamiltonian at every $k$ on a regular grid using the symmetric Jacobi eigensolver (sweeps of plane rotations annihilating the largest off-diagonal element). The matrix is strongly diagonal-dominant for typical pseudopotential strengths, so 2–3 sweeps reach the tol = 1e-10 threshold:

function jacobi(H) {
  for (let sweep = 0; sweep < 50; sweep++) {
    let off = 0;
    for (p, q above diagonal) off += H[p,q]**2;
    if (off < tol) break;
    for (p, q above diagonal) rotate(p, q);   /* annihilate H[p,q] */
  }
  return diagonal(H).sortAscending();
}

The lowest 3 eigenvalues at every grid point are stored into a 3- band tensor bandGrid[b][i,j,k]. As in cubium-dos, we use the $O_h$ symmetry (factor ~48) to evaluate the diagonalisation only on the canonical wedge of the positive octant and splat onto the permutations.

Iso-surfaces at the user-selected $E_F$ are then extracted with Marching Cubes on the wedge octant and mirrored to the full BZ. A small-but-important detail: vertex normals are computed as the numerical gradient of the band on the trilinear grid $\nabla_\mathbf{k} E_n$ via central differences, not by averaging face normals. This recovers shading nearly as smooth as the analytic- sphere case, even for the deformed (non-spherical) Fermi surfaces of Ca and Al.

References: Numerical Recipes §11.1 (Jacobi eigensolver), §3.4 (trilinear interpolation), §13.1 (FFT context for plane-wave expansion); Lorensen & Cline 1987 (op. cit.); Ashcroft, Phys. Lett. 4, 202 (1963) and Animalu & Heine, Phil. Mag. 12, 1249 (1965) for the pseudopotential coefficients.

Right-panel readouts

Per band the right panel shows:

A vertical bar with the band’s energy range $[\,E_n^{\min}, E_n^{\max}\,]$, split by the iso-energy line into an occupied (orange, $E < E_F$) lower part and an empty (blue, $E > E_F$) upper part. Bands whose maximum exceeds the visible cap (drawn at $E$ corresponding to $N = 4$ e⁻/cell) are rendered with no top border and a fade to white, marking that the band continues off-chart.
A 5 eV vertical scale bar near the band-3 column.
An in-canvas read-out of the free-electron filling $N(E) \approx \int_0^E g_{\text{free}}(E')\,dE'$.
“Magnetic” snap of the iso-energy slider to the targets $N = 1, 2, 3$ e⁻/cell during drag, so the textbook fillings are easy to land on.

Fig. 5.6 — Fermi surfaces of real metals

Source: _interactive/ch05-metals-fs.js (viewer) and _interactive/_test/gen-metals-fs.mjs (offline generator; the JSON it produces is injected into the viewer as DATA, ~1.3 MB with nine materials). A companion validator, _test/check-metals-fs.mjs, measures the physical targets quoted below.

Physical models

Li, Na, K (BCC), Ca and Al (FCC) — local empirical pseudopotential on plane waves: $H_{\mathbf g\mathbf g'} = \frac{\hbar^2}{2m}|\mathbf k+\mathbf g|^2\delta_{\mathbf g\mathbf g'} + V(|\mathbf g-\mathbf g'|)$, with two symmetric form factors per element. Al uses the canonical dHvA-fitted values $V_{111}=0.0179$ Ry, $V_{200}=0.0562$ Ry [Ashcroft, Phil. Mag. 8, 2055 (1963)]; Ca the Animalu–Heine model-potential values $V_{111}=0.0250$ Ry, $V_{200}=0.0500$ Ry (0.34/0.68 eV — the same pair used by the “Ca” preset of Fig. 5.5, Pipeline 2 — Ca / Al (15-PW Hamiltonian + marching cubes)); Na uses the dHvA-fitted $|V_{110}| = 0.0165$ Ry [Lee, Proc. Roy. Soc. A 295, 440 (1966)]; Li the positron-annihilation value $V_{110} \simeq 0.10$ Ry [Donaghy & Stewart, Phys. Rev. 164, 391 (1967)] — Li is the “least NFE” alkali; K a representative small value (its measured surface is spherical to $\sim 0.1\%$). Second form factors are minor.
Cu, Ag, Au, Pd (FCC) — minimal $s{+}d$ tight binding, $6\times6$ per $\mathbf k$: one $s$ orbital ($E_s = 0$ reference) plus the five $d$ orbitals, first-neighbour hoppings via the full Slater–Koster $dd(\sigma,\pi,\delta)$ table and $sd\sigma$ hybridisation. Hoppings start from the “modified Harrison” universal values of Shi & Papaconstantopoulos [PRB 70, 205101 (2004)] evaluated at the nearest-neighbour distance; $E_d$ (and small hopping adjustments) are then tuned — in the true semi-empirical spirit — to the experimental anchors of the Cu band structure collected by Courths & Hüfner [Phys. Rep. 112, 53 (1984)]: $d$ complex within $[-5.2,-2.0]$ eV of $E_F$, width $3.2$ eV, $L_{2'}$ at $-0.85$ eV, $\Gamma_1$ at $-8.6$ eV. Values in eV:

Table 1: $s{+}d$ parameters (eV); $E_s = 0$ is the reference.

	$E_d$	$V_{ss\sigma}$	$V_{dd\sigma}$	$V_{dd\pi}$	$V_{dd\delta}$	$V_{sd\sigma}$
Cu	−2.10	−0.68	−0.414	0.246	−0.045	−0.47
Ag	−4.15	−0.58	−0.422	0.251	−0.046	−0.42
Au	−2.55	−0.73	−0.645	0.383	−0.070	−0.52
Pd	−1.60	−0.74	−0.664	0.394	−0.072	−0.55

Algorithm

Eigensolvers. EPM: real-symmetric Householder + implicit-shift QL (tred1/tqli), ~55–65 plane waves on the path, reduced cutoff on the FS grid. TB: the complex-Hermitian $6\times6$ is embedded into the real-symmetric $12\times12$ $[\mathrm{Re},-\mathrm{Im};\mathrm{Im},\mathrm{Re}]$ form; eigenvalues come out in degenerate pairs.
Fermi energy by bisection on the grid occupancy: $E_F$ such that $2\sum_{n\mathbf k}[E_n(\mathbf k)<E_F]/N_{\mathbf k} = n_e$, with the $\mathbf k$ sum restricted to the FBZ interior.
Fermi surfaces by shared-edge marching tetrahedra: each band is sampled on the positive octant only of a $49^3$ box grid ($25^3$ points — under $O_h$, $E(\mathbf k)$ is even in every component, the same trick as the Fig. 5.5 pipeline); every cube is split into 6 tetrahedra and the iso-surface $E_n(\mathbf k)=E_F$ is triangulated by linear interpolation on the edges (16 trivial cases — no marching-cubes tables). Each iso-vertex is generated once per global grid edge (canonical node-pair key), so the mesh is watertight by construction; the diamond-shaped tetrahedral ridges are then relaxed by 10 alternated Taubin smoothing passes ($\lambda=0.5$, $\mu=-0.53$ — shrink-free), with vertices born on the mirror planes or on the box faces pinned there so the mirrored copies and the GPU clip stay seamless. Triangles lying entirely beyond the FBZ plus a one-cell margin are dropped — the exact FBZ cut is left to the viewer’s GPU clipping planes. Vertices are quantised to integers (1/250 of $2\pi/a$).
FBZ wireframe generated from the half-space set itself: vertices = triples of planes meeting inside the zone, edges = vertex pairs sharing two planes.
Viewer. Three.js: BZ wireframe and opaque MeshPhysicalMaterial sheets (metalness 0.75, roughness 0.22, clearcoat, pastel palette and ambient + hemisphere + directional lighting — exactly the Fig. 5.5 recipe, no environment map) clipped to the FBZ by per-material GPU clippingPlanes built from the shipped half-spaces (normals must be unit vectors — Plane.applyMatrix4 assumes it). The octant mesh is expanded $\times 8$ by the $O_h$ mirrors (triangle winding reversed on odd-parity octants), coincident vertices are welded, and computeVertexNormals then yields seam-free shading; the $k$-path as a tube with node spheres and crossing markers where $E_n(\mathbf k)=E_F$ (interpolated along the path), and a purple marker synchronised with the draggable cursor of the canvas band panel. Rendering is on-demand (§D-bis gating).

Verification

Al: $E_F-\Gamma_1 = 11.6$ eV (free-electron 11.7), X gap $= 1.49$ eV $\approx 2V_{200} = 1.53$ eV; two sheets with the classic second-band-hole / third-band-monster topology.
Alkali: single sphere; Li asphericity 5.0% on the mesh, matching the measured maximum anisotropy $4.8\pm0.3\%$ [dHvA/QMC, arXiv:1912.12295]; Na and K within the grid resolution of spherical. Gaps at N: 1.71/0.40/0.15 eV for Li/Na/K.
Cu: $d$ complex $[-5.0, -1.8]$ eV vs $E_F$, width 3.20 (experiment $[-5.2, -2.0]$, width 3.17); $\Gamma_1 - E_F = -9.2$ (exp $-8.6$); $sp$ level at $L$ at $-0.49$ eV (exp $-0.85$) → open necks through the hexagonal faces (mesh verified to touch them). Ag/Au analogous with their own $d$ depths (Au widest, no spin–orbit in the model).
Pd: $E_F$ falls just below the $d$-band top ($+0.23$ eV above $E_F$): four sheets — $\Gamma$-centred electron sheet, the open hole “jungle gym” reaching the zone faces, and small hole pockets — the known sheet structure of Pd [dHvA, Vuillemin & Priestley 1965].

Fig. 5.1 — Zone folding of the 2D free-electron parabola

Source: _interactive/ch05-folding2d.js (self-contained — no offline generator; everything is analytic).

Physical model

A free electron on the 2D square lattice (lattice constant $a$). The “band structure” is the single parabola replicated on the reciprocal lattice: one copy per $\mathbf G = (2\pi/a)(m,n)$,

\[E_\mathbf{G}(\mathbf k) = \frac{\hbar^2|\mathbf k - \mathbf G|^2}{2m}, \qquad\text{i.e.}\qquad \frac{E}{E_0} = (\tilde k_x - 2m)^2 + (\tilde k_y - 2n)^2,\]

with $\tilde{\mathbf k} = \mathbf k/(\pi/a)$ and $E_0 = \pi^2\hbar^2/2ma^2$ — the same units as Figure 5.2. Only the nine copies $m,n \in \{-1,0,1\}$ are kept, grouped by shell (centre / four edge / four corner copies: blue / red / green — symmetry-related copies give degenerate folded branches). The energy cap is $E_{\max} = 9\,E_0$: this is exactly the completeness threshold of the 3×3 set, since the first missing branch (from $\mathbf G = (2,0)\cdot 2\pi/a$) touches $E = 9\,E_0$ at $X$ — below the cap the folded spectrum is complete.

Scene and interaction (Three.js + Canvas 2D)

3D scene: one shared paraboloid mesh (radial × angular grid, 36×72, closed by a lid disc at the cap) instanced at the nine $\mathbf G$ centres with the house metallic-translucent material (PMREM RoomEnvironment); the $k$-plane is tiled by 3×3 Brillouin zones with the $\Gamma$–$M$–$X$–$\Gamma$ path marked in the central zone. Path and section trajectories are drawn as TubeGeometry over a custom polyline Curve with sphere joints at the ends — 1-px GL lines would be nearly invisible against the surfaces. The trajectory of each copy, $E_\mathbf{G}(\mathbf k)$ lifted onto its paraboloid above the path, is split per leg and clipped at the cap.
Copy state machine (the 3×3 button grid, arranged like the $\mathbf G$-lattice): each copy is off, included (shaded button — section curve and 3D trajectory kept, surface hidden) or active (the ONE paraboloid displayed). Activating a copy demotes the previous active one to included; clicking a shaded/active button switches it off.
Section cuts: clicking one of the three legs in the section panel highlights it and (a) raises a translucent vertical wall over that leg only, up to the cap, with vertical border tubes at its ends; (b) slices every bowl with a GPU clipping plane (renderer.localClippingEnabled): $\Gamma$–$M$ keeps $k_x<k_y$, $M$–$X$ keeps $k_x>\pi/a$, $X$–$\Gamma$ keeps $k_y<0$. The trajectory tubes are not clipped — they lie on the cutting plane and mark the cut edge, i.e. the folded band; only the selected leg’s tubes stay visible.
Right panel: the cross-section along $\Gamma$–$M$–$X$–$\Gamma$ in the same style as the band panel of Figure 5.2 (energy axis on the right), 70 samples per leg; the active copy’s branch is drawn thicker.
The OrbitControls render loop is gated on viewport and tab visibility (the §D-bis house pattern).

Verification

The folded energies at the path nodes reproduce the analytic $V_0=0$ spectrum of the cellular solver: $\Gamma$: $0, 4, 8$; $M$: $2$ (four-fold: central + two edge + one corner copy); $X$: $1, 1, 5$ (all in units of $E_0$).

Fig. 5.2 — The 2D cellular band solver

Source: _interactive/ch05-cellular2d.js (rendering + live solver) and _interactive/_test/gen-cellular2d.mjs (offline strength sweep, injected into the figure file).

Physical model

A single electron in one square cell $[0,a]^2$ with an inverted-Gaussian confining well,

\[V(\mathbf r) = V_0 E_0\left[1 - e^{-\frac12\left(x_r^2/\sigma_1^2 + y_r^2/\sigma_2^2\right)}\right],\]

where $(x_r,y_r)$ are cell coordinates rotated by the ellipse angle $\theta$ (the two draggable handles set $\sigma_1,\sigma_2,\theta$, with snaps to a circle and to the crystal axes). Bloch’s theorem enters as Floquet boundary conditions on the cell edges, $\psi(a^-,y) = e^{ik_x a}\psi(0^+,y)$ (and likewise in $y$): wrapping the finite-difference stencil with these $\mathbf k$-dependent phases turns the cell into a complex-Hermitian eigenproblem $H(\mathbf k)$ whose eigenvalues along $\Gamma$–$M$–$X$–$\Gamma$ are the bands. Natural units $\hbar = m = a = 1$; energies in $E_0 = \pi^2\hbar^2/2ma^2$.

Algorithm

Discretisation: 5-point stencil on a $16\times16$ grid ($N = 256$); the Bloch phases sit on the wrap-around links. $H\psi$ is applied matrix-free in $O(5N)$.
Eigensolver: Lanczos with full reorthogonalisation on the complex field (real/imaginary array pair), $m$ up to 140 vectors. The seed is real: at the points where $H(\mathbf k)$ is real ($\Gamma$, $X$, $M$) the states then come out exactly real, as they physically must. The Lanczos tridiagonal is solved by QL with implicit shifts — values-only (tqli) for band curves, with eigenvectors (tql2) for the two displayed wavefunctions (measured ≈ 3–4× faster than the cyclic Jacobi it replaced, which is kept as a convergence fallback). Residuals $\|H\psi - E\psi\| \sim 10^{-12}E_0$ at $m=140$, including the degenerate bands at $M$.
Free-electron limit: at $V_0 = 0$ Lanczos cannot resolve the massive degeneracies (a single-vector iteration converges one state per degenerate subspace), so the bands are taken from the exact analytic finite-difference dispersion — the grid plane waves give $E = N_{pt}^2\,(2 - \cos\beta_x - \cos\beta_y)$ with $\beta_i = (k_i + 2\pi m_i)/N_{pt}$; enumerating all $N_{pt}^2$ modes captures every degeneracy and the folded parabola stays continuous.
Precompute + live solve: the default circular well is swept offline (17 strengths × 133 $k$-points × 7 bands, ≈ 97 kB injected), so the strength slider responds instantly; symmetry-broken wells are re-solved live in the browser, with the $k$-path chunked over animation frames (≈ 12 ms per frame) so the UI never blocks.
Wavefunction continuity (fixed band index): while $k$ or the well is dragged the two displayed states keep their band index; each new state is the projection of the previously displayed one onto the numerically degenerate cluster around $E_{band}$ (single-member cluster → pure phase alignment; projection too small → clean restart). A pair guard re-orthogonalises $\psi_B$ against $\psi_A$ (Gram–Schmidt inside the cluster, still an exact eigenstate) whenever the two distinct states meet in the same degenerate cluster — without it the pair could collapse onto one state at $M$.
Display: $\mathrm{Re}\,\psi$ / $\mathrm{Im}\,\psi$ with a diverging colormap (white = 0, red > 0, blue < 0) or the full complex field by domain colouring (hue = phase, pastel saturation = $|\psi|$), tiled over $1\times1$/$3\times3$/$5\times5$ cells with the Bloch phases $e^{i\mathbf k\cdot\mathbf R}$ applied per cell.

Verification

Generator: at $V_0=0$ the lowest band spans $1.99\,E_0$ (free-electron folding); at $V_0=9\,E_0$ it narrows to $0.35\,E_0$ with a $2.82\,E_0$ gap to the next band (tight-binding limit).
$V_0=0$ analytic bands: max adjacent-$k$ jump $0.18\,E_0$ (continuous).
Degeneracy physics: at $M$ the circular well gives bands 2–3 exactly degenerate ($5.970, 5.970\,E_0$); stretching the ellipse to $\sigma_1/\sigma_2 = 1.35^2$ splits them by $0.87\,E_0$.
Displayed pair: $|\langle\psi_A|\psi_B\rangle| \sim 10^{-16}$, and it stays $< 10^{-3}$ when the pair is dragged through the $M$ degeneracy (the collapse this guards against reached overlap 1.000 before the fix).

Table 2: Cellular-solver parameters (constants in the figure source).

Parameter	Value
grid / bands kept	$16\times16$ / 7
$k$-path	$\Gamma$–$M$–$X$–$\Gamma$, 133 points (44 per leg)
strength sweep	$V_0 = 0\ldots9\,E_0$, 17 levels
default well	$\sigma_1=\sigma_2=0.14\,a$, $\sigma\in[0.09, 0.34]$
snaps	$7^\circ$ to the axes, $0.025$ to a circle
Lanczos $m$ (ψ settled / drag)	140 / 120
Lanczos $m$ (live bands, aniso / near-circular)	110 / 140
continuity	cluster width $0.006\,E_0$, restart below $\\|P\,\psi_{prev}\\|=0.35$

Figs. 5.4 / 10.1 — EPM band structures (Si, Ge, GaAs, InAs)

Source: _interactive/ch05-epm-bands.js (the band-review explorer, Fig. 5.3), _interactive/ch06-gaps.js (the two-panel Si|GaAs classic of the semiconductor chapter, Fig. 10.1) and _interactive/_test/gen-epm-classic.mjs (offline band generator, shared by both figures; four materials).

Fig. 5.3 adds an interactive front-end to the same precomputed data: a segmented selector (free / Si / Ge / GaAs / InAs), the Cohen–Bergstresser form-factor table of the selected material, and a 3Dmol view of its conventional cubic cell (diamond/zincblende, built with the same FCC-plus-tetrahedral-basis construction and ball-and-stick style as the chapter-1 structures figure; an empty cell frame for the free-electron preset), next to the band panel drawn in exactly the classic style described below. The transition between presets is a physical morph: the form factors and the lattice constant are linearly interpolated and the same 65-PW Hamiltonian is re-diagonalised in the browser on a coarse k-mesh (Householder + QL eigenvalues on the $2N$ real embedding, with the U|K and diamond-$X$ blends applied per frame), so every intermediate frame is the true spectrum of a physical Hamiltonian and anticrossings open and close correctly; the final frame snaps back to the precomputed fine mesh. For the free-electron preset the parabola minimum is pinned at $-12$ eV (commensurate with the VB bottoms of the real materials) and no occupation overlay is drawn.

Physical model

Empirical pseudopotential method (EPM) in the Cohen–Bergstresser parametrisation: the whole crystal potential of the diamond/zincblende structure is reduced to three symmetric (plus, for zincblende, three antisymmetric) Fourier coefficients on the shells $|G|^2 = 3, 8, 11$ in units of $(2\pi/a)^2$ ($|G|^2 = 4$ enters the antisymmetric part only). The Hamiltonian over plane waves $|\mathbf k+\mathbf G\rangle$ is

\[H_{GG'} = \frac{\hbar^2|\mathbf k+\mathbf G|^2}{2m}\,\delta_{GG'} + V_S(|\Delta G|^2)\cos(\Delta\mathbf G\cdot\boldsymbol\tau) + i\,V_A(|\Delta G|^2)\sin(\Delta\mathbf G\cdot\boldsymbol\tau), \qquad \boldsymbol\tau = \tfrac{a}{8}(1,1,1).\]

No spin–orbit coupling is included — this is why the figure carries no double-group labels ($\Gamma_6/\Gamma_7/\Gamma_8$) and no split-off band.

Algorithm and precompute pipeline

Basis: all $|G|^2 \le 16$, i.e. 65 plane waves. The complex-Hermitian $65\times65$ matrix is embedded as a $130\times130$ real symmetric one and diagonalised by eigenvalue-only cyclic Jacobi (every eigenvalue appears twice; even-index picks). Eight bands are kept.
Path $L\to\Gamma\to X\to U\,|\,K\to\Gamma$ (~140 points, abscissa = cumulative $|\Delta k|$). $U$ and $K$ are equivalent points ($U = K + \mathbf G$ up to a point-group operation), so the spectrum is continuous across the junction in exact arithmetic; the finite basis breaks the identity by $$0.1–0.3 eV, and the residual per-band mismatch is reabsorbed linearly along $K\to\Gamma$ (a “gauge blend”), keeping the junction exactly continuous and the $\Gamma$ endpoint untouched. The same treatment restores the diamond $X$-point band sticking (every level at $X$ is exactly two-fold degenerate in the non-symmorphic diamond structure, a pairing the fixed $|G|^2$ cutoff breaks by $\approx 0.16$ eV on the lowest Si pair): for the diamond-structure materials the bands are pair-averaged at $X$ and the correction is blended linearly to zero towards $\Gamma$ and $U$. The real zincblende $X_1/X_3$ splitting of GaAs and InAs is untouched.
Reference: all energies are shifted so that the VBM (top of the fourth band, at $\Gamma$) sits at zero.
The whole calculation runs offline in the tracked generator; the quantised band data (~20 kB) is injected into the figure file, which therefore plots instantly (a few ms) with no in-browser diagonalisation.

Verification (printed by the generator)

Si: indirect gap $0.76$ eV with the CBM on $\Delta$ at $0.86\,\Gamma X$; VB bottom $-12.7$ eV. With the smaller 51-PW basis the CBM lands at $X$ instead (0.93 eV): the $\Delta$-minimum position is a basis-convergence detail, while the gap value drifts down with basis size — the known behaviour of the local (Löwdin-uncorrected) EPM.
Ge: indirect gap $0.71$ eV with the CBM at $L$; VB bottom $-12.1$ eV.
GaAs: direct gap $1.41$ eV at $\Gamma$; VB bottom $-12.3$ eV.
InAs: direct gap $0.44$ eV at $\Gamma$; VB bottom $-11.3$ eV.

Table 3: Cohen–Bergstresser parameters used by the generator.

Parameter	Si	GaAs
$a$ (Å)	5.43	5.64
$V_S(3), V_S(8), V_S(11)$ (Ry)	$-0.21,\ 0.04,\ 0.08$	$-0.23,\ 0.01,\ 0.06$
$V_A(3), V_A(4), V_A(11)$ (Ry)	—	$0.07,\ 0.05,\ 0.01$

Reference: M. L. Cohen and T. K. Bergstresser, Phys. Rev. 141, 789 (1966).

Fig. 4.2 — Slater–Koster two-centre integrals in 2D

Interactive illustration of the angular factors of the Slater–Koster table (Equation 4.5), restricted to the $xy$ plane.

Physical model

Two sites at distance $d$ along $\hat{\mathbf n} = (\cos\theta, \sin\theta)$ carry one orbital each, taken from the 2D set

\[\psi_s = e^{-r}, \qquad \psi_{p_x} = x\,e^{-r}, \qquad \psi_{p_y} = y\,e^{-r}\]

(lengths in units of the decay length). The left panel renders the signed product $\psi_1(\mathbf r)\,\psi_2(\mathbf r - d\hat{\mathbf n})$; the right panel its integral $S(\theta) = \int \psi_1 \psi_2\, d^2r$. Rotating the $p$ orbitals into bond-parallel ($\sigma$) and bond-perpendicular ($\pi$) components — e.g. $p_x = n_x p_\sigma - n_y p_\pi$ — and using the reflection symmetry about the bond axis (cross terms vanish) reproduces the SK angular structure exactly: $S_{s,p_x} = n_x (sp\sigma)$, $S_{p_x,p_x} = n_x^2 (pp\sigma) + n_y^2 (pp\pi)$, $S_{p_x,p_y} = n_x n_y \left[(pp\sigma) - (pp\pi)\right]$. The hopping integrals $\langle \phi_1 | V_a | \phi_2 \rangle$ of the table obey the same angular law; only the radial prefactors (and their signs) differ.

Algorithm

Colorplot: the product is sampled on a $160^2$ grid over the $(2\times 5.2)^2$ window and mapped to a diverging red/blue scale with a saturating exponent $|v|^{0.55}$ (to keep the faint overlap region visible); it is redrawn live while dragging (a few ms per frame).
Curve: $S(\theta)$ is integrated numerically on a $120^2$ grid of half-width $6.5$ at $121$ angles, then cached per orbital pair, so switching pairs back and forth is instantaneous. The curve is normalised to its largest absolute value; dashed guides mark min and max.
Sync: the bond angle $\theta$ is the single state variable — dragging either the atom (left, along the circle of radius $d$) or the marker (right, along $\theta$) updates both panels.

The separation $d = 2.6$ is chosen so that the two $p$-$p$ radial integrals have opposite signs ($(pp\sigma) = -0.12$, $(pp\pi) = +0.59$): the $p_x$-$p_x$ curve then crosses zero as $\theta$ sweeps from the $\sigma$- to the $\pi$-dominated orientation.

Verification

Automated checks on the same integrator: $S_{s,s}$ isotropic to $10^{-3}$; $S_{s,p_x}(\theta)/S_{s,p_x}(0) = \cos\theta$ to better than $1\%$; $S_{p_x,p_x}(30°)$ matches $n_x^2(pp\sigma) + n_y^2(pp\pi)$ to $1\%$; $S_{p_x,p_y}$ vanishes at $0°$ and $90°$ and equals $\left[(pp\sigma)-(pp\pi)\right]/2$ at $45°$.

Table 4: Parameters of the Slater–Koster figure.

Parameter	Value
site separation $d$	2.6
display window / colorplot grid	$\pm 5.2$ / $160^2$
integration window / grid	$\pm 6.5$ / $120^2$
$\theta$ samples	121
radial integrals at $d$	$(ss\sigma)=0.56$, $(sp\sigma)=-0.73$, $(pp\sigma)=-0.12$, $(pp\pi)=+0.59$

Symbol	Value / default	Meaning
\(P\)	\(3\pi/2\)	Kronig–Penney comb strength
\(N_\text{bands}\)	6 computed	bands 1–4 in window, higher clipped
\(E_\text{TOP}\)	band-4 top \(+\,2\) (\(= 18\))	energy window (clipping plane)
\(S_X / S_Y / S_Z\)	0.5 / 1.6 / 0.5	3-D half-extents (Re k, E, Im k); footprint square \(2\pi/a\), drawn at half the \(E\) scale
line widths	2.5 px curves / 1.2 px frame	fat lines via `Line2`/`LineMaterial`
\(E_0\)	0.5	half-gap; \(E_g = 2E_0 = 1.0\) (right panel)
\(E_\text{VIS}\)	\(\pm 1.10\)	right-panel \(E\) window (gap ≈ 45 %)
widths / height	1 : 3 / half the row width	panel layout
\(X\) window	\(\pm 4.8\)	right-panel real-space (\(x/a\))
\(eEa/E_g\)	0.20 (fixed), slope \(-eE\)	field strength (right panel)