University Physics V · The Quantum Harmonic Oscillator · 7.6
Hermite–Gaussians in the Position Basis
You already have the spectrum without solving a differential equation. Here the algebra is cashed out into a picture: one first-order equation hands you the Gaussian, repeated raising manufactures every polynomial above it, and nodes, parity and lobe heights become things you predict before plotting.
Build the model
Connect the measurement to the mechanism.
The Fock vectors |n⟩ are built from algebra alone — a, a† and the positivity of a norm — and nothing in that construction mentions a coordinate. To get a function of x you must choose a basis, and the position basis is the one where x̂ acts by multiplication and p̂ by −iħ d/dx. Making that substitution turns operator statements into differential ones, and the whole economy of this topic is that the statement you translate is a|0⟩ = 0, not H|ψ⟩ = E|ψ⟩.
A first-order equation replaces a second-order one: (ξ + d/dξ)ψ₀ = 0 has a one-dimensional solution space, C exp(−ξ²/2), already normalisable, so there is no runaway branch to discard and no quantisation condition to impose — quantisation was settled upstairs by the ladder. Every state above the ground state then comes from applying (ξ − d/dξ)/√(2(n+1)) again and again, and because that map sends a polynomial f to 2ξf − f′, it is the Hermite recursion wearing a disguise. What you buy is a complete orthonormal set of real functions with exactly n nodes, parity (−1)ⁿ, and one universal Gaussian factor set by x₀ = √(ħ/mω).
What you pay is a conceptual bill students keep failing to settle: ψₙ(x) is not a new object but the component list of |n⟩ along a continuum basis, and the same vector read along |p⟩ returns the same Hermite–Gaussian up to a phase (−i)ⁿ.
- Simple definition
- The Hermite–Gaussian function ψₙ(x) = ⟨x|n⟩ is the n-th oscillator eigenvector written out in the position basis: a degree-n Hermite polynomial in ξ = x/x₀ multiplying one universal Gaussian exp(−ξ²/2), with x₀ = √(ħ/mω).
- Example
- For ⁸⁷Rb in a 150 Hz trap, x₀ = 0.8805 µm, so ψ₂ ∝ (2ξ² − 1)exp(−ξ²/2) has nodes at x = ±x₀/√2 = ±0.623 µm and outer density peaks at ±√2.5 x₀ = ±1.392 µm, standing 1.31 times taller than the central one.
Turns every ladder identity into calculus on functions, without ever writing the Schrödinger equation down.
x₀ = √(ħ/mω) in metres, ξ dimensionless, d/dξ = x₀ d/dx; neither a nor a† is Hermitian.
First order, so the solution space is one-dimensional: there is no growing branch left to discard.
ψ₀ carries units m(−1/2); for ⁸⁷Rb at 150 Hz, x₀ = 0.8805 µm and |ψ₀(0)|² = 6.41 × 10⁵ m⁻¹.
Generates H₀ = 1, H₁ = 2ξ, H₂ = 4ξ² − 2, H₃ = 8ξ³ − 12ξ by hand, with no integral table anywhere.
f is the polynomial factor, dimensionless; the map raises its degree by exactly one.
One expression for all n, with the prefactor fixed by nothing more than ⟨n|n⟩ = 1.
Hₙ takes the dimensionless ξ, never x in metres — the commonest units slip in this topic.
The route NumPy code should take: iterate on normalised values, never on Hₙ and 2ⁿn! separately.
Already normalised, so 2ⁿn! never appears; that product overflows float64 at n = 151.
All n zeros of Hₙ lie inside |ξ| < ξₜ, so every node sits in the classically allowed region.
ξₜ is the classical turning point in units of x₀; for n = 4, xₜ = 3x₀ exactly.
Reading a and a† 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 square-integrable functions decay at ±∞; on a half-line, or against a hard wall, the same manipulation is illegal. 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.
One first-order equation, and no branch to discard
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₀.
Raising is the Hermite recursion in disguise
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_(n+1) = 2ξHₙ − Hₙ′ is one standard Hermite recursion, and substituting Hₙ′ = 2nH_(n−1) turns it into the three-term form H_(n+1) = 2ξHₙ − 2nH_(n−1). Keeping the normalisation ψ_(n+1) = (ξ − 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.
Nodes, parity, and the √(2n+1) turning point
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.
Orthonormality is inherited, never re-proved
∫ψₘ*(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: a is symmetric under X ↔ P up to a factor of 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.
Building them in NumPy without overflowing
Do not code the closed form. At n = 151 the product 2ⁿn! already exceeds 10³⁰⁸ and overflows float64, and Hₙ(ξ) is comparably enormous, while their ratio stays of order one. Divide the three-term recursion by the normalisation instead: ψ_(n+1) = ξ√(2/(n+1))ψₙ − √(n/(n+1))ψ_(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 hermiteₑ 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.
Change one variable at a time
Make the relationship visible.
Step n from 0 to 4 and watch the dashed Gaussian stay exactly where it is — every new node comes from Hn, never from a wider envelope. Then move x₀: the whole picture rescales and the readout shows ω/2π falling as 1/x₀².
NODES n2
PARITY (−1)n1
TURNING POINT xₜ2.01 µm
TRAP ω/2π (87-Rb)144 Hz
Live interpretationNODES n: 2. PARITY (−1)n: 1. TURNING POINT xₜ: 2.01 µm. TRAP ω/2π (87-Rb): 144 Hz
Catch the common trap
Explain before calculating.
An oscillator eigenstate |n⟩ with n = 4 is read in the position basis, with length scale x₀ = √(ħ/mω). Which description of ψ₄(x) = ⟨x|4⟩ is correct?
Choose an answer to test the model.
Practice & worked examples
Reason from the model, then test the result.
EasyA ⁸⁷Rb atom (m = 1.443 × 10⁻²⁵ kg) sits in a harmonic optical trap with ω/2π = 150 Hz. Project a|0⟩ = 0 onto ⟨x|, solve and normalise the resulting equation, then evaluate the length scale x₀ and the peak probability density |ψ₀(0)|².
- In the position basis x̂ → x and p̂ → −iħ d/dx, so a = (ξ + d/dξ)/√2 with ξ = x/x₀ and x₀ = √(ħ/mω). The operator statement ⟨x|a|0⟩ = 0 is then the first-order equation (ξ + d/dξ)ψ₀ = 0.
- Separate variables: dψ₀/ψ₀ = −ξ dξ, so ln ψ₀ = −ξ²/2 + const and ψ₀ = C e(−ξ²/2). The equation is first order, so this one function spans the entire solution space — there is no second, growing branch to reject.
- Normalise: ∫|ψ₀|²dx = |C|²x₀∫e(−ξ²)dξ = |C|²x₀√π = 1, so C = (πx₀²)(−1/4) and ψ₀(x) = (mω/πħ)¹⁄⁴ exp(−x²/2x₀²).
- Numbers: ω = 2π × 150 = 942.5 s⁻¹, so x₀ = √(1.0546 × 10⁻³⁴ / (1.443 × 10⁻²⁵ × 942.5)) = √(7.753 × 10⁻¹³ m²) = 8.805 × 10⁻⁷ m = 0.8805 µm.
- Peak density: |ψ₀(0)|² = C² = 1/(√π x₀) = 1/(1.7725 × 8.805 × 10⁻⁷ m) = 6.41 × 10⁵ m⁻¹, a density per unit length, as m(−1/2) squared demands.
Answerψ₀(x) = (mω/πħ)¹⁄⁴ exp(−x²/2x₀²) with x₀ = 0.8805 µm, and |ψ₀(0)|² = 6.41 × 10⁵ m⁻¹.
MediumFor the same 150 Hz ⁸⁷Rb trap (x₀ = 0.8805 µm), raise twice to build ψ₂ in the position basis. Locate its nodes and its density maxima in micrometres, compare the outer lobe height with the central one, and check that everything sits inside the classical turning points.
- Raising reads ψ_(n+1) = (ξ − d/dξ)ψₙ/√(2(n+1)). Acting on f(ξ)e(−ξ²/2) it returns (2ξf − f′)e(−ξ²/2)/√(2(n+1)), so only the polynomial ever changes.
- From f₀ = 1: f₁ = 2ξ/√2 = √2 ξ. Then f₂ = (2ξ⋅√2ξ − √2)/√4 = √2ξ² − √2/2 = (2ξ² − 1)/√2. That matches (2²⋅2!)(−1/2)H₂(ξ) with H₂ = 4ξ² − 2, since (4ξ² − 2)/√8 = (2ξ² − 1)/√2.
- Nodes: 2ξ² − 1 = 0 gives ξ = ±1/√2 = ±0.7071, so x = ±0.7071 × 0.8805 µm = ±0.623 µm. Two nodes, and ψ₂ is even, as parity (−1)² requires.
- Density extrema: d/dξ[(2ξ² − 1)²e(−ξ²)] = 2ξ(2ξ² − 1)(5 − 2ξ²)e(−ξ²). Setting it to zero and discarding ξ = 0 and the two nodes leaves ξ = ±√(5/2) = ±1.5811, i.e. x = ±1.392 µm.
- Lobe ratio: |ψ₂(±√2.5 x₀)|²/|ψ₂(0)|² = (2 × 2.5 − 1)²e(−2.5)/(−1)² = 16 × 0.08209 = 1.313, so the outer lobes stand 31% above the central one.
- Turning points: (2 + ½)ħω = ½mω²xₜ² gives xₜ = √5 x₀ = 2.2361 × 0.8805 µm = 1.969 µm. Both maxima at ±1.392 µm and both nodes lie well inside it; beyond it there is only exponential tail.
Answerψ₂ ∝ (2ξ² − 1)e(−ξ²/2); nodes at ±0.623 µm, density maxima at 0 and ±1.392 µm with the outer lobes 1.31× the central one, all inside xₜ = ±1.969 µm.
HardThe atom sits in the ground state of the 150 Hz trap when the trap frequency is suddenly cut to ω′ = ω/4, too fast for the state vector to change. Find the probability of finding the atom in each eigenstate of the new trap, and the mean energy immediately after the quench.
- Sudden means the ket is unchanged and only the basis you resolve it in changes, so P(n′) = |⟨n′|0⟩|². Since x₀ = √(ħ/mω), the new scale is x₀′ = √(ħ/m(ω/4)) = 2x₀ = 1.761 µm.
- Parity first: the old ψ₀ is even and ψₙ′ has parity (−1)ⁿ, so ⟨n′|0⟩ = 0 identically for every odd n′. Half the integrals never need doing.
- Ground-to-ground overlap of two Gaussians of widths a and b: ⟨0′|0⟩ = ∫(πab)(−1/2)e(−x²(a⁻² + b⁻²)/2)dx = √(2ab/(a² + b²)). With b = 2a this is √(4/5) = 0.8944, so P(0′) = 0.800.
- The even overlaps follow the squeezed-vacuum pattern P(2k) = P(0′)⋅[(2k)!/(2k k!)²]λ(2k) with λ = (b² − a²)/(b² + a²) = 3/5. So P(2′) = 0.800 × (2/4) × 0.36 = 0.144 and P(4′) = 0.800 × (24/64) × 0.1296 = 0.0389.
- Check completeness: 0.800 + 0.144 + 0.0389 = 0.983, and P(6′) = 0.800 × (720/2304) × 0.046656 = 0.0117 takes the running sum to 0.994, converging to 1 as Σ|n′⟩⟨n′| = 1 requires.
- Energy: ⟨H′⟩ = ⟨p²⟩/2m + ½m(ω/4)²⟨x²⟩ evaluated in the old ground state, where ⟨p²⟩/2m = ħω/4 and ⟨x²⟩ = x₀²/2 = ħ/2mω. That gives ħω/4 + ħω/64 = 17ħω/64, and with ω = 4ω′ this is 1.0625 ħω′.
AnswerP(0′) = 0.800, P(2′) = 0.144, P(4′) = 0.0389, every odd n′ exactly zero by parity; ⟨H′⟩ = 17ħω/64 = 1.0625 ħω′, above the new zero-point energy 0.5ħω′.