Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section you will be able to:

Set up the time-independent Schrödinger equation for a particle in a parabolic well and explain why the parabola appears again and again in physics.
Derive the equally-spaced energy ladder $E_n = \hbar\omega(n + \tfrac{1}{2})$ using the ladder operators $\hat{a}, \hat{a}^\dagger$ .
Recognise the Hermite-Gaussian wave functions $\psi_n(x)$ and read off the number of nodes directly from $n$ .
Interpret the zero-point energy $E_0 = \tfrac{1}{2}\hbar\omega$ as a mandatory consequence of the uncertainty principle.
Compare the quantum probability density $|\psi_n|^2$ to the classical one and see the correspondence principle at work as $n \to \infty$ .
Build a coherent state in your head and explain why it is "the most classical" quantum state.
Implement the QHO twice in Python: once analytically with Hermite polynomials, once by diagonalising the discretised Hamiltonian in PyTorch.

The Big Picture

Almost any smooth potential, expanded around a stable minimum, looks like a parabola. That is the secret reason the harmonic oscillator keeps reappearing — in the rumbling of atoms in a crystal (phonons), in the rungs of a laser cavity (photons), in molecular vibrations, in every quantum field. Solve the QHO once and you have the universal building block for "small oscillations of anything."

Why this section matters: the QHO is the single most useful exactly-solvable problem in quantum mechanics. Whenever you meet a new quantum system — a vibrating molecule, a trapped atom, a mode of the electromagnetic field — your first move is to ask does this reduce to a harmonic oscillator near equilibrium? Nine times out of ten the answer is yes.

Three big surprises lie ahead, all absent from the classical version:

Energy is quantised. Not just "small," but discretely indexed by $n = 0, 1, 2, \ldots$ .
The ground state is not at rest. Even the lowest energy $E_0 = \tfrac{1}{2}\hbar\omega$ is positive — perfect stillness is forbidden.
Probability can leak past the classical turning point. The wavefunction is non-zero in regions where a classical marble simply could not reach.

The Classical Oscillator: Spring, Pendulum, Atom

A mass $m$ connected to a spring of stiffness $k$ obeys Newton's law $m\ddot{x} = -k x$ . Solutions oscillate at the angular frequency $\omega = \sqrt{k/m}$ and the total energy

E = \tfrac{1}{2} m \dot{x}^2 + \tfrac{1}{2} m \omega^2 x^2

is conserved. The two terms split smoothly into kinetic and potential energy; the second is the parabolic potential $V(x) = \tfrac{1}{2} m \omega^2 x^2$ .

Why the parabola is everywhere: Taylor-expand any smooth

V(x)

around a minimum

x_0

V(x) \approx V(x_0) + \tfrac{1}{2}V''(x_0)(x-x_0)^2

. The linear term vanishes because the minimum has zero slope; the quadratic term is exactly the spring energy with

k = V''(x_0)

. Every small oscillation in nature reduces to this same shape.

Classically, a marble released at amplitude $A$ swings between $+A$ and $-A$ , the two turning points where all energy is potential and the velocity is zero. Quantum mechanics is going to keep these turning points but change everything else.

The Quantum Hamiltonian

Replace position and momentum by operators — $\hat{x}$ acts by multiplication, $\hat{p} = -i\hbar\,d/dx$ — and the classical energy becomes the Hamiltonian operator

\hat{H} = \frac{\hat{p}^2}{2m} + \tfrac{1}{2}m\omega^2 \hat{x}^2 = -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + \tfrac{1}{2}m\omega^2 x^2

and the time-independent Schrödinger equation reads $\hat{H}\psi(x) = E\,\psi(x)$ . We are looking for the special functions $\psi$ that the Hamiltonian merely multiplies by a number — the energy eigenstates.

Why the parabola is so hard analytically: the equation

-\tfrac{\hbar^2}{2m}\psi'' + \tfrac{1}{2}m\omega^2 x^2 \psi = E\psi

