Quantum Matter

In June 1995, in a basement in Boulder, Colorado, Eric Cornell and Carl Wieman's group cooled about two thousand rubidium atoms to 170 billionths of a degree above absolute zero — then let the cloud fly apart and photographed where the atoms went, which is a way of photographing their velocities. Above a critical temperature: a smooth thermal blob, exactly what Maxwell would have drawn. Below it: a spike erupts from the blob's center — thousands of atoms with the same near-zero velocity. Not similar velocities. The same single quantum state, occupied by a macroscopic crowd, big enough to see with a camera. Bose and Einstein had predicted that spike in 1925; it took seventy years of cryogenic ingenuity to develop the film.

That photograph is one panel of a triptych this chapter assembles. Panel two hangs in every building: ordinary solid matter, whose rigidity — we will now be able to say precisely — is Pauli pressure, the exclusion principle pushing back when you squeeze 10²³ electrons. Panel three is the strangest: mercury below 4.2 kelvin, where Kamerlingh Onnes found in 1911 that electrical resistance does not get small but becomes exactly zero — currents that circulate for years without a battery. Three materials, three behaviors, one explanation: the two signs of Chapter 15, compounded by sheer number. Fermions build seas. Bosons build spikes. And — the twist the century took 46 years to find — fermions can pair up and defect.

Chapter 15 said fermions stack; now we count what the stacking costs. Put N electrons (spin gives two per state) in a box of volume V and fill plane-wave states from the bottom. Momentum states fill a sphere of radius k_F, and counting states inside it — two per (2π)³/V of k-space — fixes the radius in terms of nothing but the density n = N/V:

Averaging ε = ħ²k²/2m over the filled sphere weights the large-k shells (there are more of them), and the average lands at exactly 60% of the top:

Put numbers in for copper: n ≈ 8.5×10²⁸ m⁻³ gives E_F ≈ 7 eV — a temperature equivalent of eighty thousand kelvin. Room temperature is 300 K. So the electron gas in every wire on Earth is, for all practical purposes, at absolute zero: frozen into its ground state, each electron moving at ~10⁶ m/s because Pauli leaves it nowhere slower to sit, with only the thin skin of states within kT of the surface free to do anything (which is why electrons add almost nothing to a metal's heat capacity — a classical scandal this one formula retires). And the pressure in (16.2) is mechanical and enormous — ~10¹⁰ pascals in copper. Matter feels solid because compressing it would push the sea level up.

Push that pressure to its limit. When a star's fusion dies, gravity compresses the core until electron degeneracy pressure — $P\propto n^{5/3}$, no temperature required — carries the entire load: a white dwarf, a star-mass of matter held up by the exclusion principle alone. Squeeze harder and the electrons go relativistic, the exponent softens to $n^{4/3}$, and above ~1.4 solar masses Pauli loses to gravity (exercise 3 walks Chandrasekhar's argument). Every white dwarf in the sky is a monument to (c†)² = 0.

Free fermions fill a sea; real electrons also repel each other, and the honest Hamiltonian with e²/r₁₂ terms has no exact solution for N > 2. But we own a tool for exactly this situation. Chapter 12's variational principle says: pick a family of trial states, minimize ⟨H⟩, and you get a rigorous upper bound plus the best caricature the family can afford. Chapter 15 tells us the natural family for fermions: Slater determinants. Minimizing over every determinant of N orbitals yields the Hartree–Fock equations: each orbital obeys a Schrödinger equation in the electrostatic field of all the others' charge clouds —

— plus a second, stranger term. The exchange term has no classical picture: it is Chapter 15's exchange hole doing energetics — parallel-spin electrons are automatically kept apart by antisymmetry, so they feel less repulsion than a classical charge smear would predict, and the determinant collects that discount. (It is the same 2J that split helium in 15.3, now working wholesale.) This one approximation built twentieth-century chemistry and atomic physics; its self-consistent loop — guess orbitals, build the field, re-solve, repeat — is why the periodic table can be computed and not just observed.

Be honest about what it misses. A single determinant lets opposite-spin electrons pass through each other uncorrelated; the energy it fails to find is called, in professional confession, the correlation energy. Usually it is a small tax — a percent or so. The entire drama of Section 16.5 is a case where that “small” leftover contains the whole phenomenon.

