University Physics V · The Hydrogen Atom · 11.4

Radial Wavefunctions & the Laguerre Solution

The angular half of hydrogen is settled for every central potential; the energies live entirely on the half-line. Here you scale the radius, divide out the behaviour at both ends, and let a recursion that refuses to terminate tell you which energies are allowed — then rebuild the same spectrum as a matrix and learn where the grid lies.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

At fixed l the hydrogen problem is one ordinary differential equation on the half-line, and everything about the spectrum is decided by how its solutions behave at the two ends. Scale the radius by the decay length 1/κ and the equation loses every constant except a single dimensionless number λ = Z/(κa₀); the origin then fixes the regular exponent ρ^(l+1) and infinity fixes the decaying factor e^(−ρ/2). Divide both out and what remains obeys Kummer's equation, whose power series has coefficient ratio (j + l + 1 − λ)/[(j+1)(j+2l+2)] → 1/j, the large-j ratio belonging to e^ρ.

So a generic λ produces a perfectly smooth solution that grows as e^(+ρ/2) and is not in L²(0,∞). Only when the recursion terminates, at λ = nᵣ + l + 1 ≡ n, does the decaying tail survive, and that one arithmetical accident is the whole of Eₙ = −13.606 Z²/n² eV. The cost is that this is an exact-solubility argument rather than a method: it works because 1/r is special, and it collapses the moment you screen the charge.

The general-purpose version is the matrix one — put the same operator on a grid with u(0) = 0 built in, and diagonalise — which reproduces the levels to five figures but pays with two errors that are not physics: an O(h²) stencil bias, and a wall at radius R that squeezes exactly the diffuse, high-n states you most wanted.

Simple definition: The radial wavefunction uₙₗ = r Rₙₗ is the half-line eigenfunction of the radial Hamiltonian at fixed l: an r^(l+1) rise at the origin, an e^{(−Zr/na₀)} tail, and between them an associated Laguerre polynomial of degree n − l − 1.
Example: For n = 3, l = 1 the polynomial has degree 3 − 1 − 1 = 1, so u₃₁ ∝ r²(1 − r/6a₀)e^(−r/3a₀): exactly one interior node, at r = 6a₀ = 0.318 nm, with E₃ = −13.606/9 = −1.512 eV.

Radial equation in the scaled variabled²u/dρ² = [1/4 − λ/ρ + l(l+1)/ρ²] u

Every bound state solves this one equation. The energy is hidden inside the length 1/κ, so λ is the only knob left to turn.

ρ = 2κr with κ = √(−2μE)/ħ in m⁻¹; λ = Z/(κa₀) is dimensionless, a₀ = 4πε₀ħ²/μe²

Stripping both asymptoticsu(ρ) = ρ^(l+1) e^(−ρ/2) v(ρ)

What is left over can be a plain power series, because the two behaviours a series handles badly are already outside it.

ρ^(l+1) is the regular Frobenius root at ρ = 0; e^(−ρ/2) is the decaying root at ρ → ∞

The recursion for vc₍ⱼ₊₁₎ = (j + l + 1 − λ) / [(j+1)(j+2l+2)] · cⱼ

That limiting ratio is the one belonging to e^ρ, so an unterminated v drags u up as e^(+ρ/2) and straight out of L²(0,∞).

v = Σⱼ cⱼ ρ^j with c₀ ≠ 0; the ratio tends to 1/j as j → ∞

Termination fixes n and the spectrumλ = n = nᵣ + l + 1 → Eₙ = −(Z²/n²)⋅ħ²/2μa₀²

Quantisation here is a normalisability condition, not an extra postulate — and n is revealed as a radial count plus an angular one.

nᵣ = deg v = number of interior radial nodes; ħ²/2μa₀² = 13.606 eV, and l ≤ n − 1 follows

The normalised radial functionRₙₗ = Nₙₗ e^(−ρ/2) ρ^l L^(2l+1)₍n−l−1₎(ρ), ρ = 2Zr/na₀

The polynomial's n − l − 1 positive roots are exactly the radial nodes, so raising l at fixed n trades radial structure for angular.

Nₙₗ = √[(2Z/na₀)³ (n−l−1)! / (2n(n+l)!)], fixed by ∫₀^∞ |R|² r² dr = 1