mixes a constant-coefficient second derivative with a coefficient that grows like

x^2

. Brute-force series solutions work but generate ugly recurrence relations. The ladder operator trick below is dramatically cleaner — and gives the spectrum without solving a single differential equation explicitly.

The Parabolic Well and Its Quantized Ladder

The potential is V(x) = ½ω²x². Inside it sit equally spaced energy rungs Eₙ = ℏω(n + ½). The red dots mark the classical turning points — where a marble with that energy would stop and turn around.

Angular frequency ω1.00

Levels shown (n_max)6

Notice the rungs are equally spaced by ℏω — a signature only the parabolic well has.

The Ladder Operators — Quantum Mechanics' Cleanest Trick

Define two new operators by mixing position and momentum with carefully chosen prefactors:

\hat{a} = \sqrt{\frac{m\omega}{2\hbar}}\left(\hat{x} + \frac{i}{m\omega}\hat{p}\right), \qquad \hat{a}^\dagger = \sqrt{\frac{m\omega}{2\hbar}}\left(\hat{x} - \frac{i}{m\omega}\hat{p}\right).

The first lowers a quantum number by 1; the second raises it. The single fact you need about them is the commutator

[\hat{a}, \hat{a}^\dagger] \;=\; \hat{a}\hat{a}^\dagger - \hat{a}^\dagger \hat{a} \;=\; 1,

which follows from the canonical commutator $[\hat{x}, \hat{p}] = i\hbar$ . Rewriting the Hamiltonian in terms of these operators gives the punchline

\boxed{\,\hat{H} \;=\; \hbar\omega\left(\hat{a}^\dagger \hat{a} + \tfrac{1}{2}\right)\,}

and the operator $\hat{N} = \hat{a}^\dagger \hat{a}$ is called the number operator. Its eigenvalues are non-negative integers $n = 0, 1, 2, \ldots$ , and the Hamiltonian is just $\hbar\omega(\hat{N} + 1/2)$ .

Intuition: think of $\hat{a}^\dagger$ as a creation operator that adds one quantum of energy $\hbar\omega$ (a phonon, a photon — pick your favourite name), and $\hat{a}$ as the annihilation operator that removes one. The integer $n$ simply counts how many quanta are in the oscillator.

Energy Quantization Falls Out Automatically

Suppose $|n\rangle$ is an eigenstate of $\hat{N}$ with eigenvalue $n$ . A two-line calculation using $[\hat{N}, \hat{a}^\dagger] = \hat{a}^\dagger$ shows that $\hat{a}^\dagger|n\rangle$ is an eigenstate with eigenvalue $n+1$ , and $\hat{a}|n\rangle$ is one with eigenvalue $n-1$ . So once you find any eigenstate, you can climb up and down a discrete ladder.

The ladder must terminate at the bottom — otherwise we would generate states with negative $n$ , which would have negative $\langle\hat{a}^\dagger\hat{a}\rangle$ (impossible, it is a norm). The cleanest way to stop is the rule

\hat{a}\,|0\rangle = 0

— the ground state is annihilated by $\hat{a}$ . Acting with $\hat{a}^\dagger$ repeatedly builds every other state: $|n\rangle = \frac{1}{\sqrt{n!}}(\hat{a}^\dagger)^n|0\rangle$ . Energies follow immediately:

E_n = \hbar\omega\left(n + \tfrac{1}{2}\right), \qquad n = 0, 1, 2, \ldots

n	Eₙ (in units of ℏω)	Number of nodes	Symmetry
0	1/2	0	Even
1	3/2	1	Odd
2	5/2	2	Even
3	7/2	3	Odd
n	n + 1/2	n	(−1)ⁿ

Equal spacing is uniquely parabolic. The infinite square well gives

E_n \propto n^2

; the hydrogen atom gives

E_n \propto -1/n^2

. Only the parabolic well produces rungs separated by the same step

