Unitary evolution and the propagator

Write the law as d|Ψ⟩/dt = −(i/ħ)Ĥ|Ψ⟩ and it is dz/dt = λz with a vector for z and an operator for λ. The scalar solution is e^(λt)z(0); the operator solution is |Ψ(t)⟩ = exp(−iĤt/ħ)|Ψ(0)⟩, with the exponential defined by the same series, exp(Â) = Σₖ Âᵏ/k!, applied to  = −iĤt/ħ. Two checks make that legitimate. Differentiate the series term by term and you get −(i/ħ)Ĥ multiplying the same series, so Û really does solve the equation, while Û(0) = 𝟙 supplies the initial condition. And there is no ordering ambiguity, because every term is a power of the single operator Ĥ, which commutes with itself — the innocent-looking fact that fails as soon as Ĥ depends on time. The exponent must be dimensionless, and it is: Ĥt/ħ carries J × s ÷ (J s). For an unbounded Ĥ the series does not converge on every vector in the space, which is why the honest definition of Û is the spectral one below. The series is the motivation; the spectral theorem is the machinery.

Any function of a self-adjoint operator is defined by its action on the eigenvectors: if Ĥ|n⟩ = Eₙ|n⟩ then f(Ĥ) = Σₙ f(Eₙ)|n⟩⟨n|, so Û(t) = Σₙ e^{(−iEₙt/ħ)}|n⟩⟨n|. In the energy basis Û is a diagonal matrix of unit-modulus numbers and applying it costs one multiplication per level. That is the whole practical recipe: run scipy.linalg.eigh on the truncated Ĥ once, an O(N³) job, after which any t costs O(N²) — whereas scipy.linalg.expm(−1j*H*t/ℏ) pays O(N³) again for every t you ask for. In two dimensions you can do it in closed form. Because (n̂⋅σ̂)² = 𝟙, the series splits into even and odd parts and gives exp(−iθ n̂⋅σ̂) = cos θ 𝟙 − i sin θ n̂⋅σ̂. So the ammonia Hamiltonian Ĥ = −Δσ̂ₓ has Û(t) = cos(Δt/ħ)𝟙 + i sin(Δt/ħ)σ̂ₓ, and reading off ⟨2|Û|1⟩ = i sin(Δt/ħ) returns P₂ = sin²(Δt/ħ) with no differential equation solved at all. That is what a propagator buys: the dynamics of every initial state, in one object.

Differentiate Û†Û. From dÛ/dt = −(i/ħ)ĤÛ and its dagger dÛ†/dt = +(i/ħ)Û†Ĥ†, the product obeys d(Û†Û)/dt = (i/ħ)Û†(Ĥ† − Ĥ)Û, which vanishes when Ĥ† = Ĥ; since Û†Û = 𝟙 at t = 0, it is 𝟙 forever. Now push on what that buys. Any two states carried by the same Û keep their overlap: ⟨Φ(t)|Ψ(t)⟩ = ⟨Φ₀|Û†Û|Ψ₀⟩ = ⟨Φ₀|Ψ₀⟩. Norm conservation is only the special case Φ = Ψ. Orthogonality survives too, so two perfectly distinguishable states never blur together, and two states you cannot tell apart now you never will. That one line drives the no-cloning theorem: a copier would have to turn the overlap ⟨Φ|Ψ⟩ into ⟨Φ|Ψ⟩², and unitarity forbids any change at all unless the overlap is already 0 or 1. Reversibility comes free from the group law, Û(−t)Û(t) = Û(0) = 𝟙: the propagator has an exact inverse, so nothing is lost between measurements, and any irreversibility in physics is imported from outside this law.

The group law Û(t₁)Û(t₂) = Û(t₁ + t₂) with Û(0) = 𝟙 makes the family (Û(t)) a one-parameter group, and it is strongly continuous because Û(t)|Ψ⟩ → |Ψ⟩ as t → 0. Stone's theorem says every such group is exp(−iĤt/ħ) for exactly one self-adjoint Ĥ, and conversely that every self-adjoint Ĥ generates one. Two consequences matter here. First, expanding for small δt gives Û(δt) = 𝟙 − iĤδt/ħ + O(δt²), so Ĥ = iħ(dÛ/dt) at t = 0: the Hamiltonian is the infinitesimal generator of time translation, and calling its eigenvalues energies is a naming convention laid on top of that. Second, the theorem is strict about self-adjointness and merely symmetric will not do. Take p̂ = −iħ d/dx on the half-line x > 0: symmetric on smooth functions vanishing at the origin, but its deficiency indices are (1,0), so it has no self-adjoint extension. Physically, the group it would generate is a rigid translation, and sliding a function towards the origin pushes probability off the end of the line. Naming a domain is not pedantry — it is the difference between having dynamics and not.

It is tempting to read Û†Û = 𝟙 as a conservation law for energy. It is not. Norm conservation follows from Ĥ† = Ĥ alone. The constancy of ⟨Ĥ⟩ needs a separate and stronger fact: Û is built out of Ĥ, so [Ĥ, Û] = 0, and ⟨Ĥ⟩(t) = ⟨Ψ₀|Û†ĤÛ|Ψ₀⟩ = ⟨Ψ₀|ĤÛ†Û|Ψ₀⟩ = ⟨Ĥ⟩(0). Break that premise and the conclusion goes with it. Drive a two-level atom on resonance with Ĥ(t) = (ħω₀/2)σ̂_z + ħΩ_R cos(ω₀t)σ̂ₓ: the operator is self-adjoint at every instant, so the evolution is unitary and ⟨Ψ|Ψ⟩ never leaves 1 — yet a π-pulse carries the atom from the lower level to the upper one and ⟨Ĥ₀⟩, the mean of the static part Ĥ₀ = (ħω₀/2)σ̂_z, climbs by ħω₀. For an electron spin in 0.500 T that is 5.79 × 10⁻⁵ eV, delivered at 14.0 GHz, and it came from the drive rather than from nowhere. The general statement sits in the same commutator: d⟨Â⟩/dt = (i/ħ)⟨[Ĥ, Â]⟩ + ⟨∂Â/∂t⟩, so any observable commuting with a static Ĥ is conserved, and Ĥ itself is conserved only while it stands still.

For a time-dependent Ĥ the naive answer exp(−(i/ħ)∫₀ᵗ Ĥ(t′)dt′) averages the Hamiltonian and then exponentiates the average — legal only if [Ĥ(t₁), Ĥ(t₂)] = 0 for every pair of times. Baker–Campbell–Hausdorff prices the mistake: e^Âe^B̂ = exp(Â + B̂ + ½[Â, B̂] + …), so what you dropped is the commutator. The correct object is the time-ordered exponential, whose expansion is the Dyson series with later times always to the left. In code, chop the interval into windows short enough that Ĥ barely turns, exponentiate each with Ĥ frozen at its midpoint, and multiply them in time order; the error then falls as the square of the window length. What you must never do is truncate the exponential itself. The Euler step 𝟙 − iĤΔt/ħ multiplies an eigen-component by 1 − iα with α = EΔt/ħ, and |1 − iα| = √(1 + α²) > 1: at α = 0.1, a hundred steps inflate the norm by 64%. The Cayley form (𝟙 − iĤΔt/2ħ)(𝟙 + iĤΔt/2ħ)⁻¹ has modulus exactly 1 by construction, with a phase error of only α − 2 arctan(α/2) ≈ α³/12 per step — which is why Crank–Nicolson, not Euler, is the default in wavepacket codes.