Discretised states and matrix operators

Fix a box of length L and N interior points xⱼ = jΔx with Δx = L/(N+1). Sampling replaces |ψ⟩ by the column (ψ₁ … ψ_N), and ∫φ*ψ dx by the sum Δx Σⱼ φⱼ*ψⱼ = Δx φ†ψ, which is the trapezoidal rule with the wall samples set to zero. That weight is not decoration: ψ carries units of m^(−1/2), so ψ†ψ has units of m⁻¹ and only Δx ψ†ψ is dimensionless. Normalisation becomes Δx Σⱼ |ψⱼ|² = 1, and Born's rule for an interval becomes P(a < x < b) = Δx Σ |ψⱼ|² over the samples inside it. Two independent errors have already entered. The quadrature error in the inner product is O(Δx²) for a smooth ψ, and far smaller than that for one vanishing at both ends. The error in the operators, which is what the rest of this topic is about, is a separate account entirely. For L = 1.00 nm and N = 4, Δx = 0.20 nm and the samples sit at 0.20, 0.40, 0.60 and 0.80 nm: four complex numbers standing in for a function.

Multiplication by x acts pointwise, so on the grid it is X = diag(x₁ … x_N) and ⟨x⟩ = Δx Σⱼ xⱼ|ψⱼ|². Nothing is approximated: X is the position operator restricted to the sample points. Its eigenvectors are the unit columns eⱼ, and once rescaled by 1/√Δx they are this space's stand-in for |x⟩ — which is a genuine gain, because the continuum's position eigenkets are δ-normalised distributions outside L² while eⱼ is an honest vector in C^N with an honest eigenvalue xⱼ. It is also a genuine loss. The grid has quietly declared the resolution of a position measurement to be Δx: structure in |ψ|² finer than 0.20 nm is simply not represented, and X² = diag(xⱼ²) inherits the same blindness, so ⟨x²⟩ is a four-term sum rather than an integral.

p̂ = −iħ d/dx has no pointwise action, so it must be built from neighbours. The centred difference gives (Pψ)ⱼ = −iħ(ψⱼ₊₁ − ψⱼ₋₁)/2Δx, so Pⱼₖ = −iħ/2Δx for k = j + 1 and +iħ/2Δx for k = j − 1. Check the adjoint: (P†)ₖⱼ = (Pⱼₖ)*, which for k = j + 1 is (−iħ/2Δx)* = +iħ/2Δx, exactly the entry Pₖⱼ. Hermitian — and it is the centring that does it. Build P from the forward difference (ψⱼ₊₁ − ψⱼ)/Δx instead and the matrix carries entries above the diagonal and nothing below, so ‖P − P†‖ is of order ħ/Δx and grows as the grid is refined. On a real ψ that matrix returns a pure-imaginary ⟨p⟩, which is the fastest way to find the bug. The FFT route differs in kind: P = F†diag(ħkₙ)F with kₙ = 2πn/(NΔx) is exactly Hermitian and spectrally accurate, but it imposes periodic walls, so it solves a ring rather than a box.

With matrices in hand, every question in this unit is linear algebra: ⟨A⟩ = Δx ψ†Aψ, Var(A) = Δx ψ†A²ψ − (Δx ψ†Aψ)², and H|n⟩ = Eₙ|n⟩ is one call to numpy.linalg.eigh. One trap deserves naming. The product P@P is not the three-point kinetic matrix. Applying the centred difference twice reaches out to 2Δx, and its plane-wave symbol is (ħ sin(kΔx)/Δx)², while the three-point Laplacian gives (2ħ sin(kΔx/2)/Δx)². Expand both: the three-point form is low by (kΔx)²/12, P⋅P by (kΔx)²/3, four times worse. At Δx = 0.20 nm and k = 4.0 nm⁻¹, the exact ħ²k²/2m is 0.610 eV, the three-point matrix returns 0.578 eV and P⋅P returns 0.490 eV. Worse still, sin(kΔx) vanishes at kΔx = π, so P⋅P assigns zero kinetic energy to the alternating vector (+1, −1, +1, −1, …) — the odd-even decoupling that makes the three-point second derivative the right choice inside a Hamiltonian.

Cutting the matrix off at N rows is the statement ψ₀ = ψ_N₊₁ = 0: a Dirichlet wall you never wrote down, which raises the energy of any state with amplitude beyond the box. The spacing pushes the other way. The eigenvalues of (ħ²/2mΔx²)tridiag(−1,2,−1) are exactly (2ħ²/mΔx²)sin²(nπΔx/2L), always below ħ²kₙ²/2m, by the factor 1 − (kₙΔx)²/12 to leading order. For an electron in a 1.00 nm box with N = 4, the prefactor ħ²/2mΔx² is 0.952 eV, and the ground level comes out at 0.364 eV against the exact 0.376 eV — 3.2% low, matching the predicted −(πΔx/L)²/12 = −3.3%. The fourth level is 3.45 eV against 6.02 eV, 43% low. That is the general shape: the bottom of a grid spectrum is trustworthy and the top is not, so keep kΔx well under 1 for every level you intend to quote. Beyond the Nyquist cut |k| ≤ π/Δx — 15.7 nm⁻¹, or 9.40 eV of kinetic energy, at Δx = 0.20 nm — the grid cannot represent the state at all.

In C^N, Hermitian and self-adjoint are one condition: a matrix equal to its conjugate transpose, no domain to state, no boundary term to check. That reads as a simplification and is really a concealment. The continuum question — which self-adjoint extension of −(ħ²/2m)d²/dx² you meant, whether −iħ d/dx has one at all on a half-line — has not been answered; the discretisation answered it for you when you chose a stencil and a truncation. The sharpest symptom is the canonical commutator. Since P is nonzero only on the first off-diagonals, [X, P]ⱼₖ = (xⱼ − xₖ)Pⱼₖ gives [X, P] = (iħ/2)(S₊ + S₋), a matrix with no diagonal at all, eigenvalues iħcos(nπ/(N+1)) and trace zero. It cannot equal iħI at any N, because Tr[A, B] = 0 for any two finite matrices while Tr(iħI_N) = iħN. What does converge is the expectation on states the grid resolves: ⟨[X, P]⟩ = iħ[1 − Δx²⟨p²⟩/2ħ²].