Mathematical foundations of quantum mechanics

Heated hydrogen does not glow in all colours. It emits light at sharp wavelengths, which Johann Balmer (1885) fitted by $1/\lambda = R\,(1/4 - 1/n^2)$, $n = 3, 4, 5, \dots$ — a law no classical model of an electron orbiting a nucleus could explain, since such an electron should radiate continuously and spiral inwards. The energies of an atom are discrete, and the light it emits carries the differences.

In 1925 Werner Heisenberg replaced electron orbits by arrays of transition amplitudes; Max Born and Pascual Jordan recognised the arrays as infinite matrices and found the rule $PQ - QP = \frac{h}{2\pi i} I$. In 1926 Erwin Schrödinger found a wave equation whose eigenvalues reproduce the energy levels. The two theories looked nothing alike, until John von Neumann (1932) showed that both are statements about operators on a Hilbert space: matrix mechanics lives on $\ell^2$, wave mechanics on $L^2(\R^3)$, and the two spaces are unitarily isomorphic. Energy levels are points of the spectrum of an operator.

The commutation rule already shows that the mathematics must be infinite-dimensional and delicate. If $X$ and $P$ were $n \times n$ matrices, the trace of $XP - PX$ would be $0$, while the trace of $i\hbar I$ is $i\hbar n$. This lesson develops the framework — states, observables, spectra, unitary evolution — and the two most famous consequences: the uncertainty relation and the quantised harmonic oscillator.

The smallest example. A qubit (an electron's spin, a photon's polarisation) has state space $\C^2$. A state is a unit vector $\psi = \alpha e_0 + \beta e_1$ with $|\alpha|^2 + |\beta|^2 = 1$. Measuring "which basis vector" gives $0$ with probability $|\alpha|^2$ and $1$ with probability $|\beta|^2$. Observables are Hermitian matrices; the spin components are the Pauli matrices

Measuring $\sigma_z$ gives an eigenvalue $\pm 1$, with probability $|\langle v_\pm, \psi \rangle|^2$ for the unit eigenvectors $v_\pm$. The state $e_0$ has a sharp value of $\sigma_z$ ($+1$ with certainty) but gives $\pm 1$ with probability $\tfrac12$ each for $\sigma_x$, whose eigenvectors are $(e_0 \pm e_1)/\sqrt 2$. The matrices do not commute — $\sigma_x \sigma_z = -\sigma_z \sigma_x$ — and no state has sharp values for both. Multiplying $\psi$ by a phase $e^{i\theta}$ changes no probability, so states are really unit vectors up to phase; but relative phases matter: $(e_0 + e_1)/\sqrt 2$ and $(e_0 - e_1)/\sqrt 2$ have identical $\sigma_z$ statistics and opposite values of $\sigma_x$.

A particle on a line. Now the state is a function $\psi \in L^2(\R)$ with $\int |\psi|^2 = 1$, and $|\psi(x)|^2$ is the probability density of the position. Position acts by multiplication, $(X\psi)(x) = x \psi(x)$, and momentum by differentiation, $P\psi = -i\hbar \psi'$; integrating by parts, $\langle \phi, XP\psi \rangle - \langle \phi, PX\psi \rangle = i\hbar \langle \phi, \psi \rangle$ for nice functions: $[X, P] = i\hbar$. Neither operator is defined on all of $L^2$ — $x\psi(x)$ need not be square-integrable, and $\psi$ need not be differentiable — and that is not a technicality to be ignored, but the source of real phenomena.

For bounded operators defined everywhere, symmetric and self-adjoint mean the same thing (and are often called Hermitian). For unbounded operators they differ, and only self-adjoint operators have real spectrum and a spectral theorem. A basic criterion: a symmetric operator $A$ is self-adjoint if and only if $A \pm i$ are both surjective (Reed–Simon I, Thm VIII.3).

Grouping equal eigenvalues, $A = \sum_\lambda \lambda E_\lambda$ with $E_\lambda$ the orthogonal projection onto the $\lambda$-eigenspace: the projection-valued measure of the postulates, concentrated on finitely many points.

We do not prove this here; see Hall, Ch. 10, or Reed–Simon I, §VIII.3. The theorem says that every observable is a multiplication operator in suitable coordinates: position already is one, and momentum becomes one after the Fourier transform. Both have spectrum $\R$ and no eigenvectors in $L^2$ — a "plane wave" $e^{ikx}$ is not a state — which is why the postulates are phrased with spectral measures rather than eigenvectors. The functional calculus $f(A) = \int f \, dE_A$ defines $e^{-itH/\hbar}$, $\sqrt{A}$ and indicator functions $E_A(B) = \mathbf 1_B(A)$.

Why not allow unbounded operators everywhere on $\mathcal H$? Because there are none that are symmetric.

