Superposition and eigenbasis expansion

Nothing in the postulates says a particle is "in two places". What they say is that states live in a complex vector space and the Schrödinger equation is linear: if |ψ₁⟩ and |ψ₂⟩ both solve it, so does α|ψ₁⟩ + β|ψ₂⟩ for any complex α and β. A superposition is therefore one vector, not an alternation between two. Two constraints follow. The combination must be renormalised — for orthonormal |1⟩ and |2⟩ that means |α|² + |β|² = 1, while for non-orthogonal kets the cross term 2Re(α*β⟨1|2⟩) joins the sum. And multiplying the whole vector by e^(iθ) changes nothing measurable, so a state is a ray; only the relative phase between components survives. Keep those two phases apart and half the confusion in this topic disappears.

Pick a Hermitian  with a discrete, non-degenerate spectrum. The spectral theorem gives an orthonormal eigenbasis (|n⟩) with ⟨m|n⟩ = δₘₙ and Σₙ|n⟩⟨n| = 1̂. The derivation is then one line: |ψ⟩ = 1̂|ψ⟩ = Σₙ|n⟩⟨n|ψ⟩ = Σₙ cₙ|n⟩ with cₙ = ⟨n|ψ⟩. In the position representation that inner product is the integral cₙ = ∫ψₙ*(x)ψ(x)dx, which is a Fourier coefficient whenever the basis is sinusoidal. The step doing the work is orthonormality: acting with ⟨m| on the sum kills every term but one, so each coefficient is a single projection. Lose orthogonality and you lose that — on a non-orthogonal set the coefficients obey Gc = v, with G the Gram matrix ⟨m|n⟩ and vₘ = ⟨m|ψ⟩, and you must invert G. Orthonormalise first; it is cheaper.

Born's rule reads off the eigenbasis directly: measuring  returns aₙ with probability |cₙ|², and normalisation makes Σₙ|cₙ|² = ⟨ψ|ψ⟩ = 1, which is Parseval's identity rather than a new postulate. Because the modulus discards the argument of cₙ, no measurement of  can see any phase, and anything diagonal in this basis inherits that blindness: ⟨Â⟩ = Σₙ|cₙ|²aₙ. An operator with off-diagonal elements does not. Take |ψ⟩ = (√3/2)|1⟩ + (1/2)e^(iφ)|2⟩ in the infinite well, where ⟨n|x̂|n⟩ = L/2 and ⟨1|x̂|2⟩ = −16L/(9π²) = −0.1801L. Then ⟨x⟩ = L/2 + 2|c₁||c₂|cos φ ⟨1|x̂|2⟩ = 0.500L − 0.156L cos φ, running from 0.344L at φ = 0 to 0.656L at φ = π while ⟨H⟩ sits at 1.75E₁ throughout. Same probabilities, different states.

The resolution of the identity is an assumption about the basis, not a theorem the algebra supplies. For any orthonormal set, complete or not, Bessel's inequality gives Σₙ|cₙ|² ≤ ⟨ψ|ψ⟩; completeness is exactly the case of equality. A short basis therefore loses probability silently — the sum simply comes out under one and nothing in the arithmetic complains. Put numbers on it. The nodeless parabolic state ψ(x) = √(30/L⁵)⋅x(L−x) expands on the box eigenfunctions with cₙ = 4√60/(n³π³) for odd n and zero for even n by parity, giving |c₁|² = 0.99856, |c₃|² = 1.370×10⁻³ and |c₅|² = 6.39×10⁻⁵. One term already carries 99.856% of the probability, but stopping there leaves 1.44×10⁻³ unassigned; keeping n ≤ 3 cuts the deficit to 7.5×10⁻⁵. Auditing Σ|cₙ|² against 1 is the cheapest check you own.

When an eigenvalue a is g-fold degenerate, the label n no longer names a unique vector: any orthonormal basis (|a, k⟩) of the g-dimensional eigenspace will serve. The individual coefficients cₐₖ = ⟨a, k|ψ⟩ are then convention, and only the projector P̂ₐ = Σₖ|a, k⟩⟨a, k| is basis-independent. So the physical quantity is P(a) = ⟨ψ|P̂ₐ|ψ⟩ = Σₖ|cₐₖ|², a sum over the whole block, invariant under any unitary mixing inside it. Concretely, a square two-dimensional box has E = (nₓ² + n_y²)E₁, with E₁ = h²/(8mL²) the one-dimensional ground energy, so |1,2⟩ and |2,1⟩ share E = 5E₁, and the pair (|1,2⟩ ± |2,1⟩)/√2 spans the same eigenspace just as well. A state with c₁₂ = c₂₁ = 1/2 gives P(5E₁) = 1/4 + 1/4 = 1/2 in either basis. Lüders projection then returns P̂ₐ|ψ⟩ ⁄ ‖P̂ₐ|ψ⟩‖ — a vector in the subspace, not one basis element of it.

Discretise on a grid: the state becomes a length-N complex vector, the observable a Hermitian N×N matrix, and numpy.linalg.eigh returns eigenvalues w together with a matrix V whose columns are the |n⟩. The whole expansion is then one matrix product, c = V.conj().T @ ψ, with probabilities p = np.square(np.abs(c)) and the audit p.sum() against 1. Three warnings. eigh fixes each eigenvector only up to a phase, and that phase is library- and machine-dependent, so an individual cₙ's argument means nothing until you pin a convention — force the largest component real and positive, for instance. For a degenerate eigenvalue eigh hands back an arbitrary orthonormal pair spanning the block, so only the summed probability over that block is reproducible. And the grid basis is complete on the grid, so Σ|cₙ|² = 1 there even when the physical basis has been badly truncated: the grid hides the very deficit Bessel's inequality was warning you about.