The Fock basis and operators as matrices

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: 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.

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.

Inverting a = √(mω/2ħ)(x + ip/mω) 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.

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 route by which a dipole or polarisability quadratic in x reaches the Δn = 2 overtone. Four factors reach further still: x⁴ spans nine diagonals, of which only the five at Δn = 0, ±2, ±4 are non-zero. 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 ħω(a†a + ½) can imitate.

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 PA(1−P)BP, which bites whenever B can carry a kept state out of the retained space and A 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)x_zpf² where the operator gives (2N−1)x_zpf²; at N = 4 that is 3x_zpf² against 7x_zpf², 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.

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.