\hbar\omega

. That is why a perfect harmonic oscillator can absorb and emit light at exactly one frequency —

\omega

— over and over.

The Hermite Wavefunctions

Solving $\hat{a}|0\rangle = 0$ in the position basis gives a first-order differential equation $(\xi + d/d\xi)\psi_0(\xi) = 0$ where $\xi = x\sqrt{m\omega/\hbar}$ . The unique normalisable solution is the Gaussian

\psi_0(x) = \left(\frac{m\omega}{\pi\hbar}\right)^{1/4} e^{-m\omega x^2 / (2\hbar)}.

Applying $\hat{a}^\dagger$ repeatedly produces the full family

\psi_n(x) = \frac{1}{\sqrt{2^n n!}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/4} H_n(\xi)\, e^{-\xi^2/2}

where $H_n$ is the physicists' Hermite polynomial. The first few are

H_0 = 1, \quad H_1 = 2\xi, \quad H_2 = 4\xi^2 - 2, \quad H_3 = 8\xi^3 - 12\xi

and they satisfy the simple recurrence $H_{n+1}(\xi) = 2\xi H_n(\xi) - 2n H_{n-1}(\xi)$ — exactly what we use in code below.

Three features show up immediately in the picture:

Number of nodes equals $n$ (the polynomial degree). Each rung of the ladder adds one zero.
Parity alternates: even $n$ → even function; odd $n$ → odd function. Comes from the parity of the parabolic potential.
Gaussian envelope $e^{-\xi^2/2}$ kills the wavefunction in the classically forbidden region — but only exponentially, not abruptly, leaving a small leakage tail.

The Hermite Wave Functions — ψₙ(x) and |ψₙ(x)|²

Drag the n slider to climb the ladder. Each step adds exactly one node (yellow dot). In |ψ|² mode, compare the quantum probability to the classical one (orange dashed) — they only resemble each other for large n. This is the correspondence principle in action.

Quantum number nn = 3, E = 3.5 ℏω

Zero-Point Energy — The Lowest Rung Is Not Zero

At absolute zero a classical spring sits motionless at $x = 0$ with $E = 0$ . The quantum oscillator cannot — the rules of the game forbid it. To see why, use the uncertainty principle $\Delta x \,\Delta p \geq \hbar/2$ . For the ground state, the expectation value of the energy is

\langle \hat{H}\rangle = \frac{\langle \hat{p}^2\rangle}{2m} + \tfrac{1}{2}m\omega^2 \langle \hat{x}^2\rangle \;\geq\; \frac{(\Delta p)^2}{2m} + \tfrac{1}{2}m\omega^2(\Delta x)^2 \;\geq\; \frac{\hbar^2/4}{2m(\Delta x)^2} + \tfrac{1}{2}m\omega^2(\Delta x)^2.

Minimise the right-hand side over $\Delta x$ (take a derivative, set to zero) and you find the minimum value is exactly $\tfrac{1}{2}\hbar\omega$ . The uncertainty principle forces the spring to keep wobbling even at absolute zero.

Observable consequence: liquid helium remains liquid all the way down to $T = 0\,\mathrm{K}$ at atmospheric pressure. The zero-point motion of the He atoms is large enough to overcome the weak van-der-Waals attraction that would otherwise lock them into a crystal.

Quantum vs Classical Probability

Classically, a marble oscillating with amplitude $A$ spends most of its time near the turning points where it is moving slowly. The classical probability density is

P_{\mathrm{cl}}(x) = \frac{1}{\pi\sqrt{A^2 - x^2}}, \qquad |x| < A,

— a U-shape peaked at $\pm A$ . The quantum ground state, by contrast, is a Gaussian peaked at $x = 0$ : the opposite shape. As $n$ increases, however, the rapid oscillations of $|\psi_n|^2$ — when smoothed — bend up at the edges and start matching the classical U.

Correspondence principle: in the limit

n \to \infty