So the canonical commutation relation $[X, P] = i\hbar$ forces at least one of $X, P$ to be unbounded, and then by Hellinger–Toeplitz it cannot be defined on the whole space. Domains are unavoidable in quantum mechanics, and statements such as $[X, P] = i\hbar$ hold only on suitable dense subspaces. The Stone–von Neumann theorem (below) explains which representations of the commutation relation are the right ones.

(Proof: Hall, Ch. 10; Reed–Simon I, §VIII.4.) Stone's theorem is the reason observables must be self-adjoint, not merely symmetric: symmetries of a quantum system that depend continuously on a parameter — time translation, space translation, rotation — are strongly continuous unitary groups, and their generators — energy, momentum, angular momentum — are self-adjoint operators. A symmetric operator with no self-adjoint extension generates no dynamics at all, and one with many self-adjoint extensions (like $P_0$ above) generates many different dynamics.

Consequences: the norm $\|\psi(t)\| = 1$ is conserved, and so is the energy distribution, since $U(t)$ commutes with $E_H$. If $H\psi = E\psi$ then $\psi(t) = e^{-iEt/\hbar}\psi$, the same ray: energy eigenstates are stationary, and all their measurement statistics are constant in time. Superpositions of two eigenstates oscillate with angular frequency $(E_m - E_n)/\hbar$ — Bohr's frequency rule, which explains why spectral lines measure energy differences.

For position and momentum, $[X, P] = i\hbar$ on Schwartz functions gives the Heisenberg–Kennard inequality $\sigma_X \sigma_P \ge \hbar/2$ (Kennard, 1927); by approximation it holds for every unit $\psi \in D(X) \cap D(P)$. Equality requires equality in Cauchy–Schwarz, $(P - p_0)\psi = i\lambda (X - x_0)\psi$ with $\lambda > 0$ real, a first-order differential equation whose $L^2$ solutions are exactly the Gaussian wave packets $\psi(x) = C \exp\big(-\frac{\lambda (x - x_0)^2}{2\hbar} + \frac{i p_0 x}{\hbar}\big)$.

Domains matter here too. On $L^2(0, 2\pi)$ take the angular momentum $L = -i\frac{d}{d\theta}$ on periodic functions and the angle $\Phi$ (multiplication by $\theta$). Formally $[\Phi, L] = i$, yet the eigenstate $\psi = e^{in\theta}/\sqrt{2\pi}$ has $\sigma_L = 0$, so the naive inequality $\sigma_\Phi \sigma_L \ge \frac12$ fails. The resolution: $\Phi\psi$ is not periodic, so $\psi \notin D(L\Phi)$ and Robertson's hypothesis is violated.

What the inequality says is that no state has sharp values of two non-commuting observables: it constrains preparations, not the clumsiness of measuring devices. Relations between measurement error and disturbance are different, later theorems (Ozawa, 2003; Busch, Lahti and Werner, 2013).

The Hamiltonian of a particle in a quadratic potential is $H = \frac{P^2}{2m} + \frac12 m\omega^2 X^2$. In units with $\hbar = m = \omega = 1$ it is $H = \frac12(P^2 + X^2)$ on $L^2(\R)$. Work on the Schwartz space $\mathcal S(\R)$, which is dense and mapped into itself by $X$, $P$ and the ladder operators

Here $\langle a^\dagger\phi, \psi \rangle = \langle \phi, a\psi \rangle$ on $\mathcal S$, and $[X, P] = i$ gives $[a, a^\dagger] = \frac12\big(-i[X, P] + i[P, X]\big) = 1$ and $H = a^\dagger a + \frac12 = N + \frac12$, with the number operator $N = a^\dagger a$.

Since the eigenvectors form an orthonormal basis, $H$ extends uniquely to a self-adjoint operator (essential self-adjointness on $\mathcal S$), namely $H = \sum_n (n + \frac12) \langle h_n, \cdot \rangle h_n$, with spectrum exactly $\{\frac12, \frac32, \frac52, \dots\}$. Explicitly $h_n(x) = (2^n n! \sqrt\pi)^{-1/2} H_n(x) e^{-x^2/2}$ with the Hermite polynomials $H_0 = 1$, $H_1 = 2x$, $H_2 = 4x^2 - 2$, $H_{n+1} = 2xH_n - H_n'$. The lowest energy is not $0$ but $\frac12\hbar\omega$: the zero-point energy, forced by the uncertainty relation, since $\langle H \rangle = \frac12(\sigma_P^2 + \sigma_X^2 + \langle P \rangle^2 + \langle X \rangle^2) \ge \sigma_X\sigma_P \ge \frac12$.

Mathlib works with bounded (everywhere-defined) operators, so the formal statements below are about symmetric linear maps A : E →ₗ[ℂ] E on a complex inner product space E (by Hellinger–Toeplitz these are automatically bounded on a Hilbert space) and about elements of normed algebras. Mathlib's inner product ⟪x, y⟫_ℂ is, like the physicists', conjugate-linear in the first slot. First the basic facts about observables: real expectation values (proved here), real eigenvalues (Mathlib), orthogonality of eigenvectors (proved here, as in the lesson) and the finite-dimensional spectral theorem (Mathlib):

