University Physics V · The Quantum Harmonic Oscillator · 7.5

The Fock Basis & Operators as Matrices

Once |n⟩ is generated by a†, every operator in this course becomes a band of √n's you can type straight into NumPy without evaluating a single Hermite integral. The price is the last row: truncate the ladder and the matrix stops obeying the algebra that built it.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

The ladder hands you a basis for free. Starting from the one vector killed by a, repeated raising generates |n⟩ = (a†)ⁿ|0⟩/√(n!), and because each rung is an eigenvector of the Hermitian N = a†a with a distinct eigenvalue, the set is automatically orthonormal — no Gram–Schmidt, no Gaussian integrals. In that basis every question about the oscillator becomes matrix bookkeeping: ⟨m|a|n⟩ = √n δ(m, n−1) puts a on the first superdiagonal, a† on the first subdiagonal, H = ħω(N + ½) on the diagonal, and x and p become tridiagonal combinations with zero diagonal.

Products inherit bands rather than filling in — x² and p² are pentadiagonal, x⁴ carries nine bands — and the width of the band is a selection rule you can read off by eye. What it costs is that the space is infinite and your array is not. Truncating to N states replaces the identity by a projector P, and (PAP)(PBP) is not P(AB)P: the commutator picks up a defect −N|N−1⟩⟨N−1|, and ⟨N−1|x²|N−1⟩ comes out N−1 instead of 2N−1.

The matrix is the operator only well below the top of the basis, which is a convergence test to run, not an assumption to make.

Simple definition: The Fock basis is the orthonormal set of number eigenvectors |0⟩, |1⟩, |2⟩, … generated by applying a† to the vacuum, in which every oscillator operator is stored as an array of matrix elements ⟨m|O|n⟩ rather than as a differential operator.
Example: In it a is strictly upper triangular with √1, √2, √3 on the first superdiagonal, so the 4×4 block of x/x_zpf has rows (0, 1, 0, 0), (1, 0, 1.414, 0), (0, 1.414, 0, 1.732), (0, 0, 1.732, 0).

Ladder action and its normsa|n⟩ = √n |n−1⟩, a†|n⟩ = √(n+1) |n+1⟩

The root comes from a norm: ‖a†|n⟩‖² = ⟨n|(a†a + 1)|n⟩ = n + 1.

n = 0, 1, 2, …; both coefficients are pure numbers; a|0⟩ = 0 closes the ladder from below

Fock state built from the vacuum|n⟩ = (a†)ⁿ |0⟩ / √(n!)

Normalisation with no integral: every Fock state is fixed by the vacuum and a power of a†.

√(n!) is the running product √1⋅√2⋅…⋅√n picked up by n successive raisings

Matrix elements and where the bands sit⟨m|a|n⟩ = √n δ(m, n−1) · ⟨m|a†|n⟩ = √(n+1) δ(m, n+1)

In NumPy this is one line: a = np.diag(np.√(np.arange(1, N)), 1).

row m, column n: a is the first superdiagonal, a† the first subdiagonal, both with zero diagonal

x and p as tridiagonal matricesx = x_zpf (a + a†), p = i √(mħω/2) (a† − a), x_zpf = √(ħ/2mω)

Zero diagonal means ⟨n|x|n⟩ = 0; one band means Δn = ±1, the dipole selection rule.

x_zpf = √(ħ/2mω) in m, the momentum scale √(mħω/2) in kg m s⁻¹; both are tridiagonal

Squares are pentadiagonalx² = x_zpf²(a² + a†² + 2a†a + 1) → ⟨n|x²|n⟩ = (2n + 1) x_zpf²

Delivers ⟨x²⟩ and ΔxΔp = (n + ½)ħ by counting, plus the Δn = ±2 overtone band.

bands at Δn = 0 and ±2, with ⟨n+2|x²|n⟩ = x_zpf²√((n+1)(n+2))

What an N-state truncation breaks[a, a†]_N = 1 − N |N−1⟩⟨N−1|

Why ⟨N−1|x²|N−1⟩ returns N − 1 instead of 2N − 1: never trust the top rungs.

N × N truncation; the defect sits only in the bottom-right corner and grows with N

The basis generates itself from a|0⟩ = 0