, the quantum probability (averaged over the fine oscillations) recovers the classical prediction. Toggle the Classical PDF overlay on the explorer above and push

n

up — by

n=10

the two are already telling the same story.

Coherent States — The Most Classical Quantum State

Energy eigenstates have probability densities that are stationary — they do not move. That seems unlike a classical oscillator, which definitely sloshes back and forth. The trick is to consider a special superposition of eigenstates, the coherent state:

|\alpha\rangle = e^{-|\alpha|^2/2}\sum_{n=0}^{\infty}\frac{\alpha^n}{\sqrt{n!}}|n\rangle.

This state is the unique eigenstate of the lowering operator: $\hat{a}|\alpha\rangle = \alpha|\alpha\rangle$ . Its position-space wavefunction is a Gaussian of the same width as the ground state, but centred at $x_c = \sqrt{2\hbar/(m\omega)}\,\mathrm{Re}(\alpha)$ . Under time evolution the center oscillates classically,

\langle x(t)\rangle = \sqrt{\tfrac{2\hbar}{m\omega}}\,|\alpha|\cos(\omega t - \varphi),

while the shape of $|\psi|^2$ never changes — the wave packet does not spread, in stark contrast to the free-particle packet of section 29.3.

Coherent State — A Quantum Wave Packet That Oscillates Like a Marble

A coherent state is the closest thing quantum mechanics has to a classical oscillating particle. Its shape never changes; only its center swings back and forth at frequency ω, just like a mass on a spring. Period T = 2π/ω ≈ 6.28.

Coherent amplitude |α|2.50

Angular frequency ω1.00

⟨n⟩ = |α|² = 6.25 (mean photon / phonon number) · ⟨H⟩ = ℏω(|α|² + ½) = 6.75 ℏω · current t = 0.00

Why lasers care: a single mode of a laser cavity is a harmonic oscillator (the electromagnetic field oscillates), and the light coming out is well-described by a coherent state with very large

|\alpha|^2

. Its photon number

\langle n\rangle = |\alpha|^2

can be in the trillions, yet its quantum fluctuations behave exactly like classical noise around a sinusoid. Coherent states are the bridge between the photon picture and the classical wave picture of light.

Worked Example: A Real Carbon-Monoxide Molecule

Let's plug real numbers into the QHO. The C-O bond in carbon monoxide vibrates at wavenumber $\tilde{\nu} \approx 2143\,\mathrm{cm}^{-1}$ (measured by infrared spectroscopy). That fully determines the energy ladder — no extra inputs needed.

Click to work through the numbers by hand

Step 1 — convert wavenumber to angular frequency. By definition $\omega = 2\pi c \tilde{\nu}$ . With $c = 2.998 \times 10^{8}\,\mathrm{m/s}$ and $\tilde{\nu} = 2.143 \times 10^{5}\,\mathrm{m^{-1}}$ :

\omega = 2\pi \cdot (2.998\times 10^{8}) \cdot (2.143\times 10^{5}) \approx 4.04 \times 10^{14}\,\mathrm{rad/s}.

Step 2 — compute one quantum of energy. $\hbar\omega = (1.0546\times 10^{-34})(4.04\times 10^{14}) \approx 4.26 \times 10^{-20}\,\mathrm{J}$ . Converting to electron-volts (divide by $1.602\times 10^{-19}$ ) gives $\hbar\omega \approx 0.266\,\mathrm{eV}$ .

Step 3 — build the ladder. $E_0 = 0.133\,\mathrm{eV}$ , $E_1 = 0.399\,\mathrm{eV}$ , $E_2 = 0.665\,\mathrm{eV}$ , and so on. Each rung is exactly $0.266\,\mathrm{eV}$ above the previous.

Step 4 — predict the absorption wavelength. A photon that bumps the molecule from $n=0$ to $n=1$ must carry energy $\Delta E = \hbar\omega$ , whose wavelength is