Grid Hamiltonian for the matrix routeHⱼⱼ = ħ²/μh² + V_eff(rⱼ), H₍j, j±1₎ = −ħ²/2μh²

u(0) = 0 costs nothing — it is the omitted j = 0 row. Eigenvalues converge as h², and the wall biases any state reaching it.

rⱼ = jh for j = 1…N, Dirichlet wall at R = (N+1)h, so u₀ = u₍N+1₎ = 0

Scale the radius and one number survives

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.

Divide out both ends before writing a series

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.

The recursion, and the one thing that can go wrong

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.

What termination buys: n, nodes and degeneracy

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.

The same operator, this time as a matrix

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.

Two errors in the matrix answer that are not physics

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.

Change one variable at a time

Make the relationship visible.

Interactive model

Principal quantum number n3

Grid box radius R30 a₀

Set n = 4 and slide the wall out. At R = 40 a₀ the wall amplitude still reads 0.0737, and a grid diagonalisation there returns −0.83136 eV instead of −0.85036 eV; at R = 60 a₀ the readout falls to 0.0033 and the level lands 5 × 10⁻⁵ eV from exact. Now set n = 1: the same 40 a₀ wall reads 0.0000.

ENERGY Eₙ-1.512 eV

RADIAL NODES n − 12

MEAN RADIUS ⟨r⟩13.5 a₀

WALL AMPLITUDE |u(R)|0.0250 a₀⁻¹⁄²

Live interpretationENERGY Eₙ: −1.512 eV. RADIAL NODES n − 1: 2. MEAN RADIUS ⟨r⟩: 13.5 a₀. WALL AMPLITUDE |u(R)|: 0.0250 a₀⁻¹⁄²

Catch the common trap

Explain before calculating.

Misconception check

The recursion is not truncated for convenience, and the unterminated series is not divergent. It converges for every ρ — that is the trap. What goes wrong is what it converges to: a ratio tending to 1/j is the signature of e^ρ, so v e^(−ρ/2) behaves as e^(+ρ/2) and u leaves L²(0,∞). Read the consequence carefully. The radial ODE has a smooth solution obeying u(0) = 0 at every negative E, so quantisation is not something the differential equation produces; it is the condition at infinity, imposed afterwards, that culls a continuum of solutions down to the integers λ = l + 1, l + 2, l + 3, …, which is to say n ≥ l + 1. The matrix route makes the point from the other side: a grid diagonalisation hands you discrete levels for whatever potential you feed it, because the wall at R has already imposed normalisability — which is why levels whose amplitude reaches that wall come out wrong.

In the l = 0 hydrogen radial problem, writing u = ρ e^(−ρ/2) v(ρ) gives series coefficients obeying c₍ⱼ₊₁₎/cⱼ = (j + 1 − λ)/[(j+1)(j+2)], a ratio that tends to 1/j as j grows. What does that limit establish?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyA hydrogen electron sits in the n = 4, l = 2 state. Give the degree of the polynomial that survives the recursion, the number of interior radial nodes, the radius of that node in nm, and the energy. Take a₀ = 0.052918 nm.

Degree first: nᵣ = n − l − 1 = 4 − 2 − 1 = 1, so v is linear and u₄₂ carries exactly one interior node.
Run the recursion c₍ⱼ₊₁₎ = (j + l + 1 − n)/[(j+1)(j+2l+2)] cⱼ with l = 2, n = 4 and c₀ = 1. c₁ = (0 + 3 − 4)/(1 × 6) = −1/6, then c₂ = (1 + 3 − 4)/(2 × 7) × c₁ = 0. So v(ρ) = 1 − ρ/6.
The node is where v vanishes: ρ = 6. With ρ = 2Zr/na₀ = r/2a₀ at Z = 1, n = 4, that is r = 12a₀ = 12 × 0.052918 = 0.635 nm.
Energy depends on n alone: E₄ = −13.606/4² = −0.850 eV. Sanity check the labels: l ≤ n − 1 reads 2 ≤ 3 ✓, and ⟨r⟩ = (a₀/2)(3n² − l(l+1)) = (a₀/2)(48 − 6) = 21a₀ = 1.11 nm.