Now the bosons. In thermal equilibrium at temperature T and chemical potential μ — the energy cost to add one more particle to the system, the dial nature turns to hold the total particle number fixed — the mean occupation of a state at energy ε is

— one bracket-sign apart, exactly like the algebra that produced them (exercise 5 derives both from Chapter 15's n+1 and 1−n matrix elements; Feynman-style, the distributions are the operator algebra thermalized). For bosons the −1 puts a wall at μ = ε₀: push μ up to the lowest level and its occupation diverges. Here is the trap that becomes a phase transition. In a 3D harmonic trap, add up what all the excited states can hold at μ → ε₀, and the sum converges:

The number $\zeta (3)\approx 1.202$ is nothing to fear: it is the Riemann zeta function at 3, just the definite sum ${\sum }_{n=1}^{\infty }1/n^3=1+\frac{1}{8}+\frac{1}{27}+\cdots$ that falls out of turning the state-counting sum into an integral — a fixed constant like π, no more. What it multiplies is the point: the excited states have a finite carrying capacity. Cool the gas until N exceeds it and the surplus has exactly one place to go: the single lowest state. Setting N = N_exc^max defines the critical temperature, and below it the overflow is macroscopic:

No interaction, no attraction, no force pulls the atoms together — condensation is a bookkeeping catastrophe, statistics and nothing else.

That is the theory. Turning it into the photograph of Section 16.1 took the apparatus in Figure 16.1. A dilute gas of rubidium atoms is held in a magnetic (or optical) trap — a bowl-shaped potential that keeps the cloud off the walls — and chilled in two stages: crossed laser beams slow the atoms by photon recoil (laser cooling), then the trap edge is lowered so the hottest atoms boil away and the rest rethermalize colder ( evaporative cooling, the same physics that cools a coffee cup). To read the cloud's velocities, the trap is switched off and the atoms fly freely for a fixed time — time-of-flight — so a fast atom lands far from center and position becomes a map of velocity; a probe laser then casts a shadow (an absorption image) recorded by a camera. Above $T_c$ the image is a broad thermal blob; cool through $T_c$ and a sharp spike of near-zero-velocity atoms erupts from its center — the condensate, made visible.

Watch both species handle the same cold:

The contrast on screen is the whole argument. The fermion curve fills smoothly, one state at a time, no state ever holding more than its allotment; the boson ground-state occupation hugs zero and then, as T drops through $T_c$, turns sharply upward — a genuine phase transition with no force behind it, driven purely by the excited states running out of room. That upward kink is what Cornell and Wieman photographed, and its sharpness is also the frontier: a finite N rounds the corner, so how cleanly the kink appears is a direct readout of how close a real, countable cloud is to the thermodynamic ideal.

Real condensates interact, and there mean-field thinking returns in bosonic dress. If all N atoms share one wavefunction ψ, each feels the contact repulsion of the local crowd, g|ψ|² — giving a Schrödinger equation with the density fed back into the potential, the Gross–Pitaevskii equation:

The eigenvalue $\mu$ is the same chemical potential as before — the energy to add one atom to the condensate. When the repulsion is strong, the kinetic term $-\frac{{\hslash }^2}{2m}{\nabla }^2$ is negligible beside it and the equation collapses to $\mu =V_{\text{trap}}+g|\psi {|}^2$: the density simply mirrors the trap, $|\psi {|}^2=(\mu -V_{\text{trap}})/g$ wherever that is positive. This is the Thomas–Fermi limit, and the constant it fixes — ${\mu }_{\text{TF}}$, set by interaction and trap alone, with kinetic energy written off — is the ${\mu }_{\text{TF}}$ that exercise 6 pins down.

It is nonlinear, so Chapter 3's superposition comforts are gone — but Chapter 6's split-operator solver never needed linearity. Run it in imaginary time (t → −iτ turns oscillation into decay; every state dies at its energy's rate, so the ground state is the last one standing) and the condensate finds itself. Two exact facts referee the answer: at g = 0 it must be Chapter 7's Gaussian with E = ½ħω, and at every g the scaling identity

must hold at the true minimum (stretch ψ by λ, demand dE/dλ = 0 at λ = 1 — exercise 6). Here is that solver, running in your browser:

What the solver settles onto is not a plot but an object: the actual density profile of a trapped condensate, the same shape the absorption images of Section 16.1 trace on a camera. Turn the interaction up and watch the Gaussian swell into the flat-topped Thomas–Fermi cloud, pressure from the atoms' own repulsion flattening the peak. And the nonlinearity you had to give up Chapter 3's superposition for is exactly the ingredient that makes a condensate behave like a fluid — quantized vortices, solitons, superfluid flow all live in that $g|\psi {|}^2$ term. Exercise 9 turns the trap off entirely and sends the same solver after a bright soliton that holds itself together, the attraction standing in for the walls that are no longer there.

Back to Onnes's mercury. Zero resistance is not “very good conduction” — persistent currents in superconducting rings have been watched for years without measurable decay, and the flux they trap is quantized in units of h/2e (Chapter 13's exercise; hold that 2 — it is a fingerprint). For 46 years the effect defied every theory built by improving the electron gas perturbatively. The reason it had to fail is now in your hands: the answer is proportional to $e^{-1/\lambda }$, and every Taylor coefficient of that function at λ = 0 is zero. Perturbation theory wasn't approximating the answer; it was computing, correctly, order by order, the Taylor series of zero.

Two ingredients crack it. First, an attraction — between electrons? Yes: a passing electron pulls the positive lattice ions toward its track, and the ions, being heavy, are still leaning inward after the electron is long gone; a second electron feels that leftover ridge of positive charge. A retarded, phonon-carried handshake — weak, and buried under Coulomb repulsion at short range, but attractive in a thin shell of energies ~ħω_D around the Fermi surface.

Second — Cooper's 1956 calculation, the hinge of the whole story. Two electrons interacting above a frozen Fermi sea (Pauli blocks every state below E_F). In empty 3D space, a weak attraction binds nothing — there is a threshold. But above the sea the available states pile up at the Fermi surface, and the pairing equation binds at any coupling strength:

Here $V$ is the strength of the phonon-borne attraction and $N(0)$ is the density of states at the Fermi surface — the number of electron energy levels available per unit energy right at $E_F$ — so the dimensionless coupling $\lambda =N(0)V$ measures how many pairable states the attraction can actually reach. The Fermi sea makes its own surface unstable. Let all the electrons near the surface pair up — each pair a (k↑, −k↓) partnership, zero net momentum, roughly boson-like — and solve for the new ground state self-consistently (the BCS variational state of Bardeen, Cooper, and Schrieffer, 1957; exercise 8 sketches it) and the pairing opens an energy gap around the Fermi surface:

The gap is the whole ballgame. A current is a coherent motion of the paired sea; to degrade it, a scatterer must break a pair, and that costs at least 2Δ. Below T_c thermal kicks can't pay the toll, so the current has no channel for decay — resistance is not small but absent, the same way an atom in its ground state does not slowly radiate. The impersonation is complete: fermions, forbidden from condensing alone, condense as pairs — charge 2e, hence the flux quantum h/2e, the smoking gun Chapter 13 planted. One loose end stays honest: BCS explains mercury and aluminum; the copper-oxide superconductors found in 1986 pair by a mechanism still argued about — this chapter's story is not finished anywhere on Earth.

The lab (Rust-QP/ch16-matter) runs all four stories as numerics with exact referees: shell-filling (the closure numbers 2, 8, 20, 40, 70, 112 — exact combinatorics for the 3-D harmonic oscillator, and worth naming carefully: only the first three coincide with nuclear physics' magic numbers 2, 8, 20, 28, 50, 82, 126, because a real nucleus has a strong spin–orbit coupling this trap does not), box-filling (the 3/5 law to five digits at N = 10⁵), condensate curves at N = 10², 10³, 10⁵ (finite-size honesty: the deviation from 1 − t³ must shrink as N grows, and does), the GPE solver, and Cooper/BCS against their closed forms. The GPE referee earned its keep the usual way — by catching a real bug:

1for iter in 0..200_000 {
2    // strang splitting: half kinetic, full potential+interaction, half kinetic
3    let mut q = psi.clone();
4    // …
5    // Renormalize BEFORE the nonlinear step — a war story. Imaginary time
6    // SHRINKS the norm, and feeding g|ψ|² a decayed ψ is an O(dt) error
7    // in the coupling. We didn't spot it by staring; the VIRIAL REFEREE
8    //     2 E_kin - 2 E_pot + E_int = 0   (exact at the true minimum)
9    // did: first version |virial| = 1.9e-4. Renormalizing here cut it to
10    // 2.5e-5 — still not zero, because the remainder is the splitting's
11    // own O(dt²) fixed-point bias, which no smaller dt could chase down.
12    // That bias is what the projected-gradient polish below removes.
13    normalize(&mut q);
14    for (c, &x) in q.iter_mut().zip(&xs) {
15        let v = 0.5 * x * x + g * c.norm_sqr();
16        *c *= (-dt * v).exp();
17    }
18    // …
19}

1// The cure for the splitting's O(dt²) fixed-point bias, which no smaller
2// dt could chase down: polish with projected gradient flow,
3//     ψ ← normalize(ψ − dτ (H_GP ψ − μψ)).
4// Its FIXED POINT satisfies H_GP ψ = μψ exactly, for ANY dτ — the step
5// size sets the speed, not the answer — so the virial and μ(x)-constancy
6// referees below end up testing physics rather than dt bias. What that
7// buys, on this run's scoreboard:
8//     worst |2Ek - 2Ep + Ei|  = 1.8e-12   (was 1.9e-4)
9//     worst μ(x) spread       = 8.6e-11   (μ constant across the cloud)
10//     |E(g=0) - 1/2|          = 1.1e-16   (the oscillator, rediscovered)
11let apply_h = |p: &Vec<Complex64>| -> Vec<Complex64> {
12// …
13let dtau = 2e-4;
14for _ in 0..60_000 {
15    let hp = apply_h(&psi);
16    let mu: f64 = hp.iter().zip(&psi).map(|(h, p)| (h * p.conj()).re).sum::<f64>() * dx;
17    let mut res_norm = 0.0f64;
18    for i in 0..n {
19        let r = hp[i] - psi[i] * mu;
20        res_norm += r.norm_sqr();
21        psi[i] -= r * dtau;
22    }
23    normalize(&mut psi);
24    if (res_norm * dx).sqrt() < 1e-12 {
25        break;
26    }
27}

1/// Cooper pair binding: solve the pair equation λ ∫₀¹ dx/(2x + e) = 1 for the
2/// binding energy e (units ω_D = 1). Closed form: e = 2/(e^{2/λ} − 1).
3///
4/// WHY THIS INTEGRAND IS HARD. The integrand 1/(2x + e) has a pole at
5/// x = −e/2, sitting just left of the interval. As λ shrinks the binding
6/// e → 0 and that pole crawls toward x = 0, so the integrand spikes to 1/e at
7/// the left endpoint and collapses over a width ~ e/2. A uniform grid
8/// (midpoint or Simpson) must scatter O(1/e) nodes into that spike to resolve
9/// it — the old first-order MIDPOINT rule needed M = 200_000 nodes just to
10/// reach ~6e-7, because midpoint error is only O(h²) and here f'' ~ 1/e³.
11///
12/// TANH–SINH (DOUBLE-EXPONENTIAL) QUADRATURE. We instead substitute
13///     x(t) = ½·(1 + tanh(½π·sinh t)),   t ∈ (−∞, ∞),
14/// so dx = ½·(½π·cosh t)/cosh²(½π·sinh t) dt and ∫₀¹ f dx = ∫ f(x(t)) x'(t) dt,
15/// then apply the plain trapezoid rule in t with step h. Two things happen:
16// … (the payoff, from the same doc comment: machine precision, ~121 nodes)
17fn cooper_binding(lambda: f64, n_half: usize) -> f64 {
18// …
19        let cn = cooper_binding(lambda, 60); // 2·60+1 = 121 tanh–sinh nodes
20// …
21    // tanh–sinh reaches ~2e-15 (machine precision) for λ ≥ 0.3; the 1e-12
22    // referee keeps a few decades of headroom for cross-platform rounding
23    // while still being a 6-orders-of-magnitude tightening over the old 1e-6
24    // midpoint bound. (Below λ ≈ 0.3 the binding e^{−1/λ} dives under any grid
25    // refinement — the chapter's moral — so the numeric check stays at λ ≥ 0.3.)
26    score(&mut rs, "cooper_worst", cooper_worst, 1e-12);
27    score(&mut rs, "bcs_gap_worst", bcs_worst, 1e-12);

1. (F) Derive (16.1) and (16.2) from scratch, then put in numbers for copper (one conduction electron per atom, n ≈ 8.5×10²⁸ m⁻³): E_F in eV, T_F in kelvin, the Fermi velocity, and the degeneracy pressure in pascals. Compare that pressure with atmospheric.
2. (C) In the statistics explorer, cool slowly from T/T_c = 1.3 to 0.3 and record N₀/N at five temperatures. Fit the exponent p in 1 − (T/T_c)^p and compare with the trap value p = 3 — then explain why a box of the same N would give p = 3/2 (count the density of states).
3. (F, hard) Chandrasekhar on one page. Model a white dwarf as N_e degenerate electrons at uniform density in radius R, plus gravity −αGM²/R. Non-relativistic pressure (E_kin ∝ N_e^(5/3)/R²) gives a stable minimum with R ∝ M^(−1/3) — heavier dwarfs are smaller. Now let the electrons go ultra-relativistic (ε = ħck): show E_kin ∝ N_e^(4/3)/R, the same power of R as gravity, so stability becomes a fight of coefficients — and extract the critical mass where support fails. You should land within a factor of two of 1.4 solar masses using only ħ, c, G, and the proton mass.
4. (F) The exchange discount, wholesale: for the uniform electron gas, the HF exchange energy per particle is −(3/4)(e²k_F/π)(1/4πε₀). Combine with (16.2) to get E/N = 2.21/r_s² − 0.916/r_s Rydberg (r_s = mean spacing in Bohr radii), minimize over r_s, and compare the predicted equilibrium density and cohesive energy with sodium's (r_s ≈ 4). You have computed, crudely, why metals hold together.
5. (F) Derive both distributions (16.4) from Chapter 15's matrix elements: balance emission ∝ (n+1) against absorption ∝ n for bosons (and 1−n against n for fermions) in contact with a Boltzmann reservoir, and solve for n̄(ε). Statistics is the operator algebra, thermalized.
6. (F) Prove the virial identity (16.8): scale ψ_λ(x) = √λ ψ(λx), write E(λ) term by term, and differentiate at λ = 1. Check your signs against the lab table — every row satisfies it to ~10⁻¹². Then derive the Thomas–Fermi chemical potential μ_TF = (3g/4√2)^(2/3) and verify the g = 50 row.
7. (F, hard) Cooper's calculation, in full. Two electrons added to a frozen Fermi sea, attraction −V in the shell E_F < ε < E_F + ħω_D; write the pair state Σ a_k c†_{k↑}c†_{−k↓}|sea⟩, turn the Schrödinger equation into 1 = VΣ_k 1/(2ξ_k − E), replace the sum by N(0)∫dξ, and solve to get (16.9). Then show (a) binding occurs for any V > 0, (b) d^n E_B/dλ^n → 0 at λ = 0⁺ for all n, and (c) in free space (no sea) the same equation in 3D has a binding threshold — identify exactly what the sea changes (the effective density of states at zero pair energy).
8. (F, hard) The BCS state |BCS⟩ = Π_k (u_k + v_k c†_{k↑}c†_{−k↓})|0⟩ with u_k² + v_k² = 1. Show it does not conserve particle number — and why that is the same honesty as Chapter 7's coherent states (sharp phase, fuzzy number). Minimize ⟨H − μN⟩ over v_k to derive the gap equation, verify Δ = ħω_D/sinh(1/λ) for the shell interaction, and confirm the quasiparticle spectrum √(ξ² + Δ²) — no excitation cheaper than Δ.
9. (P, hard) Give the lab's GPE solver a stress test with an exact answer: set the trap to zero, make the interaction attractive (g < 0), and converge in a large box. The 1D GPE has an exact normalized bright-soliton ground state $\psi (x)=\sqrt{-g/4}\mathrm{sech}(-gx/2)$ with μ = −g²/8 (verify both by substitution first). Referee your numerics against the profile and μ, show the trap virial (16.8) must be replaced in free space by 2E_kin + E_int = 0 (re-run the λ-scaling with no E_pot) and check the solver satisfies it, and report how the needed box size scales with 1/|g|.