\lambda = \frac{hc}{\hbar\omega} = \frac{2\pi c}{\omega} \approx 4.67\,\mu\mathrm{m}.

That puts the line squarely in the mid-infrared — and indeed every astrophysicist looking at cold molecular clouds measures CO via this 4.67 μm transition. The QHO prediction agrees with experiment to about four significant figures.

Step 5 — sanity-check zero-point motion. Convert to position units: $\Delta x_0 = \sqrt{\hbar/(2 m_\mathrm{red} \omega)}$ with the reduced mass $m_\mathrm{red} = 1.139 \times 10^{-26}\,\mathrm{kg}$ gives $\Delta x_0 \approx 3.4\,\mathrm{pm}$ , roughly 3% of the equilibrium bond length of 113 pm. Tiny, but non-zero — the bond is permanently jittering at absolute zero.

Python Implementation: Hermite Wave Functions

Let's now build the QHO from scratch with a few lines of NumPy. The plan: define $E_n$ and $\psi_n(x)$ in code, check that each wavefunction is properly normalised, then draw them on the parabolic well to match the picture you have been looking at.

Hermite Wave Functions — Pure NumPy

🐍python

Explanation(13)

Code(43)

1Standard scientific Python imports

We use `math.factorial` for n!, NumPy for fast array math, and `hermval` from `numpy.polynomial.hermite` to evaluate the physicists' Hermite polynomial Hₙ(x) reliably at any x. `matplotlib` plots the result.

8Energy ladder in natural units

Choosing units where ℏ = m = ω = 1 strips every formula down to its mathematical bones. The full SI formula is Eₙ = ℏω(n+½), but it is exactly the function `energy(n) = n + 0.5` in our units. The reader can restore physical units by multiplying by ℏω at the end.

13Build the wavefunction ψₙ(x)

We construct ψₙ from three pieces: (i) Hermite polynomial Hₙ(x) — supplies the oscillations and the right number of nodes; (ii) Gaussian envelope exp(-x²/2) — pins the wavefunction near the well; (iii) normalization 1/√(2ⁿ n!) · π^(-1/4) — makes ∫|ψ|²dx = 1.

15Selecting a single Hermite polynomial

`hermval(x, coeff)` evaluates the polynomial c₀H₀(x) + c₁H₁(x) + … . By zeroing every coefficient except the n-th, we extract Hₙ alone. For n = 3 the coeff array is [0, 0, 0, 1] and hermval returns 8x³ − 12x at each x.

16Hermite evaluation, vectorised

`x` is an array of 4001 sample points; `hermval` returns an array of the same shape with Hₙ evaluated at each one in pure C code — no Python loop. For ψ₃ at x = 1: H₃(1) = 8·1 − 12·1 = −4. The Gaussian exp(-0.5) ≈ 0.6065 dampens it. Result before normalization: -2.426.

17Normalization constant

For n = 3: 1/√(2³·6) · π^(−1/4) = 1/√48 · 0.7511 ≈ 0.1083. Multiplying our preliminary −2.426 by this gives ψ₃(1) ≈ −0.2628 — the textbook value.

21Position grid for sampling and integration

We sample x ∈ [−8, 8] at 4001 evenly-spaced points (Δx = 0.004). The Hermite wavefunctions die exponentially outside the classically allowed region, so this window contains essentially all of the probability for n up to ~10.

22Loop over the lowest six states

We do not need to recompute Hₙ analytically; we just call psi_n(n, x). For each n we evaluate ψₙ on the grid, then verify that ∫|ψₙ|² dx = 1 using the trapezoidal rule.

24Numerical normalization check

`np.trapezoid(psi*psi, x)` approximates ∫|ψ(x)|² dx by summing the trapezoid areas between adjacent samples. For all six states we should get 1.000000 to six digits — the Hermite normalization constant is exact.

25What the print produces

