Hermite–Gaussians in the position basis

The position basis is the one where x̂ acts by multiplication and p̂ by differentiation: ⟨x|x̂|ψ⟩ = xψ(x) and ⟨x|p̂|ψ⟩ = −iħ ψ′(x). Feed that into the dimensionless quadratures X = x̂/x₀ and P = x₀p̂/ħ, with x₀ = √(ħ/mω), and a = (X + iP)/√2 becomes (ξ + d/dξ)/√2 while a† becomes (ξ − d/dξ)/√2, where ξ = x/x₀. The sign flip in a† is the formal adjoint of d/dξ under integration by parts, and the boundary term you drop vanishes only because states in the domain of p̂ vanish at ±∞. On a half-line, or against a hard wall, the ladder fails even where ψ = 0 at the wall kills that term, because d/dξ carries a state obeying the wall condition into one that does not. Nothing physical has changed. |n⟩ is still the same vector — you have only chosen a continuum basis in which to write its components, and every commutator identity proved upstairs survives as an identity between differential operators.

Project a|0⟩ = 0 onto ⟨x| and you get (ξ + d/dξ)ψ₀ = 0. Separate it: dψ₀/ψ₀ = −ξ dξ, so ψ₀ = C exp(−ξ²/2). Because the equation is first order its solution space is one-dimensional — that single decaying Gaussian is the whole of it, and it is normalisable on the spot. Contrast the textbook route. Solving −(ħ²/2m)ψ″ + ½mω²x²ψ = Eψ by series gives a two-dimensional solution space at every E, and you quantise by demanding the series terminate so the runaway exp(+ξ²/2) branch dies. The ladder did that job already, which is precisely why nothing is left to throw away here. Normalising, ∫|ψ₀|²dx = |C|²x₀√π = 1, so ψ₀(x) = (mω/πħ)¹⁄⁴ exp(−x²/2x₀²). For ⁸⁷Rb in a 150 Hz trap, x₀ = 0.8805 µm and |ψ₀(0)|² = 1/(√π x₀) = 6.41 × 10⁵ m⁻¹. Note that x₀ is the only length in the problem: every other length you quote is a pure number times x₀.

Write ψₙ = fₙ(ξ)e^(−ξ²/2) and apply the raising operator: (ξ − d/dξ)[f e^(−ξ²/2)] = (2ξf − f′)e^(−ξ²/2). The Gaussian is untouched — the entire action falls on the polynomial, as the map f ↦ 2ξf − f′. Start from f = 1 and iterate: 1 → 2ξ → 4ξ² − 2 → 8ξ³ − 12ξ → 16ξ⁴ − 48ξ² + 12. Those are H₀ through H₄. Indeed Hₙ₊₁ = 2ξHₙ − Hₙ′ is one standard Hermite recursion, and substituting Hₙ′ = 2nHₙ₋₁ turns it into the three-term form Hₙ₊₁ = 2ξHₙ − 2nHₙ₋₁. Keeping the normalisation ψₙ₊₁ = (ξ − d/dξ)ψₙ/√(2(n+1)) supplies the prefactor (2ⁿn!)^(−1/2), giving ψₙ(x) = (2ⁿn!)^(−1/2)(mω/πħ)¹⁄⁴Hₙ(ξ)e^(−ξ²/2). Each application raises the polynomial degree by exactly one, which is the structural reason state n has n zeros and not some other count.

Parity Π commutes with H, and ΠaΠ = −a, so once Π|0⟩ = |0⟩ is fixed by the even Gaussian, Π|n⟩ = (−1)ⁿ|n⟩ and ψₙ(−x) = (−1)ⁿψₙ(x). Every odd state therefore vanishes exactly at the origin — a real dark line in an absorption image, not a bookkeeping accident. The degree-n polynomial has n real simple zeros, all inside |ξ| < √(2n+1): setting (n + ½)ħω = ½mω²xₜ² gives the classical turning point xₜ = √(2n+1)x₀, so every node lies in the classically allowed region and the forbidden region carries only exponential tail. For n = 2 with x₀ = 0.8805 µm the nodes are at ±0.623 µm, the outer density maxima at ±1.392 µm, and xₜ = ±1.969 µm; the outer lobes stand 1.31 times the central one — the first hint of the classical density that piles up where the particle moves slowest. The tail is not negligible either: for n = 0, P(|x| > xₜ) = erfc(1) = 0.157.

∫ψₘ*(x)ψₙ(x)dx = ∫⟨m|x⟩⟨x|n⟩dx = ⟨m|(∫|x⟩⟨x|dx)|n⟩ = ⟨m|n⟩ = δₘₙ. The intimidating Hermite orthogonality integral ∫HₘHₙ e^(−ξ²)dξ = √π 2ⁿn! δₘₙ is therefore a corollary of the ladder's ⟨m|n⟩ = δₘₙ, not an independent fact to look up. Completeness runs the same way: Σₙ|n⟩⟨n| = 1 becomes Σₙ ψₙ(x)ψₙ(x′) = δ(x − x′), so any square-integrable ψ expands with cₙ = ∫ψₙ(x)ψ(x)dx, real coefficients for real ψ. Two payoffs. Parity kills half the integrals before you start — an even ψ has cₙ = 0 for every odd n. And the argument repeats in momentum space: in the p basis a = i(η + d/dη)/√2 with η = p/p₀, the position-basis form times i, so ⟨p|n⟩ is the same Hermite–Gaussian with x₀ replaced by p₀ = ħ/x₀ = √(mħω), times (−i)ⁿ. These functions are eigenfunctions of the Fourier transform.

Do not code the closed form. At n = 151 the product 2ⁿn! already exceeds 10³⁰⁸ and overflows float64, and Hₙ(ξ) is itself of order √(2ⁿn!) or larger, while the normalised product (2ⁿn!)^(−1/2)Hₙ(ξ)e^(−ξ²/2) stays of order one. Divide the three-term recursion by the normalisation instead: ψₙ₊₁ = ξ√(2/(n+1))ψₙ − √(n/(n+1))ψₙ₋₁, started from ψ₀ = π^(−1/4)e^(−ξ²/2) and ψ₁ = √2 ξψ₀ in units where x₀ = 1. Everything stays O(1) at every step. Two traps worth checking. NumPy ships both families — numpy.polynomial.hermite is the physicists' Hₙ and numpy.polynomial.HermiteE is the probabilists' Heₙ(ξ) = 2^(−n/2)Hₙ(ξ/√2) — and mixing them costs a factor of √2 per degree. And a finite-difference Hamiltonian on a grid reproduces the low states well but distorts the high ones, whose width outgrows the box: test it by computing ∫ψₘψₙ dx with the trapezoid rule and counting sign changes against n.