/-- **Expectation values of observables are real.** If `A` is symmetric (`⟪A x, y⟫ = ⟪x, A y⟫`),
then `⟪ψ, A ψ⟫` equals its own complex conjugate. Proof: `conj ⟪ψ, Aψ⟫ = ⟪Aψ, ψ⟫ = ⟪ψ, Aψ⟫`, by
conjugate symmetry of the inner product and then symmetry of `A`. -/
theorem expectation_real {A : E →ₗ[ℂ] E} (hA : A.IsSymmetric) (ψ : E) :
    conj ⟪ψ, A ψ⟫_ℂ = ⟪ψ, A ψ⟫_ℂ := by
  rw [inner_conj_symm, hA]

/-- **Eigenvalues of a symmetric operator are real** (`LinearMap.IsSymmetric.conj_eigenvalue_eq_self`):
measured values of an observable are real numbers. -/
theorem eigenvalue_real {A : E →ₗ[ℂ] E} (hA : A.IsSymmetric) {ev : ℂ}
    (h : Module.End.HasEigenvalue A ev) : conj ev = ev :=
  hA.conj_eigenvalue_eq_self h

/-- **Eigenvectors of an observable for different eigenvalues are orthogonal.** If `A v = λ v`
and `A w = μ w` with real `λ ≠ μ`, then `⟪v, w⟫ = 0`. Proof:
`λ ⟪v, w⟫ = ⟪A v, w⟫ = ⟪v, A w⟫ = μ ⟪v, w⟫`, so `(λ − μ) ⟪v, w⟫ = 0`. Distinct measurement
outcomes are perfectly distinguishable states. -/
theorem eigenvectors_orthogonal {A : E →ₗ[ℂ] E} (hA : A.IsSymmetric) {v w : E} {lv lw : ℝ}
    (hv : A v = (lv : ℂ) • v) (hw : A w = (lw : ℂ) • w) (hne : lv ≠ lw) : ⟪v, w⟫_ℂ = 0 := by
  have h := hA v w
  rw [hv, hw, inner_smul_left, inner_smul_right, Complex.conj_ofReal] at h
  have h' : ((lv : ℂ) - lw) * ⟪v, w⟫_ℂ = 0 := by rw [sub_mul, h, sub_self]
  rcases mul_eq_zero.mp h' with h0 | h0
  · exact absurd (by exact_mod_cast sub_eq_zero.mp h0) hne
  · exact h0

