Radial wavefunctions and the Laguerre solution

Start from −(ħ²/2μ)u″ + [ħ²l(l+1)/2μr² − Ze²/4πε₀r]u = Eu on 0 < r < ∞, with u = rR and the condition u(0) = 0 inherited from self-adjointness on the half-line. For a bound state, E < 0, so κ = √(−2μE)/ħ is real and 1/κ is a decay length; set ρ = 2κr. Every constant now collapses into u″ = [1/4 − λ/ρ + l(l+1)/ρ²]u, with the single dimensionless λ = Z/(κa₀) and a₀ = 4πε₀ħ²/μe², which is 0.052918 nm at μ = mₑ. One number carries the energy: fix λ and you have fixed E = −Z²ħ²/(2μa₀²λ²). The equation has two singular points, of different kinds. At ρ = 0 it is regular singular, with Frobenius exponents l+1 and −l. At ρ → ∞ it is irregular, with the two behaviours e^(±ρ/2). Nothing so far says λ must be an integer, and nothing will until infinity is consulted.

A power series about ρ = 0 converges everywhere here, which sounds like good news and is not. A convergent series can happily converge to the growing solution, and it gives no direct grip on which of e^(±ρ/2) you have landed on. So remove both known behaviours by hand: write u(ρ) = ρ^(l+1) e^(−ρ/2) v(ρ), keeping the regular root at the origin and the decaying root at infinity. Substituting and cancelling the common ρ^l e^(−ρ/2) leaves Kummer's equation, ρv″ + (2l + 2 − ρ)v′ + (λ − l − 1)v = 0. That is the whole point of the manoeuvre: the asymptotics are now carried by the prefactors rather than by the series, so any misbehaviour of v shows up in its coefficients instead of hiding inside a slowly converging sum.

Put v = Σⱼ cⱼ ρ^j into Kummer's equation and match powers of ρ: c₍ⱼ₊₁₎ = [(j + l + 1 − λ)/((j+1)(j+2l+2))] cⱼ. Every coefficient descends from c₀, so the regular solution is unique up to normalisation — there is no freedom left to repair anything later. Now look at large j, where the ratio tends to 1/j. That is the large-j ratio of e^ρ = Σⱼ ρ^j/j!, whose own coefficient ratio 1/(j+1) has the same limit, so a generic λ gives v ~ Ce^ρ and therefore u ~ Cρ^(l+1)e^(+ρ/2): smooth, finite at every finite ρ, obedient to u(0) = 0, and hopelessly non-normalisable. The single escape is that the numerator vanish at some j = nᵣ ≥ 0, which kills the next coefficient and every one after it, leaving v a polynomial of degree exactly nᵣ. That demands λ = nᵣ + l + 1, an integer no smaller than l + 1. Call it n.

Setting λ = n fixes κ = Z/(na₀) and hence Eₙ = −(Z²/n²)(ħ²/2μa₀²) = −13.606 Z²/n² eV with μ ≈ mₑ; it also fixes the scaled variable after the fact to ρ = 2Zr/na₀, the form the polynomials are always quoted in. Three readings follow immediately. First, n is not primitive: it is a radial node count plus an angular one, n = nᵣ + l + 1, so l ≤ n − 1 arrives as a consequence rather than a rule to memorise. Second, the surviving polynomial of degree nᵣ is the associated Laguerre polynomial L^(2l+1)₍n−l−1₎(ρ), whose nᵣ positive roots are exactly the interior nodes of u. Third, E depends on n alone, so the sum over l = 0…n−1 of (2l+1) = n² states share one energy — a degeneracy this derivation exhibits but does not explain, which is what the Runge–Lenz vector is for.

None of the above generalises: it is 1/r that makes the recursion two-term. Screen the charge and the series solution dies, so the working method for any other V(r) is diagonalisation. Put rⱼ = jh for j = 1…N with a wall at R = (N+1)h, approximate u″ by (u₍ⱼ₊₁₎ − 2uⱼ + u₍ⱼ₋₁₎)/h², and read off a symmetric tridiagonal matrix: Hⱼⱼ = ħ²/μh² + V_eff(rⱼ), H₍j, j±1₎ = −ħ²/2μh². The domain condition is free — dropping the j = 0 row is precisely u(0) = 0. In atomic units at h = 0.01a₀ with R = 40a₀, the l = 0 eigenvalues come out −13.60535, −3.40140 and −1.51170 eV against the exact −13.60569, −3.40142 and −1.51174 eV. Note the conditioning: the diagonal carries 1/h² = 2.72 × 10⁵ eV, four orders of magnitude above even the ground state, so every level you want is a cancellation between large numbers — reach for a tridiagonal solver, which costs O(N) per eigenvalue and never forms the dense matrix.

The three-point stencil is second order, so every eigenvalue carries an O(h²) bias. Halving h from 0.02a₀ to 0.01a₀ shrinks the ground-state error from 1.36 × 10⁻³ eV to 3.40 × 10⁻⁴ eV, a clean factor of four, and Richardson extrapolation, (4E₍h/2₎ − Eₕ)/3, returns −13.60569 eV. The wall is the other error, and a finer grid does nothing for it. At R = 40a₀ the fourth s level lands at −0.83136 eV against the exact −0.85036 eV: 2.2% too high, and high rather than low, which is what squeezing does. Its mean radius is ⟨r⟩ = (a₀/2Z)(3n² − l(l+1)) = 24a₀, so that wall stood at only 1.7⟨r⟩. Move it to 80a₀ on the same grid and the level becomes −0.850354 eV. Keep R ≳ 3⟨r⟩ for the highest state you want, confirm the node count is n − l − 1, and throw away every eigenvalue near or above zero — those are box states, a chopped-up continuum, not bound levels.