Orbitals, radial densities and probability clouds

Three densities live inside one state and answer three different questions. |ψ(r)|² = |Rₙₗ|²|Yₗ^m|² is probability per unit volume — what a point detector of volume dV registers. Integrating over angles gives P(r) = r²|Rₙₗ(r)|² = |u(r)|², probability per unit radius, which is what a spherical shell of thickness dr collects. The r² is a Jacobian, not a fudge: the shell at 2a₀ has four times the volume of the shell at a₀. For 1s, |ψ|² = e^(−2r/a₀)/πa₀³ falls monotonically from 2.15 × 10³⁰ m⁻³ at the nucleus, while P(r) = (4/a₀³)r²e^(−2r/a₀) starts at zero, peaks at exactly a₀, and encloses 1 − 5e⁻² = 32.3% of the probability inside a₀ and 1 − 13e⁻⁴ = 76.2% inside 2a₀. Both densities are physical, and the value at the origin is not an artefact: |ψ(0)|² is precisely what the Fermi contact term of hyperfine structure integrates against, and it is non-zero only for l = 0.

The qualitative curve follows from the labels. Write Rₙₗ(r) = r^l · L(r) · e^{(−Zr/na₀)}, with L an associated Laguerre polynomial of degree n − l − 1. The exponential never vanishes and r^l vanishes only at the origin, so every zero at finite r > 0 comes from L: exactly n − l − 1 radial nodes. Concretely, with ρ = r/a₀ and Z = 1: R₃₀ ∝ (27 − 18ρ + 2ρ²)e^(−ρ/3) has roots at ρ = (9 ± 3√3)/2 = 1.90 and 7.10; R₃₁ ∝ ρ(6 − ρ)e^(−ρ/3) has one, at ρ = 6; R₃₂ ∝ ρ²e^(−ρ/3) has none. The angular count needs care about which basis you are in. The polar factor Pₗ^|m|(cos θ) has l − |m| zeros in 0 < θ < π, which are nodal cones; the azimuthal factor e^(imφ) of a complex Yₗ^m has no zeros at all, so |2,1,±1⟩ carries no nodal surface whatever. The familiar |m| planes through the z-axis belong to the REAL combinations, where cos(mφ) and sin(mφ) genuinely vanish; there the angular count is |m| planes plus l − |m| cones = l, and the textbook total n − 1 holds across a shell. The radial count, which is basis-independent, is also the cheapest numerical diagnostic available. Diagonalise the radial Hamiltonian on a grid with u(0) = 0 and count sign changes in each eigenvector: the k-th eigenvector at fixed l, counting from k = 0, must show exactly k of them, and if it does not, the grid is too coarse or the box too small.

A density has no single size. ⟨1/r⟩ = Z/n²a₀ follows in one line from Hellmann–Feynman: ∂H/∂Z = −e²/4πε₀r, so differentiating E = −Z² × 13.606 eV/n² with respect to Z returns ⟨1/r⟩ directly, and no l appears anywhere. ⟨r⟩ = (a₀/2Z)[3n² − l(l+1)] does depend on l: for n = 3 it is 13.5a₀ = 714 pm, 12.5a₀ = 661 pm and 10.5a₀ = 556 pm for l = 0, 1, 2 — a 22% spread across states of identical energy, because the centrifugal barrier ħ²l(l+1)/2μr² keeps high-l clouds out of the small-r region while low-l states put weight there and pay for it further out. Meanwhile 1/⟨1/r⟩ = n²a₀/Z = 9a₀ for all three, equal to ⟨r⟩ for none of them; Jensen's inequality guarantees ⟨1/r⟩ ≥ 1/⟨r⟩ for any spread density. Higher inverse moments are more l-sensitive still: ⟨1/r²⟩ = Z²/[a₀²n³(l + ½)] and ⟨1/r³⟩ = Z³/[a₀³n³l(l + ½)(l + 1)], the latter defined only for l ≥ 1 because the true integral ∫|R|²r⁻¹dr diverges logarithmically when R(0) ≠ 0. Be careful what those moments buy you. ⟨1/r³⟩ is what the spin–orbit term is built from, and that term is l-dependent piece by piece, but in hydrogen it combines with the relativistic kinetic correction so that the total order-α⁴ fine structure depends only on n and j: 2s₁/₂ and 2p₁/₂ stay degenerate at that order, and only the Lamb shift separates them. What actually lifts the l-degeneracy in a many-electron atom is screening, not spin–orbit.