Start from the one vector annihilated by a and let a† do the rest. The coefficient is not a convention pulled from nowhere: ‖a†|n⟩‖² = ⟨n|aa†|n⟩ = ⟨n|(a†a + 1)|n⟩ = n + 1, so a†|n⟩ = √(n+1)|n+1⟩ once the phase is chosen real and positive. Likewise ‖a|n⟩‖² = ⟨n|a†a|n⟩ = n gives a|n⟩ = √n|n−1⟩, and at n = 0 that returns the zero vector, which is what closes the ladder from below. Iterating the raising rule accumulates the product √1⋅√2⋅…⋅√n, so (a†)ⁿ|0⟩ = √(n!)|n⟩ and |n⟩ = (a†)ⁿ|0⟩/√(n!). Orthogonality is then free: N = a†a is Hermitian and its eigenvalues are distinct integers, so ⟨m|n⟩ = δ(m, n) without a single Gaussian integral being done. What the algebra does not supply is that each rung is one-dimensional; non-degeneracy is a separate input, and it is what lets the single label n name a basis vector.

Matrix multiplication is the completeness relation

An operator becomes a matrix the moment you write Oₘₙ = ⟨m|O|n⟩, row m and column n. For a, the delta δ(m, n−1) kills everything except the first superdiagonal, whose entries are √1, √2, √3, …; a† is its transpose, on the first subdiagonal. In a four-state block a has rows (0, 1, 0, 0), (0, 0, √2, 0), (0, 0, 0, √3), (0, 0, 0, 0) — strictly upper triangular, hence nilpotent, so a⁴ = 0 there. That nilpotency is the truncation talking, not the operator. Meanwhile H = ħω(a†a + ½) is diagonal with entries ħω(n + ½), so the eigenvalue problem is already solved in this basis and nothing needs diagonalising. The reason products of matrices compute products of operators is completeness: ⟨m|AB|n⟩ = Σₖ ⟨m|A|k⟩⟨k|B|n⟩ is the insertion of 1 = Σₖ |k⟩⟨k| written as a row-times-column sum. Hold on to that identity — it is exactly the step a truncation breaks.

x and p are one band wide, and that is a selection rule

Inverting a = (X + iP)/√2 gives x = x_zpf(a + a†) and p = i√(mħω/2)(a† − a), with x_zpf = √(ħ/2mω). Both are tridiagonal with a zero diagonal: ⟨n|x|n⟩ = 0 in every Fock state, so no stationary state of the oscillator sits off centre, and the only surviving elements are ⟨n+1|x|n⟩ = √(n+1) x_zpf and ⟨n−1|x|n⟩ = √n x_zpf. One band means Δn = ±1, which is the dipole selection rule: a strictly harmonic mode absorbs at ω and nowhere else. Put numbers on it. The CO stretch at 2143 cm⁻¹ has ω = 4.04 × 10¹⁴ rad s⁻¹ and reduced mass μ = 1.138 × 10⁻²⁶ kg, so x_zpf = 3.39 pm — about 3% of the 113 pm bond. Then ⟨1|x|0⟩ = 3.39 pm and ⟨2|x|1⟩ = √2 × 3.39 = 4.79 pm, while ⟨2|x|0⟩ = 0 exactly, which is why the overtone stays dark until anharmonicity or a quadratic dipole term supplies it.

Products widen the band; sparsity is what you buy

Squaring keeps the structure but doubles the reach. x² = x_zpf²(a² + a†² + 2a†a + 1) is pentadiagonal: a diagonal ⟨n|x²|n⟩ = x_zpf²(2n + 1), which hands you ⟨x²⟩ = (ħ/mω)(n + ½) and, with the matching p² result, the uncertainty product ΔxΔp = (n + ½)ħ by counting rather than integrating; plus a Δn = ±2 band with entries x_zpf²√((n+1)(n+2)), the overtone Raman scattering sees. Four factors reach further still: x⁴ has nine bands, at Δn = 0, ±2, ±4. Bandwidth is the point. An anharmonic Hamiltonian H = ħω(a†a + ½) + λx⁴ is banded, not full, so storing it costs O(N) and scipy.sparse.linalg.eigsh can extract the lowest levels without touching the rest. But banded is still not diagonal: with ε = λx_zpf⁴/ħω the first-order spacing becomes ħω[1 + 12ε(n + 1)], a drift no ħω(N + ½) can imitate.

Truncation: where the matrix stops being the operator