/-- **The spectral theorem in finite dimensions.** A symmetric operator on an `n`-dimensional
complex inner product space has an orthonormal basis `e₀, …, e_{n−1}` of eigenvectors with real
eigenvalues `λᵢ`: `A eᵢ = λᵢ eᵢ`. In a state `ψ`, a measurement of `A` gives `λᵢ` with probability
`|⟪eᵢ, ψ⟫|²` (Born's rule). This is Mathlib's `LinearMap.IsSymmetric.apply_eigenvectorBasis`. -/
theorem spectral_theorem_finite [FiniteDimensional ℂ E] {A : E →ₗ[ℂ] E} (hA : A.IsSymmetric)
    {n : ℕ} (hn : Module.finrank ℂ E = n) (i : Fin n) :
    A (hA.eigenvectorBasis hn i) = (hA.eigenvalues hn i : ℂ) • hA.eigenvectorBasis hn i :=
  hA.apply_eigenvectorBasis hn i

The uncertainty relation, proved from scratch exactly as in the human proof: symmetry turns $\langle \psi, [A, B]\psi \rangle$ into $\langle u, w \rangle - \langle w, u \rangle$, then the triangle and Cauchy–Schwarz inequalities finish. The second version allows arbitrary real shifts $a, b$; taking them to be the expectation values gives $\sigma_A \sigma_B \ge \frac12 |\langle [A, B] \rangle|$.

/-- **Robertson's inequality, core estimate.** For symmetric `A, B` and any vector `ψ`,
`|⟪ψ, [A, B] ψ⟫| ≤ 2 ‖Aψ‖ ‖Bψ‖`, where `[A, B] = AB − BA`. Proof: by symmetry
`⟪ψ, ABψ⟫ = ⟪Aψ, Bψ⟫` and `⟪ψ, BAψ⟫ = ⟪Bψ, Aψ⟫`; then the triangle inequality and the
Cauchy–Schwarz inequality `|⟪u, v⟫| ≤ ‖u‖ ‖v‖`. -/
theorem robertson {A B : E →ₗ[ℂ] E} (hA : A.IsSymmetric) (hB : B.IsSymmetric) (ψ : E) :
    ‖⟪ψ, (A ∘ₗ B - B ∘ₗ A) ψ⟫_ℂ‖ ≤ 2 * ‖A ψ‖ * ‖B ψ‖ := by
  have h1 : ⟪ψ, (A ∘ₗ B - B ∘ₗ A) ψ⟫_ℂ = ⟪A ψ, B ψ⟫_ℂ - ⟪B ψ, A ψ⟫_ℂ := by
    rw [LinearMap.sub_apply, inner_sub_right, LinearMap.comp_apply, LinearMap.comp_apply,
      ← hA ψ (B ψ), ← hB ψ (A ψ)]
  rw [h1]
  calc ‖⟪A ψ, B ψ⟫_ℂ - ⟪B ψ, A ψ⟫_ℂ‖ ≤ ‖⟪A ψ, B ψ⟫_ℂ‖ + ‖⟪B ψ, A ψ⟫_ℂ‖ := norm_sub_le _ _
    _ ≤ ‖A ψ‖ * ‖B ψ‖ + ‖B ψ‖ * ‖A ψ‖ :=
        add_le_add (norm_inner_le_norm _ _) (norm_inner_le_norm _ _)
    _ = 2 * ‖A ψ‖ * ‖B ψ‖ := by ring

/-- **The Robertson–Heisenberg uncertainty relation.** For symmetric `A, B`, any `ψ` and any real
numbers `a, b`: `½ |⟪ψ, [A, B] ψ⟫| ≤ ‖(A − a)ψ‖ · ‖(B − b)ψ‖`. For a unit vector `ψ` and
`a = ⟨A⟩ = ⟪ψ, Aψ⟫`, `b = ⟨B⟩`, the right side is the product of standard deviations `σ_A σ_B`.
Proof: shifting by real scalars does not change the commutator, and the shifted operators are
still symmetric; the cross terms `a ⟪ψ, Bψ⟫` and `b ⟪Aψ, ψ⟫` cancel in pairs, leaving
`⟪Aψ, Bψ⟫ − ⟪Bψ, Aψ⟫`, which is bounded as in `robertson`. -/
theorem uncertainty {A B : E →ₗ[ℂ] E} (hA : A.IsSymmetric) (hB : B.IsSymmetric) (ψ : E)
    (a b : ℝ) :
    ‖⟪ψ, (A ∘ₗ B - B ∘ₗ A) ψ⟫_ℂ‖ / 2 ≤ ‖A ψ - (a : ℂ) • ψ‖ * ‖B ψ - (b : ℂ) • ψ‖ := by
  set u := A ψ - (a : ℂ) • ψ
  set w := B ψ - (b : ℂ) • ψ
  have hAψ : ⟪A ψ, ψ⟫_ℂ = ⟪ψ, A ψ⟫_ℂ := hA ψ ψ
  have hBψ : ⟪B ψ, ψ⟫_ℂ = ⟪ψ, B ψ⟫_ℂ := hB ψ ψ
  have h1 : ⟪ψ, (A ∘ₗ B - B ∘ₗ A) ψ⟫_ℂ = ⟪u, w⟫_ℂ - ⟪w, u⟫_ℂ := by
    simp only [u, w, LinearMap.sub_apply, LinearMap.comp_apply, inner_sub_right, inner_sub_left,
      inner_smul_left, inner_smul_right, Complex.conj_ofReal, ← hA ψ (B ψ), ← hB ψ (A ψ), hAψ,
      hBψ]
    ring
  rw [h1]
  have h2 : ‖⟪u, w⟫_ℂ - ⟪w, u⟫_ℂ‖ ≤ 2 * (‖u‖ * ‖w‖) :=
    calc ‖⟪u, w⟫_ℂ - ⟪w, u⟫_ℂ‖ ≤ ‖⟪u, w⟫_ℂ‖ + ‖⟪w, u⟫_ℂ‖ := norm_sub_le _ _
      _ ≤ ‖u‖ * ‖w‖ + ‖w‖ * ‖u‖ := add_le_add (norm_inner_le_norm _ _) (norm_inner_le_norm _ _)
      _ = 2 * (‖u‖ * ‖w‖) := by ring
  linarith

The ladder-operator algebra holds in any ring (noncomm_ring verifies the commutator identities), the eigenvector-shifting property in any complex vector space, and positivity of the number operator in any inner product space where $a^\dagger$ is the adjoint of $a$:

/-- **The number operator raises with the creation operator.** In any ring, if `a a† − a† a = 1`
and `N = a† a`, then `N a† − a† N = a†`. Proof: `N a† − a† N = a† (a a† − a† a) = a†`. -/
theorem number_raise {R : Type*} [Ring R] (a ad : R) (h : a * ad - ad * a = 1) :
    (ad * a) * ad - ad * (ad * a) = ad :=
  calc (ad * a) * ad - ad * (ad * a) = ad * (a * ad - ad * a) := by noncomm_ring
    _ = ad := by rw [h, mul_one]

/-- **…and lowers with the annihilation operator:** `N a − a N = −a`, because
`N a − a N = −(a a† − a† a) a`. -/
theorem number_lower {R : Type*} [Ring R] (a ad : R) (h : a * ad - ad * a = 1) :
    (ad * a) * a - a * (ad * a) = -a :=
  calc (ad * a) * a - a * (ad * a) = -((a * ad - ad * a) * a) := by noncomm_ring
    _ = -a := by rw [h, one_mul]

/-- **Ladder operators move between eigenvectors.** If `a a† − a† a = 1` on a complex vector space
and `N v = λ v` for `N = a† a`, then `N (a† v) = (λ + 1) a† v` and `N (a v) = (λ − 1) a v`. Each
application of `a†` adds one quantum `ℏω` of energy and each application of `a` removes one. -/
theorem raise_eigenvector {V : Type*} [AddCommGroup V] [Module ℂ V] (a ad : Module.End ℂ V)
    (h : a * ad - ad * a = 1) {v : V} {l : ℂ} (hv : (ad * a) v = l • v) :
    (ad * a) (ad v) = (l + 1) • ad v ∧ (ad * a) (a v) = (l - 1) • a v := by
  have hv' : ad (a v) = l • v := by simpa only [Module.End.mul_apply] using hv
  have hr := congrArg (fun T : Module.End ℂ V => T v) (number_raise a ad h)
  have hl := congrArg (fun T : Module.End ℂ V => T v) (number_lower a ad h)
  simp only [LinearMap.sub_apply, Module.End.mul_apply, LinearMap.neg_apply] at hr hl
  rw [hv', map_smul] at hr hl
  simp only [Module.End.mul_apply]
  rw [sub_eq_iff_eq_add] at hr hl
  constructor
  · rw [hr, add_smul, one_smul]
    abel
  · rw [hl, sub_smul, one_smul]
    abel

/-- **The number operator is positive, so the ladder must stop.** If `a†` is the adjoint of `a`
(`⟪a† x, y⟫ = ⟪x, a y⟫`), then `⟪v, a† a v⟫ = ‖a v‖²`. Hence an eigenvalue `λ` of `N = a† a` with
eigenvector `v ≠ 0` satisfies `λ ‖v‖² = ‖a v‖² ≥ 0`: repeated lowering must end at a vector with
`a v₀ = 0`, the ground state, and the eigenvalues of `N` are `0, 1, 2, …`. -/
theorem number_nonneg (a ad : E →ₗ[ℂ] E) (hadj : ∀ x y, ⟪ad x, y⟫_ℂ = ⟪x, a y⟫_ℂ) (v : E) :
    ⟪v, ad (a v)⟫_ℂ = ((‖a v‖ ^ 2 : ℝ) : ℂ) := by
  rw [← inner_conj_symm, hadj, inner_self_eq_norm_sq_to_K]
  simp

The two no-go theorems: the trace argument for matrices, Wielandt's proof for any nontrivial normed algebra (induction for the commutator with $p^{n+1}$, induction for $p^n \ne 0$, the norm estimate, and the Archimedean property), and its corollary for bounded operators on a nonzero normed space:

/-- **No matrices satisfy `[X, P] = c · 1` with `c ≠ 0`.** Taking traces: `tr(XP) = tr(PX)`, so the
left side has trace `0`, while `tr(c · 1) = c n ≠ 0` for `n ≥ 1`. Hence position and momentum
(`[X, P] = iℏ`) cannot be represented on a finite-dimensional Hilbert space. -/
theorem no_finite_ccr {n : ℕ} (hn : 0 < n) (X P : Matrix (Fin n) (Fin n) ℂ) {c : ℂ}
    (hc : c ≠ 0) : X * P - P * X ≠ c • (1 : Matrix (Fin n) (Fin n) ℂ) := by
  intro h
  have ht := congrArg Matrix.trace h
  rw [Matrix.trace_sub, Matrix.trace_mul_comm, sub_self, Matrix.trace_smul, Matrix.trace_one,
    Fintype.card_fin, smul_eq_mul] at ht
  have hn' : (n : ℂ) ≠ 0 := by exact_mod_cast hn.ne'
  exact mul_ne_zero hc hn' ht.symm

/-- **The Wintner–Wielandt theorem.** In a nontrivial normed algebra (for example the bounded
operators on a nonzero Banach space), no elements satisfy `x p − p x = c · 1` with `c ≠ 0`.
Wielandt's proof: by induction `x pⁿ⁺¹ − pⁿ⁺¹ x = (n + 1) c pⁿ`; hence no power of `p` vanishes
(otherwise the next lower one would too, down to `p⁰ = 1`); taking norms,
`(n + 1) |c| ‖pⁿ‖ ≤ 2 ‖x‖ ‖p‖ ‖pⁿ‖`, so `(n + 1) |c| ≤ 2 ‖x‖ ‖p‖` for every `n`, which is absurd. -/
theorem wintner_wielandt {A : Type*} [NormedRing A] [NormedAlgebra ℂ A] [Nontrivial A]
    (x p : A) {c : ℂ} (hc : c ≠ 0) : x * p - p * x ≠ c • (1 : A) := by
  intro h
  -- The commutator with a power of `p`.
  have key : ∀ n : ℕ, x * p ^ (n + 1) - p ^ (n + 1) * x = (((n : ℂ) + 1) * c) • p ^ n := by
    intro n
    induction n with
    | zero => simp [h]
    | succ n ih =>
      have e : x * p ^ (n + 1 + 1) - p ^ (n + 1 + 1) * x =
          (x * p ^ (n + 1) - p ^ (n + 1) * x) * p + p ^ (n + 1) * (x * p - p * x) := by
        rw [pow_succ p (n + 1)]
        noncomm_ring
      rw [e, ih, h, smul_mul_assoc, mul_smul_comm, mul_one, ← pow_succ, ← add_smul]
      congr 1
      push_cast
      ring
  -- No power of `p` is zero.
  have hne : ∀ n : ℕ, p ^ n ≠ 0 := by
    intro n
    induction n with
    | zero => simp
    | succ n ih =>
      intro hp
      have hk := key n
      rw [hp, mul_zero, zero_mul, sub_zero] at hk
      have hcoef : ((n : ℂ) + 1) * c ≠ 0 := mul_ne_zero (by exact_mod_cast n.succ_ne_zero) hc
      exact ih ((smul_eq_zero.mp hk.symm).resolve_left hcoef)
  -- The norm estimate.
  have hbound : ∀ n : ℕ, ((n : ℝ) + 1) * ‖c‖ ≤ 2 * ‖x‖ * ‖p‖ := by
    intro n
    have hpos : 0 < ‖p ^ n‖ := norm_pos_iff.mpr (hne n)
    have hp1 : ‖p ^ (n + 1)‖ ≤ ‖p ^ n‖ * ‖p‖ := by rw [pow_succ]; exact norm_mul_le _ _
    have hlhs : ‖(((n : ℂ) + 1) * c) • p ^ n‖ = ((n : ℝ) + 1) * ‖c‖ * ‖p ^ n‖ := by
      rw [norm_smul, norm_mul]
      congr 2
      exact_mod_cast Complex.norm_natCast (n + 1)
    have hrhs : ‖x * p ^ (n + 1) - p ^ (n + 1) * x‖ ≤ 2 * ‖x‖ * ‖p‖ * ‖p ^ n‖ :=
      calc ‖x * p ^ (n + 1) - p ^ (n + 1) * x‖
          ≤ ‖x * p ^ (n + 1)‖ + ‖p ^ (n + 1) * x‖ := norm_sub_le _ _
        _ ≤ ‖x‖ * ‖p ^ (n + 1)‖ + ‖p ^ (n + 1)‖ * ‖x‖ :=
            add_le_add (norm_mul_le _ _) (norm_mul_le _ _)
        _ ≤ ‖x‖ * (‖p ^ n‖ * ‖p‖) + (‖p ^ n‖ * ‖p‖) * ‖x‖ := by gcongr
        _ = 2 * ‖x‖ * ‖p‖ * ‖p ^ n‖ := by ring
    rw [key n, hlhs] at hrhs
    exact le_of_mul_le_mul_right hrhs hpos
  -- Contradiction with the Archimedean property.
  have hcpos : 0 < ‖c‖ := norm_pos_iff.mpr hc
  obtain ⟨n, hn⟩ := exists_nat_gt (2 * ‖x‖ * ‖p‖ / ‖c‖)
  rw [div_lt_iff₀ hcpos] at hn
  have := hbound n
  nlinarith

/-- **Position and momentum cannot both be bounded.** On a nonzero complex normed space there
are no bounded operators `X, P` with `XP − PX = iℏ · 1` for `ℏ ≠ 0`. This is Wintner–Wielandt
applied to the normed algebra of bounded operators; so at least one of `X, P` is unbounded, and
then (by the Hellinger–Toeplitz theorem) it cannot be defined on the whole Hilbert space. -/
theorem no_bounded_ccr {F : Type*} [NormedAddCommGroup F] [NormedSpace ℂ F] [Nontrivial F]
    (X P : F →L[ℂ] F) {ℏ : ℝ} (hℏ : ℏ ≠ 0) :
    X * P - P * X ≠ (Complex.I * ℏ) • (1 : F →L[ℂ] F) :=
  wintner_wielandt X P (mul_ne_zero Complex.I_ne_zero (by exact_mod_cast hℏ))

Finally, the bounded case of unitary evolution: in a C*-algebra, $\exp(itH)$ for self-adjoint $H$ is unitary and the evolution operators form a one-parameter group.

/-- **Time evolution is a one-parameter unitary group** (bounded case). For a self-adjoint
element `H` of a C*-algebra (e.g. a Hermitian matrix), `U(t) = exp(i t H)` is unitary
(`selfAdjoint.expUnitary`) and `U(s + t) = U(s) U(t)`, because `sH` and `tH` commute
(`Commute.expUnitary_add`). The physicists' `e^{−itH/ℏ}` is `U(−t/ℏ)`. Stone's theorem — every
strongly continuous unitary group arises this way from a possibly unbounded self-adjoint
generator — is not yet in Mathlib. -/
theorem evolution_group {A : Type*} [NormedRing A] [NormedAlgebra ℂ A] [StarRing A]
    [ContinuousStar A] [CompleteSpace A] [StarModule ℂ A] (H : selfAdjoint A) (s t : ℝ) :
    selfAdjoint.expUnitary ((s + t) • H) =
      selfAdjoint.expUnitary (s • H) * selfAdjoint.expUnitary (t • H) := by
  rw [add_smul]
  exact Commute.expUnitary_add (((Commute.refl (H : A)).smul_left s).smul_right t)

Not formalised (and not yet in Mathlib): unbounded operators with domains, self-adjointness in the unbounded sense, the spectral theorem for unbounded self-adjoint operators, Stone's theorem, position and momentum on $L^2(\R)$, and the completeness of the Hermite functions. The formal development covers the algebra and the bounded-operator theorems on which those analytic results build — and, in the Wintner–Wielandt theorem, proves formally why the unbounded theory is unavoidable.

The first experiment replaces $H = -\frac12 \frac{d^2}{dx^2} + \frac12 x^2$ by a finite-difference matrix on a grid in $[-8, 8]$ and computes its lowest eigenvalues. They approach $n + \frac12$, and halving the grid spacing divides the error by about $4$ (a second-order method). The ground state's uncertainty product comes out as $\frac12$.

"""The harmonic oscillator spectrum from a finite-difference Hamiltonian.

In units hbar = m = omega = 1 the Hamiltonian is H = -(1/2) d^2/dx^2 + x^2/2.
We replace it by a symmetric tridiagonal matrix on a grid of N points in
[-L, L] (Dirichlet boundary) and compute its lowest eigenvalues with numpy.
They approach n + 1/2, and halving the grid spacing divides the error by
about 4 (a second-order method). From the ground-state eigenvector we also
compute the uncertainty product sigma_x sigma_p, which is 1/2: the minimum
allowed by the Heisenberg relation.
"""

import numpy as np

def hamiltonian(n_points, half_width):
    """Finite-difference matrix for -(1/2) d^2/dx^2 + x^2/2 and its grid."""
    if n_points < 3:
        raise ValueError("need at least 3 grid points")
    x = np.linspace(-half_width, half_width, n_points)
    h = x[1] - x[0]
    main = 1.0 / h ** 2 + 0.5 * x ** 2
    off = np.full(n_points - 1, -0.5 / h ** 2)
    return np.diag(main) + np.diag(off, 1) + np.diag(off, -1), x

def lowest_levels(n_points, half_width, count):
    """The `count` smallest eigenvalues of the discretised Hamiltonian."""
    matrix, _ = hamiltonian(n_points, half_width)
    return np.linalg.eigvalsh(matrix)[:count]

def ground_state_uncertainty(n_points, half_width):
    """sigma_x * sigma_p for the lowest eigenvector (grid approximation)."""
    matrix, x = hamiltonian(n_points, half_width)
    h = x[1] - x[0]
    _, vecs = np.linalg.eigh(matrix)
    psi = vecs[:, 0] / np.sqrt(h)
    mean_x = h * np.sum(x * psi ** 2)
    var_x = h * np.sum((x - mean_x) ** 2 * psi ** 2)
    dpsi = np.gradient(psi, h)
    var_p = h * np.sum(dpsi ** 2)
    return float(np.sqrt(var_x * var_p))

def main():
    half_width, count = 8.0, 6
    print("Lowest eigenvalues of the discretised oscillator (exact: n + 1/2):")
    print("     N       h       E_0       E_1       E_2       E_5")
    errors = []
    for n_points in [161, 321, 641]:
        levels = lowest_levels(n_points, half_width, count)
        h = 2 * half_width / (n_points - 1)
        print(f"  {n_points:>4}   {h:.3f}   {levels[0]:.5f}   {levels[1]:.5f}"
              f"   {levels[2]:.5f}   {levels[5]:.5f}")
        errors.append(abs(levels[5] - 5.5))
    print("error in E_5 for the three grids:",
          ", ".join(f"{e:.2e}" for e in errors))
    print("ratios:", ", ".join(f"{a / b:.2f}" for a, b in zip(errors, errors[1:])),
          "(about 4: the error is O(h^2))")
    print()
    u = ground_state_uncertainty(641, half_width)
    print(f"Ground state: sigma_x * sigma_p = {u:.5f}; the exact value is the")
    print("Heisenberg minimum 1/2 (the difference is discretisation error).")
    print("The spacing E_(n+1) - E_n = 1 = hbar omega is the energy quantum,")
    print("and E_0 = 1/2 is the zero-point energy.")
    print()
    print("Epistemic status: floating-point eigenvalues of a finite matrix")
    print("approximating an unbounded operator on a truncated interval. The")
    print("agreement with n + 1/2 is numerical evidence; the exact spectrum is")
    print("the theorem proved in the lesson.")

if __name__ == "__main__":
    main()

The second truncates the ladder operators to $N$ levels. The commutator $[a, a^\dagger]$ is then the identity except in the last diagonal entry, $1 - N$, so that its trace is $0$ — the trace argument made visible: truncation must break the commutation relation somewhere. The number operator still has eigenvalues $0, \dots, N - 1$, and Robertson's inequality holds for all states, with equality for the ground state.

"""Ladder operators on N levels: where the commutation relation must break.

On the span of the first N oscillator states, the annihilation operator a is
the N x N matrix with sqrt(1), ..., sqrt(N - 1) above the diagonal. We compute
the commutator [a, a^dagger] exactly (its entries are integers): it is the
identity except for the last diagonal entry, 1 - N, which makes the trace 0,
as it must be for any commutator of matrices. With X = (a + a^dagger)/sqrt 2
and P = (a - a^dagger)/(i sqrt 2) we check the Robertson inequality
sigma_X sigma_P >= |<[X, P]>| / 2 for number states and random states.
"""

import math
import random

import numpy as np

def ladder_commutator_diagonal(n_levels):
    """Diagonal of [a, a^dagger] computed exactly: entries (k+1) - k, last 1-N."""
    if n_levels < 1:
        raise ValueError("need at least one level")
    a_ad = [k + 1 if k < n_levels - 1 else 0 for k in range(n_levels)]
    ad_a = list(range(n_levels))
    return [p - q for p, q in zip(a_ad, ad_a)]

def annihilation(n_levels):
    """The N x N annihilation matrix a (complex)."""
    return np.diag(np.sqrt(np.arange(1, n_levels)), 1).astype(complex)

def position_momentum(n_levels):
    """X = (a + a^dagger)/sqrt(2), P = (a - a^dagger)/(i sqrt(2))."""
    a = annihilation(n_levels)
    ad = a.conj().T
    return (a + ad) / math.sqrt(2), (a - ad) / (1j * math.sqrt(2))

def robertson_sides(state, x_op, p_op):
    """(sigma_X * sigma_P, |<[X, P]>| / 2) for a normalised state vector."""
    def mean(op):
        return np.vdot(state, op @ state)
    var_x = (mean(x_op @ x_op) - mean(x_op) ** 2).real
    var_p = (mean(p_op @ p_op) - mean(p_op) ** 2).real
    comm = mean(x_op @ p_op - p_op @ x_op)
    return math.sqrt(max(var_x, 0) * max(var_p, 0)), abs(comm) / 2

def main():
    n = 8
    diag = ladder_commutator_diagonal(n)
    print(f"N = {n}: diagonal of [a, a^dagger] = {diag}")
    print(f"trace = {sum(diag)} (a commutator of matrices always has trace 0)")
    a = annihilation(n)
    number = np.linalg.eigvalsh(a.conj().T @ a)
    print("eigenvalues of a^dagger a:",
          " ".join(f"{float(v):.3f}" for v in number))
    x_op, p_op = position_momentum(n)
    print()
    print("Number states |k>: sigma_X sigma_P versus |<[X,P]>|/2")
    for k in range(n):
        state = np.zeros(n, dtype=complex)
        state[k] = 1
        lhs, rhs = robertson_sides(state, x_op, p_op)
        print(f"  k = {k}: {lhs:.4f} >= {rhs:.4f}")
    print("For k < N - 1 these are k + 1/2 >= 1/2, as for the true oscillator;")
    print("the top level feels the truncation: <[X, P]> = i(1 - N) there.")
    rng = random.Random(1925)
    worst = math.inf
    for _ in range(200):
        v = np.array([complex(rng.gauss(0, 1), rng.gauss(0, 1)) for _ in range(n)])
        v /= np.linalg.norm(v)
        lhs, rhs = robertson_sides(v, x_op, p_op)
        worst = min(worst, lhs - rhs)
    print(f"200 random states: min of sigma_X sigma_P - |<[X,P]>|/2 = {worst:.4f}")
    print()
    print("Epistemic status: the commutator diagonal is exact integer")
    print("arithmetic; the eigenvalues and uncertainty products are floating")
    print("point for one N and finitely many states. Robertson's inequality")
    print("and the trace argument are proved in the lesson and in Lean.")

if __name__ == "__main__":
    main()