In a stationary state d⟨r⋅p⟩/dt = 0, and evaluating the commutator gives 2⟨T⟩ = ⟨r⋅∇V⟩. The Coulomb potential is homogeneous of degree −1, so r⋅∇V = −V and 2⟨T⟩ = −⟨V⟩. With E = ⟨T⟩ + ⟨V⟩ this fixes both pieces from the eigenvalue alone: ⟨T⟩ = −Eₙ = +13.606 Z²/n² eV and ⟨V⟩ = 2Eₙ = −27.211 Z²/n² eV, with no integral performed. The cross-check is ⟨V⟩ = −(Ze²/4πε₀)⟨1/r⟩ = −27.211 Z²/n² eV, the same number. Two payoffs follow. First ⟨p²⟩ = 2mₑ⟨T⟩ = ħ²Z²/n²a₀², so √⟨p²⟩ = ħZ/na₀, and for 1s the radial spread is Δr = √(⟨r²⟩ − ⟨r⟩²) = (√3/2)a₀, giving a product 0.87ħ — the atom's size is visibly an uncertainty-principle compromise rather than an orbit. (Note this is an estimate built from √⟨p²⟩ and a coordinate spread, not the canonical Δr Δpᵣ relation: radial momentum has no self-adjoint extension on the half-line and is not an observable.) Second vᵣₘₛ/c = √(2⟨T⟩/mₑ)/c = Zα/n, which is 7.30 × 10⁻³ for hydrogen's ground state: that is why the non-relativistic treatment works, and why corrections enter at relative order (Zα/n)² = 5.33 × 10⁻⁵ — a scale of α²|E₁| = 7.2 × 10⁻⁴ eV, of which the actual 1s fine-structure shift is about a quarter, −1.81 × 10⁻⁴ eV. For a 1s electron in gold, Zα = 0.58, and the same estimate says this model has run out.

The textbook pₓ, p_y and d_xy pictures are not extra states. Inside the degenerate 2p triple, |pₓ⟩ = (|1,−1⟩ − |1,1⟩)/√2 and |p_y⟩ = i(|1,−1⟩ + |1,1⟩)/√2 with |p_z⟩ = |1,0⟩ — a unitary change of basis in a three-dimensional subspace, so three states before and three after. Each remains an eigenstate of H and of L² with unchanged eigenvalues, and each has the identical radial density P(r), because P depends on n and l only. What changes is L_z: pₓ is not an L_z eigenstate, and measuring L_z on it returns +ħ or −ħ with equal probability, giving ⟨L_z⟩ = 0 and ΔL_z = ħ. The physical consequence is a current. For |n l m⟩ the probability current is azimuthal, j_φ = (ħmₗ/μr sinθ)|ψ|², carrying an orbital moment μ_z = −mₗμ_B; the real orbitals, being real up to a global phase, carry no current and no moment. Which basis is right is decided by whatever breaks the degeneracy — a field along z makes the m states the eigenstates of the perturbation, a bond or a crystal field along x makes the lobes.

Read the picture carefully. |ψ|² is the density for a position measurement on an ensemble of identically prepared atoms: the stippled cloud is ten thousand atoms measured once each, not one electron smeared into a fog. Yet −e|ψ|² is also exactly the mean charge density that sources the Hartree potential in a self-consistent field and that X-ray scattering measures as a form factor, so the charge-cloud reading is legitimate for mean quantities and wrong for fluctuations. And nothing in the picture moves: ∂|ψ|²/∂t = 0 for any eigenstate, which is why an atom in an orbital does not radiate, while a superposition of two levels has a density oscillating at ω = ΔE/ħ, and that is what makes a dipole radiate. Numerically the whole topic is one script: build the radial Hamiltonian on a grid with u(0) = 0, take eigenvectors from scipy.linalg.eigh_tridiagonal, normalise by trapezoid, and plot u². Then check three diagnostics — node count n − l − 1, ⟨1/r⟩ = Z/n²a₀, and ⟨T⟩ = −E. A box that is too small truncates the tail, which raises E and shrinks ⟨r⟩, and it bites the high-n states first.