Your array has N rows; the space does not. Keeping n = 0 … N−1 replaces the identity by the projector P = Σ |n⟩⟨n| over kept rungs, and multiplying two truncated matrices evaluates PAPBP, not P(AB)P. The difference is A(1−P)B, which bites whenever A can climb out of the retained space and B can bring it back. Two symptoms, both exact and both checkable. The commutator loses a whole rung in one corner: [a, a†] built from N × N matrices equals 1 − N|N−1⟩⟨N−1|, so its last diagonal entry is −(N−1), not +1 — at N = 4 the diagonal reads (1, 1, 1, −3). And ⟨N−1|x²|N−1⟩ comes out N−1 where the operator gives 2N−1; at N = 4 that is 3 against 7, low by 57%. The counting rule is simple: a product of k ladder factors returning to level n passes through intermediates as high as n + ⌈k/2⌉, so the entry is exact only when N ≥ n + ⌈k/2⌉ + 1. For x² that reads N ≥ n + 2; for x⁴, N ≥ n + 3.

Building it in NumPy, and the tests it must pass

Four lines: a = np.diag(np.√(np.arange(1, N)), 1), then ad = a.T, x = xzpf*(a + ad), H = ℏ*w*(ad @ a + 0.5*np.eye(N)). Everything else is a matrix product. Then test rather than assume. First, np.allclose(a @ ad − ad @ a, np.eye(N)) will fail, and it should: check instead that the identity holds on the leading (N−1)×(N−1) block and that the last entry is −(N−1). Second, np.linalg.eigvalsh(H) must return ħω(n + ½) to machine precision, since H is already diagonal — that only checks your indexing, nothing physical. Third, the test that matters: double N and keep only the levels whose eigenvalues move by less than your tolerance; for a quartic term at ε = 0.002 the lowest six settle at N of a few tens, not at N = 8. Finally, cross-check against a real-space grid diagonalisation. Both methods degrade high on the ladder but for unrelated reasons — a finite box and a coarse mesh there, a missing top rung here — so agreement between them is evidence, while agreement of one method with itself is not.

Change one variable at a time

Make the relationship visible.

Interactive model

Fock states kept N4 states

Level read off n3

Set N = 4 and read n = 3: the bar stops at 3 where the operator says 7, because the rung to |4⟩ that x² needs is not in the space. Every bar below it is exact. Raise N and the short bar just moves to the new top row — truncation error never spreads downward.

EXACT 2n + 17 x²zpf

N×N MATRIX PRODUCT3 x²zpf

SHORTFALL57 %

TOP EXACT LEVEL2

Live interpretationEXACT 2n + 1: 7 x²zpf. N×N MATRIX PRODUCT: 3 x²zpf. SHORTFALL: 57 %. TOP EXACT LEVEL: 2

Catch the common trap

Explain before calculating.

You keep four Fock states (n = 0…3) and build X = a + a† as a 4×4 matrix in units of x_zpf, then form X² by matrix multiplication. Set against the operator's own ⟨n|x²|n⟩ = (2n + 1) x_zpf², what does the diagonal of that product read?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyUsing only a|n⟩ = √n |n−1⟩ and a†|n⟩ = √(n+1) |n+1⟩, evaluate a†a|2⟩ and aa†|2⟩, confirm that [a, a†]|2⟩ = |2⟩, and find the matrix element ⟨4|(a†)²|2⟩.

Lower first: a|2⟩ = √2 |1⟩. Raise that: a†(√2|1⟩) = √2 · √2 |2⟩ = 2|2⟩. So a†a|2⟩ = 2|2⟩ — the number operator reads the label back, as N|n⟩ = n|n⟩ requires.
Now the other order: a†|2⟩ = √3 |3⟩, and a(√3|3⟩) = √3 · √3 |2⟩ = 3|2⟩. So aa†|2⟩ = 3|2⟩.
Subtract: (aa† − a†a)|2⟩ = (3 − 2)|2⟩ = |2⟩. The commutator acts as the identity on this vector, as it must on every |n⟩ in the untruncated space.
For the matrix element, raise twice: (a†)²|2⟩ = a†(√3|3⟩) = √3 · √4 |4⟩ = √12 |4⟩. Projecting onto ⟨4| gives ⟨4|(a†)²|2⟩ = √12 = 2√3 = 3.4641.
Cross-check against the general rule ⟨n+r|(a†)ʳ|n⟩ = √((n+r)!/n!): with n = 2 and r = 2 that is √(4!/2!) = √12, the same number.

Answera†a|2⟩ = 2|2⟩, aa†|2⟩ = 3|2⟩, [a, a†]|2⟩ = |2⟩, and ⟨4|(a†)²|2⟩ = 2√3 ≈ 3.464. Every number in the Fock basis is a product of square roots of integers.

