# The Radial Equation & the Hydrogen Atom
> The angular problem is solved once and forever; what remains is one radial equation per
> potential. A single substitution, $u = rR$, turns it into a one-dimensional Schrödinger equation
> — every tool from Course P.5 applies — and for the Coulomb potential the series solution
> terminates only at special energies: $-13.6\ \text{eV}/n^2$, Bohr's lucky guess, now *derived*.
> This is the calculation that convinced physics that wave mechanics was true. Approach it the way
> a cat approaches a full food bowl: unhurried, but missing nothing.
## Learning Objectives
After this lesson you will be able to:
1. Reduce a central-potential problem to a 1D equation for $u(r) = rR(r)$ with the **effective potential** $V_{\text{eff}}(r)$ and its **centrifugal barrier**.
2. Apply the correct **boundary conditions** on $u$ and solve the infinite spherical well for $\ell = 0$.
3. Carry out the **full hydrogen solution**: asymptotics, series, the **recursion relation**, and termination.
4. Derive the **Bohr energies** $E_n = -13.6\ \text{eV}/n^2$, the **Bohr radius** $a_0$, and the $n^2$ **degeneracy**.
5. Work with the radial functions $R_{n\ell}$ and the **radial probability density** $P(r) = r^2|R|^2$, computing most probable and mean radii.
6. Interpret the quantum numbers $(n,\ell,m)$, use **spectroscopic notation**, and contrast the Bohr and Schrödinger pictures of the atom.
7. Solve the **finite spherical well** for $\ell = 0$ and explain why, unlike in one dimension, a three-dimensional well needs a minimum depth to bind.
8. Separate the **two-body problem** into centre-of-mass and relative motion, and correct hydrogen's spectrum with the **reduced mass**.
---
## Intuition
Classically, conservation of angular momentum flattens any central-force problem into
one-dimensional radial motion in an *effective* potential $V(r) + L^2/2mr^2$ — the orbiting piece
of the kinetic energy masquerades as a repulsive barrier
([P.1.3 Lagrangian Mechanics](../c01-theoretical-mechanics/03-lagrangian-mechanics.md)). Quantum
mechanics repeats the maneuver verbatim, with $L^2 \to \hbar^2\ell(\ell+1)$: the electron with
angular momentum is pushed away from the origin by the same centrifugal effect that keeps a
whirling ball on a string taut.
So the hydrogen atom is, at heart, a Course P.5 problem: a particle in a 1D well (Coulomb
attraction plus centrifugal repulsion), with bound states at discrete energies. The magic is in the
numbers that come out — energies matching every observed spectral line, and a natural length scale,
the Bohr radius, assembled from $\hbar$, $m$, and $e$ with no free parameters.
---
## Theory
### From 3D to effective 1D: the radial equation
Separation in [P.6.1](01-central-potentials-angular-momentum.md) left the radial equation
$$
\frac{d}{dr}\!\left(r^2\frac{dR}{dr}\right) - \frac{2mr^2}{\hbar^2}\bigl[V(r) - E\bigr]R = \ell(\ell+1)\,R .
$$
Substitute $u(r) \equiv r\,R(r)$. Then $R = u/r$, $\dfrac{dR}{dr} = \dfrac{u'}{r} - \dfrac{u}{r^2}$, so
$r^2\dfrac{dR}{dr} = r u' - u$ and $\dfrac{d}{dr}\bigl(r u' - u\bigr) = r u''$. The equation becomes
$r u'' - \frac{2mr}{\hbar^2}(V - E)\,u = \ell(\ell+1)\,\frac{u}{r}$, i.e.
$$
-\frac{\hbar^2}{2m}\frac{d^2u}{dr^2} + \underbrace{\left[V(r) + \frac{\hbar^2\ell(\ell+1)}{2mr^2}\right]}_{V_{\text{eff}}(r)} u = E\,u .
$$
This is *exactly* a 1D Schrödinger equation, with $V$ replaced by the **effective potential**: the
true potential plus the **centrifugal barrier** $\hbar^2\ell(\ell+1)/2mr^2$ — the same
$L^2/2mr^2$ term as in the classical central-force problem
([P.1.3](../c01-theoretical-mechanics/03-lagrangian-mechanics.md)), quantized. Normalization also
cooperates: $\int|\psi|^2 d^3r = \int_0^\infty |R|^2 r^2 dr = \int_0^\infty |u|^2 dr = 1$.
### Boundary conditions: the 1D toolkit applies
Two conditions replace the free constants of the 1D problems:
- **At the origin: $u(0) = 0$.** If $u(0)\ne 0$ then $R \sim u(0)/r$, and $\nabla^2(1/r) = -4\pi\delta^3(\mathbf r)$
produces a delta function no finite $V$ can cancel — not a solution. (For $\ell > 0$ the
centrifugal barrier enforces even faster vanishing, as we'll see.)
- **At infinity: $u \to 0$** for a bound state ($E < V(\infty)$), just as in
[P.5.2 Step, Well & Barrier](../c05-bound-states/02-step-well-barrier.md).
With $u$ vanishing at both ends, all the machinery of Course P.5 — node counting, matching,
shooting-method numerics — transfers wholesale to three dimensions.
### Warm-up: the infinite spherical well
Let $V(r) = 0$ for $r < a$ and $\infty$ outside. For $\ell = 0$ the radial equation is just
$u'' = -k^2 u$ with $k = \sqrt{2mE}/\hbar$ — the 1D box. The condition $u(0) = 0$ selects
$u = A\sin(kr)$, and $u(a) = 0$ quantizes $k = n\pi/a$:
$$
E_n = \frac{n^2\pi^2\hbar^2}{2ma^2}, \qquad n = 1, 2, 3, \dots
$$
with $R_{n0} = A\sin(n\pi r/a)/r$, perfectly finite at the origin. For $\ell > 0$ the solutions
regular at $r=0$ are the spherical Bessel functions, $R \propto j_\ell(kr)$, and the energies come
from their zeros — Exercise E5 works out $\ell = 1$.
### The finite spherical well: a minimum depth to bind
Now let the walls be finite: $V(r) = -V_0$ for $r < a$ and $0$ outside, with $V_0 > 0$, and look for
$\ell = 0$ bound states, $-V_0 < E < 0$. Inside, $u'' = -k^2u$ with $k = \sqrt{2m(E + V_0)}/\hbar$, and
$u(0) = 0$ again selects $u = A\sin kr$. Outside, $u'' = \kappa^2u$ with $\kappa = \sqrt{-2mE}/\hbar$,
and normalizability selects $u = Be^{-\kappa r}$. Matching $u'/u$ at $r = a$ gives the quantization
condition $k\cot ka = -\kappa$.
That is exactly the condition for the *odd* states of the 1D finite well of half-width $a$
([P.5.2](../c05-bound-states/02-step-well-barrier.md)): with $u(0) = 0$, the $\ell = 0$ problem is a 1D
problem on the half-line with a node at the origin. And it changes the physics of binding. In one
dimension every attractive well, however shallow, holds at least one bound state, the even ground
state; the odd states need a minimum depth. Here only the odd ones exist. Writing
$k_0 = \sqrt{2mV_0}/\hbar$ for the largest possible $k$, a solution needs $\cot ka < 0$, so $ka$ must
exceed $\pi/2$; the first bound state appears only when $k_0a > \pi/2$, that is when
$V_0a^2 > \pi^2\hbar^2/8m$. Each further multiple of $\pi$ in $k_0a$ adds one more state, so the number
of $\ell = 0$ bound states is $N = \lfloor k_0a/\pi + \tfrac12\rfloor$. A three-dimensional well has to
be deep enough, or wide enough, to bind anything at all; for $\ell > 0$ the centrifugal barrier makes it
harder still.
The *deuteron*, the proton–neutron bound state, is the classic case. Model the nuclear attraction as a
square well of radius $a \approx 2$ fm. The threshold depth is then a few tens of MeV, and fitting the
measured binding energy of 2.22 MeV puts the well only modestly above it, about 34 MeV: deep enough for
one bound state, and nowhere near the depth a second would need, nine times the threshold. That is why the
deuteron is so weakly bound and has no excited states (E7).
### The hydrogen atom I: setting up
Now the main event. The proton's Coulomb potential is
$$
V(r) = -\frac{e^2}{4\pi\epsilon_0}\frac1r ,
$$
so the radial equation reads
$$
-\frac{\hbar^2}{2m}\frac{d^2u}{dr^2} + \left[-\frac{e^2}{4\pi\epsilon_0 r} + \frac{\hbar^2\ell(\ell+1)}{2mr^2}\right]u = E\,u .
$$
Bound states have $E < 0$. States with $E > 0$ are not bound: the electron can reach infinity, any
positive energy is allowed, and the spectrum there is *continuous*. This ionization continuum begins
exactly at the series limits of [P.2.2](../c02-radiation-old-quantum/02-atomic-models-and-spectra.md).
We want the discrete part. Define
$$
\kappa \equiv \frac{\sqrt{-2mE}}{\hbar}, \qquad \rho \equiv \kappa r, \qquad
\rho_0 \equiv \frac{me^2}{2\pi\epsilon_0\hbar^2\kappa} ,
$$
($\kappa$ is real; $\rho$ is dimensionless radius; $\rho_0$ packages the strength of the Coulomb
attraction relative to the energy). Dividing the equation by $E$ and tidying up:
$$
\frac{d^2u}{d\rho^2} = \left[1 - \frac{\rho_0}{\rho} + \frac{\ell(\ell+1)}{\rho^2}\right]u .
$$
### The hydrogen atom II: asymptotics and the series
**Large $\rho$:** the bracket $\to 1$, so $u'' \approx u$ and $u \sim Ae^{-\rho} + Be^{\rho}$;
normalizability kills $B$. **Small $\rho$:** the centrifugal term dominates,
$u'' \approx \dfrac{\ell(\ell+1)}{\rho^2}u$, with solutions $\rho^{\ell+1}$ and $\rho^{-\ell}$; the
second violates $u(0) = 0$ (and normalizability), so $u \sim C\rho^{\ell+1}$. Peel off both
behaviors:
$$
u(\rho) = \rho^{\ell+1}\,e^{-\rho}\,v(\rho) .
$$
Substituting (two applications of the product rule) turns the equation for $u$ into
$$
\rho\,\frac{d^2v}{d\rho^2} + 2(\ell + 1 - \rho)\,\frac{dv}{d\rho} + \bigl[\rho_0 - 2(\ell+1)\bigr]v = 0 .
$$
Try a power series $v(\rho) = \sum_{j=0}^\infty c_j\rho^j$. Then
$\rho\,v'' = \sum_j j(j+1)c_{j+1}\rho^j$, $2(\ell+1)v' = \sum_j 2(\ell+1)(j+1)c_{j+1}\rho^j$, and
$-2\rho v' = -\sum_j 2j\,c_j\rho^j$ (each after shifting the summation index). Collecting the
coefficient of $\rho^j$:
$$
j(j+1)c_{j+1} + 2(\ell+1)(j+1)c_{j+1} - 2j\,c_j + \bigl[\rho_0 - 2(\ell+1)\bigr]c_j = 0 ,
$$
which solves to the **recursion relation**
$$
c_{j+1} = \frac{2(j + \ell + 1) - \rho_0}{(j+1)(j + 2\ell + 2)}\;c_j .
$$
One constant $c_0$ (fixed later by normalization) generates the whole series.
### The hydrogen atom III: quantization
Suppose the series *never* terminates. For large $j$ the recursion is
$c_{j+1} \approx \dfrac{2}{j+1}c_j$, whose solution is $c_j \approx \dfrac{2^j}{j!}c_0$ — the
series of $c_0e^{2\rho}$. Then $u \approx c_0\,\rho^{\ell+1}e^{-\rho}e^{2\rho} = c_0\,\rho^{\ell+1}e^{+\rho}$:
the asymptotic disease we amputated grows back. The **only escape is termination**: some maximal
index $j_{\max}$ must give $c_{j_{\max}+1} = 0$, i.e.
$$
2(j_{\max} + \ell + 1) = \rho_0 .
$$
Define the **principal quantum number** $n \equiv j_{\max} + \ell + 1 \in \{1, 2, 3, \dots\}$, so
$\rho_0 = 2n$. Unpacking $\rho_0$'s definition, $\kappa = \dfrac{me^2}{4\pi\epsilon_0\hbar^2 n}$,
and $E = -\dfrac{\hbar^2\kappa^2}{2m}$ gives
$$
E_n = -\left[\frac{m}{2\hbar^2}\left(\frac{e^2}{4\pi\epsilon_0}\right)^2\right]\frac{1}{n^2}
= -\frac{13.6\ \text{eV}}{n^2},
\qquad n = 1, 2, 3, \dots
$$
— the **Bohr formula**, exactly, but derived from a real dynamical theory rather than postulated.
The natural length scale also assembles itself: $\kappa = 1/(a_0 n)$ with
$$
a_0 = \frac{4\pi\epsilon_0\hbar^2}{me^2} = 0.529\ \text{Å} = 5.29\times10^{-11}\ \text{m},
$$
the **Bohr radius**. Two structural facts come free. Since $j_{\max} = n - \ell - 1 \ge 0$,
$$
\ell \le n - 1 \qquad (\ell = 0, 1, \dots, n-1),
$$
and counting all states at energy $E_n$ ($2\ell+1$ values of $m$ for each allowed $\ell$):
$$
d(n) = \sum_{\ell=0}^{n-1}(2\ell + 1) = 2\cdot\frac{(n-1)n}{2} + n = n^2 .
$$
The $m$-degeneracy is guaranteed by spherical symmetry; the extra $\ell$-degeneracy (2s degenerate
with 2p) is an "accident" special to the pure $1/r$ potential.
### Two bodies, not one: the reduced mass
Everything above treated the proton as a fixed centre. It is not: electron and proton both move,
about their common centre of mass. The repair is exact. For two particles of masses $m_1, m_2$
interacting through $V(\lvert\mathbf r_1 - \mathbf r_2\rvert)$, switch to the *centre-of-mass* and
*relative* coordinates
$$
\mathbf R = \frac{m_1\mathbf r_1 + m_2\mathbf r_2}{M}, \qquad \mathbf r = \mathbf r_1 - \mathbf r_2,
\qquad M = m_1 + m_2 .
$$
The chain rule turns the kinetic energy into
$-\frac{\hbar^2}{2m_1}\nabla_1^2 - \frac{\hbar^2}{2m_2}\nabla_2^2 = -\frac{\hbar^2}{2M}\nabla_R^2 - \frac{\hbar^2}{2\mu}\nabla_r^2$
(E6), with the *reduced mass* $\mu = m_1m_2/(m_1 + m_2)$. Since the potential depends on $\mathbf r$
alone, the equation separates: $\Psi = \phi(\mathbf R)\,\psi(\mathbf r)$ with total energy
$E_{\text{cm}} + E$, where
$$
-\frac{\hbar^2}{2M}\nabla_R^2\,\phi = E_{\text{cm}}\,\phi,
\qquad
\Big[-\frac{\hbar^2}{2\mu}\nabla_r^2 + V(r)\Big]\psi = E\,\psi .
$$
The atom as a whole moves freely, as a plane wave of mass $M$, and the internal motion is that of a
single fictitious particle of mass $\mu$ in the potential $V$. Every result above therefore holds with
$m \to \mu$. For hydrogen $\mu = m_e/(1 + m_e/m_p) = 0.99946\,m_e$: the energies shrink by 0.054%, to
$E_1 = -13.598$ eV, the Bohr radius grows by the same fraction, and the Rydberg constant becomes
$R_H = R_\infty/(1 + m_e/m_p)$, the correction met in
[P.2.2](../c02-radiation-old-quantum/02-atomic-models-and-spectra.md). Spectroscopy resolves it easily:
deuterium's heavier nucleus moves its H-α line by 0.18 nm, which is how deuterium was discovered. In
positronium, an electron bound to a positron ([P.3.1](../c03-wave-particle-duality/01-light-as-particles.md)),
$\mu = m_e/2$: every energy halves, to $E_1 = -6.8$ eV, and every length doubles.
### The wavefunctions
The terminating polynomials $v(\rho)$ are (up to normalization) the **associated Laguerre
polynomials** $L_{n-\ell-1}^{2\ell+1}(2\rho)$, with $\rho = r/na_0$, so
$\psi_{n\ell m} = R_{n\ell}(r)\,Y_\ell^m(\theta,\varphi)$ with $R_{n\ell} = \rho^{\ell+1}e^{-\rho}v(\rho)/r$. The first few, normalized:
$$
R_{10} = \frac{2}{a_0^{3/2}}\,e^{-r/a_0}, \qquad
R_{20} = \frac{1}{\sqrt2\,a_0^{3/2}}\left(1 - \frac{r}{2a_0}\right)e^{-r/2a_0}, \qquad
R_{21} = \frac{1}{\sqrt{24}\,a_0^{3/2}}\,\frac{r}{a_0}\,e^{-r/2a_0}.
$$
The physically meaningful radial quantity is the **radial probability density**
$$
P(r) = r^2\,|R(r)|^2 ,
$$
the probability per unit radius of finding the electron in the spherical shell $[r, r+dr]$ (the
$r^2$ is the shell's growing area). For the ground state, $P(r) = \frac{4}{a_0^3}r^2e^{-2r/a_0}$;
setting $dP/dr = 0$ gives $2re^{-2r/a_0}(1 - r/a_0) = 0$, so the **most probable radius is exactly
$a_0$** — Bohr's orbit radius reborn as the peak of a distribution. The *mean* radius is larger,
because the distribution has a long tail:
$$
\langle r\rangle = \int_0^\infty r\,P(r)\,dr = \frac{4}{a_0^3}\int_0^\infty r^3 e^{-2r/a_0}dr
= \frac{4}{a_0^3}\cdot 3!\left(\frac{a_0}{2}\right)^4 = \frac32\,a_0 .
$$
> **Caution.** The electron is not "at the Bohr radius," and it is not on a trajectory at all —
> $P(r)\,dr$ is the probability of *finding* it in a shell upon measurement, per the Born rule
> ([P.4.1](../c04-schrodinger-equation/01-schrodinger-equation-born-rule.md)). And do not confuse
> $|R|^2$ with $r^2|R|^2$: the ground-state *density* $|\psi_{100}|^2 \propto e^{-2r/a_0}$ peaks at
> $r = 0$, while $P(r)$ peaks at $a_0$. Both statements are true; they answer different questions.
### Quantum numbers, notation, and Bohr vs Schrödinger
| Symbol | Name | Allowed values | Determines |
| --- | --- | --- | --- |
| $n$ | principal | $1, 2, 3, \dots$ | energy $E_n = -13.6\ \text{eV}/n^2$ |
| $\ell$ | orbital (azimuthal) | $0, 1, \dots, n-1$ | shape; $\lvert\mathbf L\rvert = \hbar\sqrt{\ell(\ell+1)}$ |
| $m$ | magnetic | $-\ell, \dots, +\ell$ | orientation; $L_z = m\hbar$ |
**Spectroscopic notation** labels $\ell = 0,1,2,3,\dots$ as $s, p, d, f,\dots$ prefixed by $n$: the
ground state is 1s; the $n=2$ level holds 2s and 2p; and so on.
| | Bohr model ([P.2.2](../c02-radiation-old-quantum/02-atomic-models-and-spectra.md)) | Schrödinger theory |
| --- | --- | --- |
| Energies | $-13.6\ \text{eV}/n^2$ (right, luckily) | same — now derived |
| Angular momentum | $L = n\hbar$; ground state $L = \hbar$ | $\lvert\mathbf L\rvert = \hbar\sqrt{\ell(\ell+1)}$; ground state $\ell = 0$, $L = 0$! |
| Picture | planetary orbits of definite radius | stationary probability clouds ("orbitals, not orbits") |
| States per level | one orbit | $n^2$ degenerate states $(n,\ell,m)$ |
The **Rydberg formula**, fitted to data in
[P.2.2](../c02-radiation-old-quantum/02-atomic-models-and-spectra.md), is now a theorem: a photon
emitted in the transition $n_i \to n_f$ carries $h\nu = E_{n_i} - E_{n_f}$, so
$$
\frac{1}{\lambda} = R_H\left(\frac{1}{n_f^2} - \frac{1}{n_i^2}\right), \qquad
R_H = \frac{m}{4\pi c\hbar^3}\left(\frac{e^2}{4\pi\epsilon_0}\right)^2 = 1.097\times10^7\ \text{m}^{-1}.
$$
Not every pair of levels actually radiates, though: coupling to the electromagnetic field enforces
the selection rule $\Delta\ell = \pm1$ (2p → 1s shines; 2s → 1s is forbidden and that state is
metastable). Deriving this requires time-dependent perturbation theory, beyond our scope here.
The Schrödinger hydrogen atom is exact for the Hamiltonian written down, but that Hamiltonian leaves
three things out. It has no **spin**, which [P.6.3](03-magnetic-fields-stern-gerlach-spin.md) adds and
which doubles every count to $2n^2$. It is non-relativistic: relativity and the coupling of spin to the
orbital motion split each level by a relative amount of order $\alpha^2 \approx 5\times10^{-5}$, the
**fine structure** that Sommerfeld fitted
([P.2.3](../c02-radiation-old-quantum/03-correspondence-and-limits.md)) and Dirac's equation explains.
And it treats the electromagnetic field classically: quantizing the field gives the excited states
their finite lifetimes (spontaneous emission) and adds the Lamb shift, about a gigahertz in the $n = 2$
level. Each is a small correction; the $-13.6\ \text{eV}/n^2$ skeleton derived here survives them all.
---
## Worked Examples
### Example 1 — Lyman-α from first principles
The $n = 2 \to 1$ transition emits $\Delta E = E_2 - E_1 = -3.40 - (-13.6) = 10.2\ \text{eV}$. The
wavelength: $\lambda = hc/\Delta E = (1240\ \text{eV·nm})/10.2\ \text{eV} = 121.6\ \text{nm}$ —
deep ultraviolet, matching the measured Lyman-α line ($121.567\ \text{nm}$) to our precision. The
same arithmetic on $3\to2$, $4\to2$, ... reproduces the visible Balmer series a century of
spectroscopists catalogued.
### Example 2 — Where is the ground-state electron?
What is the probability of finding it *beyond* the Bohr radius? With $x = 2r/a_0$:
$$
P(r > a_0) = \frac{4}{a_0^3}\int_{a_0}^\infty r^2 e^{-2r/a_0}\,dr
= \frac12\int_2^\infty x^2 e^{-x}\,dx
= \frac12\Bigl[-(x^2 + 2x + 2)e^{-x}\Bigr]_2^\infty
= \frac{10\,e^{-2}}{2} = 5e^{-2} \approx 0.677 .
$$
Two-thirds of measurements find the electron outside the "orbit" it supposedly occupies. The Bohr
radius marks the peak of $P(r)$, not a boundary — the atom has no edge, only an exponential fade.
---
## Hands-on (Python)
**Shooting method** on $u(r)$, in Bohr-radius units where the equation becomes
$u'' = [\ell(\ell+1)/x^2 - 2/x - \varepsilon]u$ with $\varepsilon = E/(13.6\ \text{eV})$:
```python
import numpy as np
def shoot(eps, l, x_max=40.0, n=8000):
"""Integrate u'' = k(x) u outward; return u(x_max). Eigenvalue <=> u(x_max) ~ 0."""
x = np.linspace(1e-4, x_max, n); h = x[1] - x[0]
k = l*(l + 1)/x**2 - 2.0/x - eps
u = np.zeros(n); u[0], u[1] = x[0]**(l + 1), x[1]**(l + 1) # u ~ x^{l+1} near 0
for i in range(1, n - 1):
u[i+1] = 2*u[i] - u[i-1] + h*h*k[i]*u[i]
return u[-1]
def eigenvalue(l, lo, hi):
for _ in range(60): # bisect on the sign of u(x_max)
mid = 0.5*(lo + hi)
if shoot(lo, l)*shoot(mid, l) < 0: hi = mid
else: lo = mid
return 0.5*(lo + hi)
for l, brackets in [(0, [(-1.5, -0.6), (-0.35, -0.15), (-0.14, -0.09)]),
(1, [(-0.35, -0.15), (-0.14, -0.09)])]:
print(l, [round(eigenvalue(l, *b)*13.6057, 2) for b in brackets])
# Expect: l=0 -> [-13.61, -3.40, -1.51]; l=1 -> [-3.40, -1.51]
# E depends only on n -- the accidental Coulomb degeneracy, found numerically.
```
**Radial probability densities** $P(r) = r^2R^2$ in units of $a_0$:
```python
import matplotlib.pyplot as plt
r = np.linspace(0, 25, 1000) # r in units of a0
R = {(1,0): 2*np.exp(-r),
(2,0): (1 - r/2)*np.exp(-r/2)/np.sqrt(2),
(2,1): r*np.exp(-r/2)/np.sqrt(24),
(3,0): (2/np.sqrt(27))*(1 - 2*r/3 + 2*r**2/27)*np.exp(-r/3)}
for (n, l), Rnl in R.items():
plt.plot(r, r**2*Rnl**2, label=f"$P_{{{n}{l}}}$")
plt.xlabel("$r/a_0$"); plt.ylabel("$P(r)\\,a_0$"); plt.legend(); plt.show()
# Expect: P_10 peaks at r = a0; P_20 has a node (2 bumps); P_21 one bump peaking
# at 4 a0; P_30 has two nodes. Check each curve integrates to 1 with np.trapezoid.
```
**Monte Carlo orbital clouds** — rejection-sample $|\psi|^2$ and scatter the points:
```python
rng = np.random.default_rng(7)
def psi2_210(x, y, z): # |psi_210|^2 ∝ e^{-r} z^2 (a0 = 1)
return np.exp(-np.sqrt(x*x + y*y + z*z))*z**2
def psi2_320(x, y, z): # |psi_320|^2 ∝ r^4 e^{-2r/3}(3cos^2θ - 1)^2
r = np.sqrt(x*x + y*y + z*z)
return r**4*np.exp(-2*r/3)*(3*(z/np.maximum(r, 1e-12))**2 - 1)**2
def sample(pdf, half_box, f_max, n_keep=4000):
kept = np.empty((0, 3))
while len(kept) < n_keep: # rejection sampling in a cube
p = rng.uniform(-half_box, half_box, (100_000, 3))
acc = p[rng.uniform(0, f_max, 100_000) < pdf(*p.T)]
kept = np.vstack([kept, acc])
return kept[:n_keep]
fig = plt.figure(figsize=(11, 5))
for i, (pdf, box, fmax, name) in enumerate(
[(psi2_210, 15, 0.6, "2p$_z$"), (psi2_320, 30, 95.0, "3d$_{z^2}$")]):
pts = sample(pdf, box, fmax)
ax = fig.add_subplot(1, 2, i + 1, projection='3d')
ax.scatter(*pts.T, s=1, alpha=0.3); ax.set_title(name)
plt.show()
# Expect: the 2p_z dumbbell along z; the 3d_z2 dumbbell wearing a donut at the
# waist -- the textbook orbital pictures, generated from |psi|^2 itself.
```
---
## Exercises
**E1 (easy).** Compute the wavelength of the Balmer H-α line ($n = 3 \to 2$). Which part of the spectrum is it in?
Solution
$\Delta E = 13.6\left(\frac14 - \frac19\right) = 13.6\cdot\frac{5}{36} = 1.89\ \text{eV}$;
$\lambda = 1240/1.89 = 656\ \text{nm}$ — red, visible; the famous H-α line of astronomy.
Solution
$\int_0^\infty \frac{4}{a_0^3}e^{-2r/a_0}r^2dr = \frac{4}{a_0^3}\cdot\frac{2!}{(2/a_0)^3} = \frac{4}{a_0^3}\cdot\frac{2a_0^3}{8} = 1$
✓ (using $\int_0^\infty r^ne^{-br}dr = n!/b^{n+1}$).
Solution
$P(r) \propto r^2R_{21}^2 \propto r^4e^{-r/a_0}$. Then $dP/dr \propto (4r^3 - r^4/a_0)e^{-r/a_0} = 0$
gives $r = 4a_0$ — exactly Bohr's $n^2a_0$ for $n = 2$. The circular Bohr orbit survives as the
peak of the maximal-$\ell$ radial distribution; small-$\ell$ states (like 2s) have no classical
counterpart at all.
Solution
$\langle 1/r\rangle = \frac{4}{a_0^3}\int_0^\infty r\,e^{-2r/a_0}dr = \frac{4}{a_0^3}\cdot\frac{a_0^2}{4} = \frac{1}{a_0}$.
So $\langle V\rangle = -\frac{e^2}{4\pi\epsilon_0}\langle 1/r\rangle = -\frac{e^2}{4\pi\epsilon_0 a_0} = -27.2\ \text{eV} = 2E_1$,
and $\langle T\rangle = E_1 - \langle V\rangle = -13.6 + 27.2 = +13.6\ \text{eV} = -E_1$ ✓ — the
quantum virial theorem for a $1/r$ potential, $2\langle T\rangle = -\langle V\rangle$.
Solution
With $u = \sin(kr)/(kr) - \cos(kr)$ (this is $kr\,j_1(kr)$), differentiating twice gives
$u'' = -\frac{k\sin kr}{r} - \frac{2\cos kr}{r^2} + \frac{2\sin kr}{kr^3} + k^2\cos kr$, while
$\bigl(k^2 - \frac{2}{r^2}\bigr)u = \frac{k\sin kr}{r} - k^2\cos kr - \frac{2\sin kr}{kr^3} + \frac{2\cos kr}{r^2}$;
the sum vanishes term by term, so $u'' + \bigl[k^2 - 2/r^2\bigr]u = 0$ — the $\ell = 1$ radial
equation with $V = 0$. As $r\to0$, $\sin(kr)/(kr) \to 1 - (kr)^2/6$ and $\cos(kr)\to 1 - (kr)^2/2$,
so $u \to (kr)^2/3 \to 0$ ✓ (indeed $u\sim r^{\ell+1}$). The wall imposes $u(a) = 0$:
$\tan(ka) = ka$, whose smallest positive root is $ka \approx 4.4934$. Hence
$E_{1,\ell=1} = (4.4934)^2\hbar^2/2ma^2 \approx 20.2\,\hbar^2/2ma^2$, about $2.05\times$ the
$\ell = 0$ ground energy $\pi^2\hbar^2/2ma^2 \approx 9.87\,\hbar^2/2ma^2$ — the centrifugal barrier
raises every $\ell = 1$ level, and unlike hydrogen there is no accidental degeneracy.
Solution
(a) $\mathbf r_1$ enters $\mathbf R$ with weight $m_1/M$ and $\mathbf r$ with weight $+1$, so the chain
rule gives $\nabla_1 = \frac{m_1}{M}\nabla_R + \nabla_r$; likewise $\nabla_2 = \frac{m_2}{M}\nabla_R - \nabla_r$.
Squaring, $\frac{1}{m_1}\nabla_1^2 = \frac{m_1}{M^2}\nabla_R^2 + \frac2M\nabla_R\cdot\nabla_r + \frac{1}{m_1}\nabla_r^2$
and $\frac{1}{m_2}\nabla_2^2 = \frac{m_2}{M^2}\nabla_R^2 - \frac2M\nabla_R\cdot\nabla_r + \frac{1}{m_2}\nabla_r^2$.
The cross terms cancel, and adding gives $\frac{1}{M}\nabla_R^2 + \big(\frac1{m_1} + \frac1{m_2}\big)\nabla_r^2$,
with $\frac1{m_1} + \frac1{m_2} = \frac1\mu$ ✓. (b) Lyman-α carries $\tfrac34$ of the Rydberg energy times
$\mu/m_e$. Hydrogen: $\tfrac34\times13.606\times(1 + 1/1836.2)^{-1} = 10.199$ eV, so
$\lambda = 1239.84/10.199 = 121.57$ nm, the measured value; with a fixed nucleus it would be 121.50 nm.
Positronium: $\mu/m_e = \tfrac12$, so 5.102 eV and $\lambda = 243.0$ nm, in the near ultraviolet.
Solution
(a) $V_0 > \pi^2(\hbar c)^2/8\mu c^2a^2 = 9.87\times(197.3)^2/(8\times469.5\times4.41) = 23$ MeV.
(b) $\kappa = \sqrt{2\mu c^2B}/\hbar c = \sqrt{2\times469.5\times2.22}/197.3 = 0.232\ \mathrm{fm^{-1}}$
and $k = \sqrt{2\mu c^2(V_0 - B)}/\hbar c$. Try $V_0 = 34$ MeV: $k = 0.875\ \mathrm{fm^{-1}}$,
$ka = 1.838$, $k\cot ka = -0.240\ \mathrm{fm^{-1}}$, slightly too negative. At $V_0 = 33.7$ MeV:
$k = 0.871$, $ka = 1.830$, $k\cot ka = -0.231\ \mathrm{fm^{-1}}$ ✓. So $V_0 \approx 34$ MeV, about 1.5
times the threshold. Note $ka \approx 1.83$, just past $\pi/2 = 1.57$: the wavefunction barely turns
over inside the well, and most of the probability lies *outside* it, in the tail $e^{-\kappa r}$ of
range $1/\kappa \approx 4.3$ fm. (c) A second state needs $k_0a > 3\pi/2$, i.e. nine times the
threshold depth, about 210 MeV. The real deuteron has no excited state, consistent with a well of
about 34 MeV.
Answers
1. $V_{\text{eff}} = V(r) + \hbar^2\ell(\ell+1)/2mr^2$. The barrier is the angular kinetic energy: classically $L^2/2mr^2$ from conservation of $\mathbf L$; quantum mechanically the $\hat L^2/2mr^2$ term acting on $Y_\ell^m$.
2. Otherwise $R\sim1/r$ and $\nabla^2\psi$ contains $-4\pi\delta^3(\mathbf r)$, which nothing in the equation can balance; the wavefunction would not solve the Schrödinger equation at the origin.
3. Asymptotics force $u = \rho^{\ell+1}e^{-\rho}v$; the series for $v$ obeys $c_{j+1} = \frac{2(j+\ell+1)-\rho_0}{(j+1)(j+2\ell+2)}c_j$; a non-terminating series grows like $e^{2\rho}$ and ruins normalizability; termination forces $\rho_0 = 2n$, which converts to $E_n \propto -1/n^2$.
4. $d(n) = \sum_{\ell=0}^{n-1}(2\ell+1) = n^2$. The $2\ell+1$ $m$-degeneracy follows from rotational symmetry; the degeneracy *across different $\ell$* is accidental, special to the $1/r$ potential.
5. They answer different questions: $|\psi|^2$ is probability per unit *volume* (largest at the origin), $P(r) = 4\pi r^2|\psi|^2$ is probability per unit *radius*, and the growing shell area $\propto r^2$ moves its peak out to $a_0$.
6. For $\ell = 0$, $u = rR$ must vanish at the origin, so the radial problem is the 1D well restricted to *odd* solutions. The even 1D ground state, which exists for any depth, is not allowed; the first odd state needs $k_0a > \pi/2$, i.e. $V_0a^2 > \pi^2\hbar^2/8m$.
7. The two-body problem separates into free centre-of-mass motion and relative motion of mass $\mu = m_1m_2/(m_1 + m_2)$, so $m_e \to \mu$ everywhere: hydrogen's energies shrink by 0.054% and $R_H = R_\infty/(1 + m_e/m_p)$. Positronium has $\mu = m_e/2$, so its energies are halved ($E_1 = -6.8$ eV) and its lengths doubled.