Output: n=0 E=0.5 1.000000; n=1 E=1.5 1.000000; … n=5 E=5.5 1.000000. Every wavefunction has unit total probability and energy E_n = n + ½ — the QHO ladder confirmed numerically.

28Plot the parabolic well

The bare parabola V(x) = x²/2 is the stage on which everything else lives. Energy on the y-axis, position on the x-axis. This is the picture you build all intuition from.

31Stack wavefunctions on their energy levels

Textbook trick: instead of plotting ψₙ in isolation, we plot 0.7·ψₙ(x) + Eₙ so the wavefunction wiggles around the horizontal line corresponding to its own energy. The dashed line at Eₙ is added separately. This is exactly the layout used in the interactive viewer above.

34Why the constant 0.7 is just visual

We rescale ψ by 0.7 only so consecutive curves do not overlap on the figure. The physics — the shape, the number of nodes, the sign — is unchanged. Try setting it to 0.1 and the curves shrink to thin lines on top of the energy levels.

30 lines without explanation

1import math
2import numpy as np
3from numpy.polynomial.hermite import hermval
4import matplotlib.pyplot as plt
5
6# ---------- 1. Energy ladder ------------------------------------------------
7# In natural units (hbar = m = omega = 1), the QHO has
8#     E_n = n + 1/2,   n = 0, 1, 2, ...
9def energy(n):
10    return n + 0.5
11
12# ---------- 2. Eigenfunctions ----------------------------------------------
13# psi_n(x) = (1/sqrt(2^n n!)) * pi^{-1/4} * H_n(x) * exp(-x^2/2)
14def psi_n(n, x):
15    coeff = np.zeros(n + 1)
16    coeff[n] = 1.0                       # selects Hermite polynomial H_n
17    H = hermval(x, coeff)                # evaluate H_n(x) at every x
18    norm = 1.0 / np.sqrt(2.0**n * math.factorial(n)) * np.pi**-0.25
19    return norm * H * np.exp(-0.5 * x * x)
20
21# ---------- 3. Normalization check -----------------------------------------
22x = np.linspace(-8, 8, 4001)
23print("n  E_n      integral |psi|^2 (should be 1)")
24for n in range(6):
25    psi = psi_n(n, x)
26    overlap = np.trapezoid(psi * psi, x)
27    print(f"{n}  {energy(n):.1f}    {overlap:.6f}")
28
29# ---------- 4. Plot the first four states on top of the parabola -----------
30fig, ax = plt.subplots(figsize=(8, 5))
31ax.plot(x, 0.5 * x**2, "k", lw=1.5, label="V(x) = x²/2")
32
33for n in range(4):
34    En = energy(n)
35    psi = psi_n(n, x)
36    # Shift each wavefunction up to its energy level for the textbook layout
37    ax.plot(x, 0.7 * psi + En, lw=2, label=f"ψ_{n}+E_{n}")
38    ax.axhline(En, color="gray", lw=0.5, linestyle="--")
39
40ax.set_xlim(-5, 5); ax.set_ylim(-0.2, 5)
41ax.set_xlabel("x"); ax.set_ylabel("Energy / Wavefunction")
42ax.legend(loc="upper right", fontsize=8)
43plt.tight_layout(); plt.show()

Running this prints six normalisation integrals each equal to 1.000000, and produces the textbook layered plot of $\psi_0$ through $\psi_3$ drawn on the parabola. The Python calculation matches the analytic formulas exactly — no discretisation, no approximation.

PyTorch: Diagonalise the Hamiltonian Directly

The Hermite approach is elegant, but it works only because we already know the answer. Let's flip the problem around: discretise the Hamiltonian on a grid, treat it as a matrix, and ask PyTorch's eigensolver for its eigenvalues and eigenvectors. If we have not made an arithmetic mistake, the lowest eigenvalues had better be $(n+\tfrac{1}{2})\hbar\omega$ .

QHO as a Matrix Eigenvalue Problem

🐍python

Explanation(11)

Code(41)

