Radial Hamiltonian and effective potential

Write ψ = R(r)Yₗₘ(θ, φ). The radial part of the Laplacian is R″ + (2/r)R′, and the norm drags the Jacobian along: ∫|R|²r²dr. One substitution disposes of both. Put R = u/r, so R′ = u′/r − u/r² and R″ = u″/r − 2u′/r² + 2u/r³; adding (2/r)R′ = 2u′/r² − 2u/r³ cancels those two correction terms in pairs and leaves R″ + (2/r)R′ = u″/r exactly. The norm follows suit, ∫|R|²r²dr = ∫|u|²dr, so u lives in L²(0,∞) with the flat measure and obeys −(ℏ²/2μ)u″ + V_eff u = Eu. The word to keep is lives. The line has been cut in half: r ∈ (0,∞) has an endpoint, the substitution has made that endpoint singular, and everything difficult about this topic is concentrated there.

L² never acts as a differential operator on the radial factor: on Yₗₘ it returns the number ℏ²l(l+1), and that number is what enters the radial equation as ℏ²l(l+1)/(2μr²). Nothing is rotating; the barrier is bookkeeping for angular momentum the state already carries. It is also universal — the same term appears for a square well, a harmonic trap or a Yukawa potential, so it is not a Coulomb effect. Written in eV with r in Bohr radii, hydrogen's radial problem reads V_eff = (−27.2114/r + 13.6057 l(l+1)/r²) eV. Put l = 1 and r = a₀: the centrifugal coefficient 2 × 13.6057 = 27.2114 exactly matches the Coulomb one, so V_eff(a₀) = 0 identically, not approximately. The floor sits at r = l(l+1)a₀ = 2a₀ with depth −13.6057/[l(l+1)] = −6.803 eV, and the 2p level at −3.401 eV floats halfway between that floor and threshold.

Set Eₙ = V_eff(r) and solve. In Bohr units the quadratic Eₙ r² + 27.2114r − 13.6057 l(l+1) = 0 has roots r± = n²[1 ± √(1 − l(l+1)/n²)] a₀. For 2p that is 1.17 a₀ and 6.83 a₀, a band 5.66 a₀ wide holding ⟨r⟩ = 5a₀; for 3d it runs from 3.80 a₀ to 14.20 a₀. The square root is real only if l(l+1) ≤ n², and for integers that is exactly l ≤ n − 1. Read it off the well: raise l at fixed n and the floor −13.6057/[l(l+1)] eV rises until the level has nowhere to sit. A 2d state would have to live at −3.401 eV in a well whose floor is −2.268 eV. Truncation of the Laguerre series proves the rule; the effective potential is what makes it obvious.

Near r = 0 the centrifugal term beats everything else in the equation: it is O(r⁻²) while the Coulomb term is O(r⁻¹) and E is O(1). The leading balance is therefore u″ = l(l+1)u/r². Substituting u = r^s gives s(s−1)r^(s−2) = l(l+1)r^(s−2), so the indicial equation is s(s−1) = l(l+1), with roots s = l+1 and s = −l. The regular root gives R = u/r ~ r^l, the familiar behaviour of a p or d orbital near the nucleus; the singular root gives R ~ r^(−l−1). Carrying the series one term further fixes the first correction, u = r^(l+1)[1 − Zr/((l+1)a₀) + …], and because it comes from the recursion — and because E first enters two orders later — it is the same for every state of that l: for l = 2 it is 1 − r/3a₀, whether the state is 3d, 4d or 5d.

For l ≥ 1 the argument is normalisation: u ~ r^(−l) makes ∫|u|²dr ~ ∫r^(−2l)dr diverge at the origin, and nothing more need be said — in Weyl's language the operator is limit point at r = 0 and needs no boundary condition at all. The dividing line is sharp: the limit-point test on the half-line asks whether l(l+1) ≥ 3/4, which l ≥ 1 clears and l = 0 fails. For l = 0 the singular root is u → constant, R ~ 1/r, and ∫|R|²r²dr = ∫dr converges. Normalisation cannot reject it. What rejects it is ∇²(1/r) = −4πδ³(r), which means R ~ 1/r satisfies the eigenvalue equation only with an unadvertised point source at the nucleus. Formally the operator is limit circle at r = 0 and admits a one-parameter family of self-adjoint extensions, the Robin conditions cos α u(0) + sin α u′(0) = 0, each of which kills the boundary term of formula five. Dirichlet is the member a point Coulomb nucleus selects: give the nucleus a radius R_N, solve inside and out, then let R_N → 0, and u(0) = 0 is what survives.

Numerically the domain condition is free if the grid is built correctly. Take rᵢ = ih for i = 1…N in a box of radius R = (N+1)h and never create a node at i = 0: the missing row is Dirichlet. In atomic units the matrix is tridiagonal — diagonal 1/h² − Z/rᵢ + l(l+1)/(2rᵢ²), off-diagonal −1/(2h²) — and one symmetric tridiagonal eigensolver call returns the spectrum. Two errors then compete. The three-point second difference is O(h²): for l = 0 in a 40a₀ box the ground state comes out −0.498756 Ha at h = 0.1a₀ and −0.499688 Ha at h = 0.05a₀, and Richardson extrapolation [4E(h/2) − E(h)]/3 = −0.499998 Ha lands on −13.6057 eV. The box is the other error and it bites the top of the spectrum: shrink the box to 10a₀ and 1s is untouched at −0.499687 Ha while 2s reads −0.1128 Ha instead of −0.125 Ha and 3s is pushed above threshold as a pseudo-state.