MediumKeep four Fock states (n = 0…3) and write X = x/x_zpf = a + a† as a 4×4 matrix. Form X² by matrix multiplication, then compare its diagonal and its Δn = 2 band with the exact ⟨n|x²|n⟩ = (2n + 1) x_zpf² and ⟨n+2|x²|n⟩ = √((n+1)(n+2)) x_zpf².

Only ⟨n+1|X|n⟩ = √(n+1) and its transpose are non-zero, so the three superdiagonal entries are √1 = 1, √2 = 1.4142 and √3 = 1.7321, and X is symmetric with a zero diagonal.
Rows of X: (0, 1, 0, 0), (1, 0, 1.4142, 0), (0, 1.4142, 0, 1.7321), (0, 0, 1.7321, 0).
Diagonal of the product, (X²)ₙₙ = Σₖ Xₙₖ Xₖₙ. n = 0: 1² = 1. n = 1: 1² + 1.4142² = 3. n = 2: 1.4142² + 1.7321² = 2 + 3 = 5. n = 3: 1.7321² = 3 alone, because the entry X₃₄ does not exist in a 4×4 array.
The operator's values are 2n + 1 = 1, 3, 5, 7. Three agree exactly; the top one reads 3 against 7, low by 4/7 = 57%.
The Δn = 2 band survives untouched: (X²)₂₀ = X₂₁ X₁₀ = 1.4142 × 1 = √2, and (X²)₃₁ = X₃₂ X₂₁ = 1.7321 × 1.4142 = 2.4495 = √6, both matching √((n+1)(n+2)) at n = 0 and n = 1.
The rule behind it: a product of k ladder factors reaches intermediates as high as n + ⌈k/2⌉, so ⟨n|x²|n⟩ needs N ≥ n + 2. At N = 4 that permits n ≤ 2 — exactly what the arithmetic shows.

Answerdiag(X²) = (1, 3, 5, 3) against the operator's (1, 3, 5, 7): only the top rung is wrong, and by 57%. The Δn = 2 entries √2 and √6 come out exact.

HardA mode is modelled as H = ħω(a†a + ½) + λx⁴ with x = x_zpf(a + a†), and ε = λx_zpf⁴/ħω = 0.002. Name the bands of x⁴, give ⟨n|x⁴|n⟩, then find the first-order E₀ and E₅, the spacings E₁ − E₀ and E₅ − E₄, and the smallest truncation N for which ⟨5|x⁴|5⟩ comes out exact.

x⁴ = x_zpf⁴(a + a†)⁴ is four ladder factors, so the reachable changes are Δn = 0, ±2, ±4 — nine bands, each entry a sum of products of square roots.
Collecting terms with equal numbers of a and a† gives ⟨n|x⁴|n⟩ = x_zpf⁴(6n² + 6n + 3). At n = 0 that is 3x_zpf⁴, three times the ⟨0|x²|0⟩² = x_zpf⁴ a naive factorisation would give — the Gaussian moment relation ⟨x⁴⟩ = 3⟨x²⟩².
First order: Eₙ = ħω[(n + ½) + ε(6n² + 6n + 3)]. For n = 0 the second term is 3, so E₀ = ħω(0.5 + 0.006) = 0.506 ħω. For n = 5 it is 150 + 30 + 3 = 183, so E₅ = ħω(5.5 + 0.366) = 5.866 ħω.
The neighbours: n = 1 gives 15 and n = 4 gives 123, so E₁ = ħω(1.5 + 0.030) = 1.530 ħω and E₄ = ħω(4.5 + 0.246) = 4.746 ħω.
Spacings: E₁ − E₀ = 1.024 ħω while E₅ − E₄ = 1.120 ħω, matching the general result ħω[1 + 12ε(n + 1)]. The ladder is no longer evenly spaced — the signature no ħω(N + ½) can fake.
Truncation: four ±1 steps returning to n = 5 must pass through intermediates as high as 5 + 2 = 7, so the basis has to contain |7⟩, i.e. N ≥ 8. At N = 7 the entry silently returns 141 instead of 183, and the 12% spacing shift would be misread as a smaller λ.

AnswerE₀ = 0.506 ħω and E₅ = 5.866 ħω; E₁ − E₀ = 1.024 ħω against E₅ − E₄ = 1.120 ħω, the spacing growing as 1 + 12ε(n + 1). ⟨5|x⁴|5⟩ needs N ≥ 8.