1PyTorch for a quantum eigenvalue problem

PyTorch is not just for neural networks — its dense tensor math is fast and differentiable, which makes it perfect for small physics problems where we may later want to optimise things like the potential. `torch.linalg.eigh` is a symmetric-matrix eigensolver.

6Discretise space

Quantum eigenstates live on the continuous real line, but our computer only knows arrays. We sample x ∈ [−8, 8] at N = 600 points so Δx ≈ 0.0267. The two derivatives in Ĥ are approximated by finite differences — a method as old as Newton.

10Diagonal potential matrix

In the position basis, V̂ is multiplication by V(x). So as a matrix in our N-dimensional space it is diagonal: Vᵢᵢ = ½xᵢ². The line `0.5 * x * x` is a vector of length N storing that diagonal.

14Three-point kinetic-energy stencil

Taylor expansion gives the centred second derivative: ψ″(x) ≈ (ψ(x+Δx) − 2ψ(x) + ψ(x−Δx))/Δx². Multiplying by −½ (since T̂ = −½ ∂²/∂x²) gives the matrix entries (Tᵢᵢ = 1/Δx², Tᵢ,ᵢ±₁ = −½/Δx²) that we are assembling on the next line.

17Assemble the tridiagonal kinetic operator

`torch.diag(main)` places `1/Δx²` on the main diagonal, while the two `torch.diag(off, ±1)` calls place `−½/Δx²` on the super- and sub-diagonals. The result is the N × N matrix representation of T̂ in the position basis. For N = 600 this is a (600 × 600) tridiagonal matrix.

20Build the full Hamiltonian

Adding `torch.diag(V)` adds the diagonal potential entries to T̂, producing Ĥ = T̂ + V̂. This is just `H[i,i] = 1/Δx² + ½xᵢ²`, `H[i,i±1] = −½/Δx²`, zero elsewhere. The whole matrix has roughly 3N non-zeros.

21Diagonalise

`torch.linalg.eigh` returns eigenvalues in ascending order (`eigvals`) and the corresponding eigenvectors as columns of `eigvecs`. Because Ĥ is Hermitian, eigh is the right call — guaranteed real eigenvalues and orthogonal eigenvectors.

22Compare to the analytic ladder

The printed line gives e.g. `[0.4999, 1.4996, 2.4990, 3.4980, …]`. These are within 0.001 of (n+½) — the only difference is the finite-difference discretisation error. Increasing N and L drives the gap to zero, confirming the exact QHO ladder Eₙ = ℏω(n+½).

25Restore the proper normalisation

`eigh` returns unit-norm vectors in the discrete inner product Σψᵢ². To get ∫|ψ|²dx = 1 we must multiply by 1/√Δx, because Σψᵢ² · Δx ≈ ∫|ψ|² dx. After this rescaling, psi[:, n] is a true probability amplitude.

28Plot eigenstates on their energy levels

Same textbook layout as before — the parabola is drawn first, then each wavefunction shifted up to its own energy. Compare the curves to the analytic Hermite functions: indistinguishable to the eye.

32What the comparison teaches

We never assumed Hermite polynomials. Yet the eigenvalues are (n+½) and the eigenvectors look exactly like Hₙ(x)·exp(−x²/2). This is reassuring: the analytic solution and the brute-force diagonalisation describe the same physical system, two languages for one truth.

30 lines without explanation