AnswerDegree 1, one interior node at r = 12a₀ = 0.635 nm, E₄ = −0.850 eV, with ⟨r⟩ = 21a₀ = 1.11 nm.

MediumRun the recursion for the l = 0, λ = 3 state of hydrogen. Show that it terminates, write u₃₀(r) up to normalisation, and locate both radial nodes in units of a₀ and in nm.

With l = 0 the recursion collapses to c₍ⱼ₊₁₎ = (j + 1 − λ)/[(j+1)(j+2)] cⱼ. Put λ = 3 and c₀ = 1.
c₁ = (0 + 1 − 3)/(1 × 2) = −1. c₂ = (1 + 1 − 3)/(2 × 3) × (−1) = +1/6. c₃ = (2 + 1 − 3)/(3 × 4) × 1/6 = 0 — the ladder stops at degree 2, so nᵣ = 2 and n = nᵣ + l + 1 = 3.
So v(ρ) = 1 − ρ + ρ²/6 with ρ = 2Zr/na₀ = 2r/3a₀, giving u₃₀ ∝ r(1 − 2r/3a₀ + 2r²/27a₀²)e^(−r/3a₀), the standard 3s radial function.
Nodes: multiply v by 6 to get ρ² − 6ρ + 6 = 0, so ρ = 3 ± √3 = 1.2679 and 4.7321.
Convert with r = (3a₀/2)ρ: r = 1.902a₀ and 7.098a₀, i.e. 0.1006 nm and 0.3756 nm. Two nodes, exactly as n − l − 1 = 2 demands, and E₃ = −13.606/9 = −1.512 eV.

Answerv = 1 − ρ + ρ²/6 terminates at degree 2, so λ = n = 3 and E₃ = −1.512 eV; the nodes sit at r = 1.902a₀ = 0.1006 nm and r = 7.098a₀ = 0.3756 nm.

HardDiagonalise the l = 0 radial Hamiltonian on a uniform grid rⱼ = jh, j = 1…N, with a Dirichlet wall at R = (N+1)h. In atomic units (ħ = m = e = 1, 1 Ha = 27.211386 eV) a run at h = 0.01a₀ and R = 40a₀ returns −13.60535, −3.40140, −1.51170 and −0.83136 eV. Say which of these you trust, and repair the rest.

The matrix is symmetric tridiagonal: Hⱼⱼ = 1/h² − 1/rⱼ and H₍j, j±1₎ = −1/(2h²) in hartree. At h = 0.01 the diagonal carries 1/h² = 10⁴ Ha = 2.72 × 10⁵ eV, four orders of magnitude above even the ground state, so every level is a cancellation between large numbers: use scipy.linalg.eigh_tridiagonal rather than a dense solver.
Compare with Eₙ = −13.60569/n²: −13.60569, −3.40142, −1.51174, −0.85036 eV. The first three agree to 3.4 × 10⁻⁴ eV or better. The fourth is out by 0.019 eV — 2.2%, and less bound, which is the direction a squeezing wall pushes.
Separate the two error sources. The stencil is O(h²): rerunning at h = 0.02a₀ gives E₁ = −13.60433 eV, an error four times larger (1.36 × 10⁻³ against 3.40 × 10⁻⁴ eV). Richardson extrapolation gives (4 × (−13.60535) + 13.60433)/3 = −13.60569 eV, matching the exact value to six figures.
The n = 4 error is not h. Its mean radius is ⟨r⟩ = (a₀/2)(3n² − l(l+1)) = 24a₀, so a wall at 40a₀ stands at just 1.7⟨r⟩ and the true u still has amplitude 0.0737 a₀^(−1/2) there. Rerun at R = 80a₀ with the same h: E₄ = −0.850354 eV, now 2 × 10⁻⁶ eV from exact.
Rule of thumb: keep R ≳ 3⟨r⟩ ≈ 4.5n²a₀/Z for the highest level you want, verify each state has n − l − 1 nodes, and discard every eigenvalue near or above zero — those are box states, a discretised continuum, not bound levels.

Answern = 1–3 are converged; E₄ = −0.83136 eV is a box artefact, 2.2% high. Richardson removes the O(h²) bias (−13.60569 eV), and moving the wall to R = 80a₀ fixes E₄ to −0.85035 eV.