Orthogonality, normalisation and basis expansion

A physical state is a ray: |ψ⟩ and z|ψ⟩ describe the same system for any nonzero complex z. Demanding ⟨ψ|ψ⟩ = 1 spends that freedom down to a single unobservable global phase, and it is what licenses reading squared moduli as probabilities. The mechanics is one inner product: |φ⟩ = |1⟩ + 2i|2⟩ has ⟨φ|φ⟩ = 1⋅1 + (−2i)(2i) = 1 + 4 = 5 — conjugate in the first slot, or the crossed i's give −4 and an impossible negative norm — so the normalised state is (|1⟩ + 2i|2⟩)/√5. On L² the same step fixes dimensions: for ψ(x) = A sin(πx/L) on [0, L], ∫|A|² sin²(πx/L) dx = |A|²L/2 = 1 gives A = √(2/L), and ψ carries units of m⁻¹⁄² so that |ψ|² dx is a pure number. Choosing A real and positive is convention, not physics: √(2/L)⋅e^(iα) is exactly as normalised.

Expand |ψ⟩ = Σₙ cₙ|n⟩ and apply ⟨m|: linearity gives ⟨m|ψ⟩ = Σₙ cₙ⟨m|n⟩, and if ⟨m|n⟩ = δₘₙ the sum collapses to the single term cₘ. That collapse is the whole convenience of an orthonormal basis — each coefficient is extracted independently, by one integral or one dot product, with no reference to the others. Remove orthogonality and the extraction entangles: for an independent but non-orthogonal set (|vₙ⟩), the same move gives ⟨vₘ|ψ⟩ = Σₙ Gₘₙcₙ with Gₘₙ = ⟨vₘ|vₙ⟩, a linear system to solve. Try |v₁⟩ = |1⟩, |v₂⟩ = (|1⟩ + |2⟩)/√2 and |ψ⟩ = |1⟩: the raw projections are 1 and 1/√2, but solving with G = [[1, 1/√2], [1/√2, 1]] returns c = (1, 0) — the second vector contributes nothing, even though ψ has a healthy overlap with it. Overlap and coefficient only coincide when the basis is orthonormal.

Compute ⟨ψ|ψ⟩ from the expansion and both slots bring a sum: ⟨ψ|ψ⟩ = Σₘ Σₙ cₘ* cₙ ⟨m|n⟩, an N²-term double sum. Orthonormality deletes every off-diagonal term and leaves Parseval's identity, Σₙ|cₙ|² = ⟨ψ|ψ⟩ = 1. This is the bookkeeping behind the Born rule: for a non-degenerate eigenvalue aₙ the numbers P(aₙ) = |cₙ|² are non-negative and now provably sum to one, so they form a probability distribution over measurement outcomes, with ⟨A⟩ = Σₙ aₙ|cₙ|² following at once. Read in reverse it is a diagnostic: probabilities that refuse to sum to one mean the basis is not orthonormal, the state is not normalised, or terms are missing. And it degrades gracefully — for any orthonormal subset, Σ|cₙ|² ≤ 1 (Bessel's inequality), the deficit equalling the probability of finding the system outside the spanned subspace. That inequality becoming equality for every state is exactly what completeness means.

Hermitian operators hand you orthogonality only across distinct eigenvalues; inside a degenerate eigenspace, the eigenvectors that fall out of solving (A − a)v = 0 by hand are merely independent. Gram–Schmidt repairs them. Normalise the first: |e₁⟩ = |v₁⟩/‖v₁‖. From the second, subtract its component along e₁: |w₂⟩ = |v₂⟩ − ⟨e₁|v₂⟩|e₁⟩, so ⟨e₁|w₂⟩ = ⟨e₁|v₂⟩ − ⟨e₁|v₂⟩ = 0 by construction; then normalise. Each later vector is stripped of every direction already banked before it is admitted. A zero remainder ‖wₙ‖ = 0 is information, not failure: vₙ was dependent on its predecessors and drops. Because every |eₖ⟩ is a combination of the degenerate eigenvectors, the new basis never leaves the eigenspace — orthonormalisation costs nothing in eigenvalue terms. Run on 1, x, x² with the inner product ∫₋₁¹ f*g dx, the same algorithm manufactures the Legendre polynomials up to normalisation; the ordering matters, and each ordering yields a different, equally legitimate basis.

No amount of orthonormality certifies that a set spans the space. Coefficients against an incomplete orthonormal set are each individually correct — still ⟨n|ψ⟩ — but Σ cₙ|n⟩ rebuilds only the projection of ψ onto the span, and the missing probability 1 − Σ|cₙ|² vanishes without an error message. In practice completeness is imported as a theorem — the spectral theorem in finite dimensions, Sturm–Liouville theory for eigenfunctions on an interval — and what you actually check is the truncation. How fast that audit converges is set by the boundary conditions: the flat state ψ(x) = 1/√L on [0, L] is square-integrable but non-zero at walls where every sine eigenfunction vanishes, so its odd-n coefficients cₙ = 2√2/(nπ) decay only as 1/n. The ground state captures |c₁|² = 8/π² = 0.811, adding n = 3 reaches 0.901, and Σ 8/(n²π²) over odd n crawls to 1 with a Gibbs overshoot at each wall — convergence in the mean-square norm, never pointwise. A smooth state that respects the walls converges orders of magnitude faster, and parity halves the work either way: a state even about the well centre has no even-n coefficients at all.

Stack the vectors as columns of A; then Q, R = np.linalg.qr(A) returns in Q an orthonormal basis for the same span, with R recording the projections Gram–Schmidt would have subtracted. Use it instead of coding the textbook loop: classical Gram–Schmidt in floating point loses orthogonality in proportion to κ(A)²ε for nearly dependent columns, the modified variant to κ(A)ε, while Householder QR keeps ‖Q†Q − I‖ at machine level regardless. Two habits keep grid work honest. First, the inner product must carry the measure: ⟨f|g⟩ ≈ Σⱼ f*(xⱼ)g(xⱼ)Δx, so vectors orthonormal under NumPy's plain dot product represent functions mis-scaled by √Δx — build the weight in before comparing with analytic states. Second, audit rather than trust: np.max(np.abs(Q.conj().T @ Q − np.eye(k))) should sit near 10⁻¹⁵, and Parseval on a known state should return 1 to the same tolerance, before any spectrum computed in that basis is believed.