1import torch
2import matplotlib.pyplot as plt
3
4# Build the Hamiltonian H = T + V on a finite grid and diagonalize it.
5# We use natural units (hbar = m = omega = 1).
6
7N = 600                         # number of grid points
8L = 8.0                         # half-width of the box [-L, L]
9x = torch.linspace(-L, L, N)
10dx = (x[1] - x[0]).item()
11
12# --- Potential: V(x) = 1/2 * x^2 ----------------------------------------
13V = 0.5 * x * x                 # shape (N,)
14
15# --- Kinetic operator using second-order finite differences -------------
16# T psi(x) ~ -1/2 * (psi(x+dx) - 2 psi(x) + psi(x-dx)) / dx^2
17main = torch.full((N,), 1.0 / dx**2)                # diagonal entries
18off  = torch.full((N - 1,), -0.5 / dx**2)           # super- and sub-diagonal
19T = torch.diag(main) + torch.diag(off, 1) + torch.diag(off, -1)
20
21# --- Full Hamiltonian and its eigen-decomposition -----------------------
22H = T + torch.diag(V)
23eigvals, eigvecs = torch.linalg.eigh(H)             # ascending order
24print("first 6 numerical energies:", eigvals[:6].tolist())
25print("expected (n+1/2):           ", [n + 0.5 for n in range(6)])
26
27# --- Normalise the eigenvectors so ∫|psi|^2 dx = 1 -----------------------
28psis = eigvecs / torch.sqrt(torch.tensor(dx))       # shape (N, N)
29
30# --- Plot the lowest four states ----------------------------------------
31fig, ax = plt.subplots(figsize=(7, 4.5))
32ax.plot(x, V, "k", lw=1.5, label="V(x)")
33for n in range(4):
34    ax.plot(
35        x,
36        eigvals[n].item() + 0.7 * psis[:, n].numpy(),
37        lw=2, label=f"ψ_{n} (E={eigvals[n].item():.3f})",
38    )
39ax.set_xlim(-5, 5); ax.set_ylim(-0.2, 4)
40ax.set_xlabel("x"); ax.legend(fontsize=8)
41plt.tight_layout(); plt.show()

Why this approach generalises: swap

V(x) = \tfrac{1}{2}x^2

for any other smooth potential and the same eight lines find its bound states. This is the finite-difference workhorse used in computational quantum chemistry, condensed matter, and Section 29.8 on numerical methods.

Why It Appears Everywhere

Once you have the QHO ladder, you have the language for an astonishing range of physics:

System	Excitation quantum	Role of ℏω
Vibrating diatomic molecule	Vibron	Spacing of IR absorption lines (CO: 0.27 eV)
Atoms in a crystal lattice	Phonon	Sets specific-heat law, sound velocity
Single mode of the EM field	Photon	ℏω is the photon energy E = hν
Trapped ion in a Paul trap	Motional quantum	Underlies trapped-ion quantum computers
Quantum field (free boson)	Particle	Each k-mode is a QHO; n_k is occupation number
Variational autoencoder latent	—	Gaussian prior shares math with the ground state

The deep reason: the harmonic oscillator is the only non-trivial system whose Hamiltonian is a sum of squares of two operators that obey a canonical commutation relation. Whenever a physical system has those ingredients — and they appear constantly — the QHO machinery is the right tool.

Test Your Understanding

Question 1 of 6

Why are the energy levels of the quantum harmonic oscillator equally spaced by ℏω?

Summary

Drop a quantum particle into a parabolic well and three things happen that have no classical analogue: its energy is quantised in equal steps of $\hbar\omega$ , its lowest rung sits at a non-zero $\tfrac{1}{2}\hbar\omega$ , and its wavefunction leaks gently past the classical turning points. The ladder operators $\hat{a}, \hat{a}^\dagger$ give the whole spectrum almost without computation, the Hermite polynomials give the wavefunctions explicitly, and a small PyTorch script recovers both numerically from the bare Hamiltonian. Coherent states sit on top of this structure as the maximally classical superposition — the same mathematical object that describes laser light, microwave cavities, and many other near-classical quantum systems. Because every smooth potential looks parabolic near its minimum, this one solvable problem is the seed crystal for understanding small oscillations across the whole of physics.

Next up (Section 29.6): we will keep the same Schrödinger machinery but swap the parabola for a finite barrier. The wavefunction leakage we just glimpsed becomes tunnelling — the engine behind alpha decay, scanning tunnelling microscopes, and the working of every flash drive.