Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Explain why the hydrogen atom is the only atom for which the Schrödinger equation can be solved exactly — and why this one solution sets the language for all of chemistry.
Derive the separation $\psi(r,\theta,\varphi) = R(r)\, Y_\ell^m(\theta,\varphi)$ from spherical symmetry of the Coulomb potential.
Identify the three quantum numbers $(n, \ell, m)$ and their allowed ranges, and explain what each one physically "counts".
Compute hydrogen energy levels $E_n = -\dfrac{13.6\,\mathrm{eV}}{n^2}$ from first principles and connect them to the Lyman, Balmer, and Paschen spectral series.
Visualise the shapes of atomic orbitals (1s, 2s, 2p, 3d, …) and predict their number of radial and angular nodes.
Implement a numerical eigenvalue solver for the radial Schrödinger equation in PyTorch and verify that it reproduces the Bohr formula.

Why Hydrogen is the Rosetta Stone

"The hydrogen atom is the most important problem in modern physics. Solve it, and you understand every atom." — paraphrasing Wolfgang Pauli

Hydrogen is the simplest possible atom: one proton sitting at the origin, and one electron bound to it by the electric force. Two particles. That is it. There are no electron–electron repulsions to worry about, no inner shells to model, no shielding to approximate. Strip everything else away and the only physics left is one charge orbiting another.

Yet this single "trivial" example is the only atom in nature whose Schrödinger equation can be solved exactly, by pencil and paper, in closed form. Every other atom in the periodic table is solved numerically by perturbing around hydrogen, by combining hydrogen-like orbitals, or by approximating with hydrogen-shaped basis sets. The structure of the periodic table itself — the s, p, d, f blocks; the row lengths 2, 8, 18, 32; the "diagonal" chemistry of elements — is just hydrogen's quantum-number rules played out under the influence of many electrons.

The intuition in one sentence

Hydrogen is a one-body problem in a central potential; spherical symmetry collapses a 3D PDE into a 1D ODE, and the special algebra of the

1/r

Coulomb force gives discrete bound states whose energies depend only on a single integer

n

Historical Context

Hydrogen's spectrum was the puzzle that made quantum mechanics. The story is worth telling because it shows how an embarrassing pile of empirical observations forced physicists into a completely new view of nature.

1885: Balmer's mysterious formula

A Swiss school teacher, Johann Balmer, noticed that the four visible lines of hydrogen (red, cyan, blue-violet, deep violet) fit the empirical formula $\lambda = \dfrac{B\,n^2}{n^2 - 4}$ with $n = 3, 4, 5, 6$ . He had no idea why.

1913: Bohr's "atom of the impossible"

Niels Bohr fused Rutherford's nuclear model with Planck's quantum hypothesis to claim — without justification — that electrons only travel in orbits whose angular momentum is a multiple of $\hbar$ . Out fell the formula $E_n = -13.6\,\mathrm{eV}/n^2$ , which explained every hydrogen line. It was a miracle that worked for the wrong reasons.

1926: Schrödinger's right reasons

Erwin Schrödinger applied his new wave equation to the Coulomb problem and recovered Bohr's answer exactly — but now as the eigenvalues of a differential operator. Bohr's quantised orbits became standing waves on a 3D potential well. The answer was the same; the understanding was new.

Everything after

From hydrogen's solution flowed: the periodic table's block structure (Pauli + hydrogen-like orbitals), all of molecular orbital theory, the laser, magnetic resonance imaging, every density-functional calculation in modern chemistry. Hydrogen is literally the "Hello World" of quantum mechanics.

The Coulomb Potential

Two charges $+e$ and $-e$ separated by $r$ feel an electric potential energy

\displaystyle V(r) \;=\; -\frac{k_e\, e^2}{r}

where $k_e = 1/(4\pi\varepsilon_0)$ is the Coulomb constant. The minus sign is everything: the electron is bound, sitting in an infinitely deep well. The walls slope gently — they fall off like $1/r$ , which is slow enough that infinitely many bound states fit inside the well, getting denser and denser as they approach the top.

Intuition: a sloped funnel

Imagine a smooth, infinitely deep funnel whose walls flatten out as you climb. A marble at the bottom is held tightly. A marble higher up can still slosh around bound, but the rims are so close together near the top that an infinite number of distinct "sloshing modes" squeeze in just below the rim. Those modes are the bound states of hydrogen, and the rim is

E = 0

— the ionisation threshold.

The visualisation below plots $V(r)$ and overlays the discrete energy ladder $E_n = -13.6\,\mathrm{eV}/n^2$ . Slide n and watch the highlighted rung; the yellow dot marks the classical turning point $r_n = 2n^2 a_0$ where a classical particle of energy $E_n$ would stop.

Loading Coulomb potential viewer…

Read off the pattern

The rungs thin out exponentially toward E = 0 because of the 1/n² spacing.
Higher rungs live further from the proton: $r_n \propto n^2$ .
There is no lowest energy below $-13.6$ eV — quantum mechanics forbids the electron from spiralling into the proton.

Spherical Symmetry and Separation of Variables

Write the time-independent Schrödinger equation for one electron in the proton's field:

\displaystyle -\frac{\hbar^2}{2 m_e}\nabla^2 \psi(\mathbf{r}) - \frac{k_e e^2}{r}\,\psi(\mathbf{r}) \;=\; E\, \psi(\mathbf{r})

Because $V(r)$ depends only on the distance from the origin — not on direction — the natural coordinates are spherical: $(r, \theta, \varphi)$ . In those coordinates the Laplacian splits into a radial piece and an angular piece:

\displaystyle \nabla^2 = \underbrace{\frac{1}{r^2}\frac{\partial}{\partial r}\!\left(r^2 \frac{\partial}{\partial r}\right)}_{\text{radial}} + \underbrace{\frac{1}{r^2}\,\nabla^2_{\theta,\varphi}}_{\text{angular}}.

Because the potential talks only to $r$ and the angular Laplacian talks only to $(\theta, \varphi)$ , we can try a separable ansatz:

\displaystyle \psi(r, \theta, \varphi) \;=\; R(r)\, Y(\theta, \varphi)

Plug it in, divide by $R\,Y/r^2$ , and the equation cleanly separates into two pieces: an angular eigenproblem whose eigenvalue we'll call $\ell(\ell+1)$ , and a radial ODE that inherits that eigenvalue as a centrifugal term:

\displaystyle \boxed{\;-\nabla^2_{\theta,\varphi} Y \;=\; \ell(\ell+1)\,Y\;}

\displaystyle \boxed{\;-\frac{\hbar^2}{2 m_e}\frac{1}{r^2}\frac{d}{dr}\!\left(r^2 \frac{dR}{dr}\right) + \!\left[\,V(r) + \frac{\hbar^2 \ell(\ell+1)}{2 m_e r^2}\right]\! R \;=\; E\, R\;}

The bracketed quantity is the effective radial potential: the true Coulomb attraction plus a centrifugal repulsion that grows like $1/r^2$ at small $r$ . The centrifugal term is the quantum echo of the classical fact that a planet with non-zero angular momentum cannot fall straight into the sun.

The trick to remember

Separation of variables is just a rephrasing of the symmetry: since V respects rotations, the eigenfunctions can be chosen to also respect rotations. The angular eigenfunctions are universal — they appear in any 3D central-force problem (gravity, harmonic oscillator in 3D, hydrogen). Only the radial equation knows about the specific

1/r

The Angular Part: Spherical Harmonics

The angular eigenproblem $-\nabla^2_{\theta,\varphi}\, Y = \ell(\ell+1)\, Y$ has solutions called spherical harmonics $Y_\ell^m(\theta, \varphi)$ . They are labelled by two integers:

Quantum number	Allowed values	What it counts
ℓ (orbital)	0, 1, 2, 3, …	Total angular momentum: \|L\| = ℏ √ℓ(ℓ+1)
m (magnetic)	−ℓ, −ℓ+1, …, +ℓ (2ℓ+1 values)	z-component of angular momentum: L_z = ℏ m

The first few real spherical harmonics have shapes you have probably seen in chemistry class:

ℓ	m	Y (proportional to)	Picture
0	0	1	Uniform sphere (s-orbital)
1	0	cos θ (= z/r)	p_z — two lobes along z-axis
1	±1	sin θ cos φ, sin θ sin φ (= x/r, y/r)	p_x, p_y — two lobes along x or y
2	0	3 cos²θ − 1	d_{z²} — donut + lobes along z
2	±1, ±2	various combinations	d_{xy}, d_{xz}, d_{yz}, d_{x²−y²} — four-lobed cloverleaves

The number of angular nodes (flat planes or cones where $Y = 0$ ) equals $\ell$ . So an s-orbital has zero angular nodes (it's round), a p-orbital has one nodal plane (the equator for $p_z$ ), a d-orbital has two nodal surfaces, and so on. This is exactly what your 3D viewer below will show.

The Radial Equation

With $\ell$ chosen, the radial ODE becomes a 1D eigenvalue problem for $R(r)$ . After substituting $u(r) \equiv r\, R(r)$ (a standard trick that kills the awkward $r^2$ factor) it takes the deceptively familiar form

\displaystyle -\frac{\hbar^2}{2 m_e}\,\frac{d^2 u}{dr^2} \;+\; \underbrace{\!\left[\,V(r) + \frac{\hbar^2 \ell(\ell+1)}{2 m_e r^2}\right]\!}_{V_{\text{eff}}(r)}\, u \;=\; E\, u

— a 1D Schrödinger equation in $u(r)$ with boundary conditions $u(0) = 0$ (the wavefunction must be finite at the origin since the volume element is $r^2 dr$ ) and $u(r) \to 0$ as $r \to \infty$ (a bound state must die off).

The bound-state solutions are $u_{n,\ell}(r) = r\, R_{n,\ell}(r)$ with $R_{n,\ell}$ built from exponentials and Laguerre polynomials. The closed-form expression (with $r$ in units of the Bohr radius $a_0$ ) is

\displaystyle R_{n\ell}(r) \;=\; \sqrt{\left(\frac{2}{n}\right)^{3}\!\frac{(n - \ell - 1)!}{2n\,(n + \ell)!}}\; e^{-r/n}\,\left(\frac{2r}{n}\right)^{\!\ell}\, L^{2\ell+1}_{n - \ell - 1}\!\!\left(\frac{2r}{n}\right).

Don't panic at the size of this formula — every piece earns its keep:

$e^{-r/n}$ — the exponential decay at infinity. Higher $n$ means slower decay, larger orbital.
$r^\ell$ — suppresses the wavefunction near the origin for $\ell \ge 1$ ; this is the centrifugal barrier saying "no angular momentum near r = 0".
$L^{2\ell+1}_{n-\ell-1}$ — a polynomial of degree $n - \ell - 1$ ; its zeros are exactly the radial nodes of the wavefunction.
The square-root prefactor enforces $\int_0^\infty r^2 |R|^2 \, dr = 1$ .

Use the explorer below to feel how these pieces combine. Drag $n$ and $\ell$ ; watch the yellow line track the most-probable radius:

Loading radial wavefunction explorer…

Pattern to look for

For every state, the radial probability

P(r) = r^2|R_{n\ell}|^2

has

n - \ell

bumps. The outermost bump always sits near

r \approx n^2 a_0

— confirming Bohr's original intuition about orbital size.

Energy Levels: Where the −13.6 eV Comes From

The miracle of the $1/r$ potential is that the radial equation's eigenvalues are dictated by a single integer:

\displaystyle E_n \;=\; -\frac{m_e\, k_e^2\, e^4}{2\,\hbar^2}\,\frac{1}{n^{2}} \;=\; -\frac{13.6057\;\text{eV}}{n^{2}}, \qquad n = 1, 2, 3, \ldots

The combination $R_H \equiv m_e k_e^2 e^4 / (2 \hbar^2) \approx 13.606\;\text{eV}$ is the Rydberg energy. It is the ionisation energy of the hydrogen atom — the work you must do to pull the ground-state electron infinitely far from the proton.

Where do the units come from?

R_H

is the only energy you can build from the fundamental constants of the problem (electron mass

m_e

, the electromagnetic coupling

k_e e^2

, and Planck's constant

\hbar

) by dimensional analysis. There is nothing arbitrary about "13.6 eV" — change any of those constants and you get a different Rydberg.

Equally important is the Bohr radius, the natural length scale of the problem:

\displaystyle a_0 \;=\; \frac{\hbar^2}{m_e\, k_e\, e^2} \;\approx\; 0.5292\;\text{Å} \;=\; 5.292\times 10^{-11}\;\text{m}.

$a_0$ sets the size of the ground state; the n-th excited state grows as $n^2 a_0$ . A Rydberg atom with $n = 100$ is $10{,}000$ times bigger than the ground state — comparable to the size of a bacterial cell!

Why does E only depend on n?

The fact that $E$ ignores $\ell$ and $m$ — so e.g. 2s and 2p have the same energy — is called accidental degeneracy. It isn't really accidental: the $1/r$ Coulomb potential has an extra hidden symmetry (the SO(4) Runge–Lenz symmetry) that ties radial and angular pieces together. In any other central potential — say, a screened $e^{-r}/r$ — this degeneracy disappears and states with the same $n$ but different $\ell$ split apart. That splitting is exactly why the periodic table fills s < p < d rather than all-at-once.

Spectral Series and the Photon

When an electron drops from state $n_i$ to state $n_f$ the energy difference comes out as a single photon:

\displaystyle \hbar\omega \;=\; E_{n_i} - E_{n_f} \;=\; R_H\!\left(\frac{1}{n_f^{2}} - \frac{1}{n_i^{2}}\right).

Photons grouped by the same $n_f$ form a spectral series:

Series	n_f	Wavelength region	Famous lines
Lyman	1	UV (91–122 nm)	Lyman-α at 121.6 nm — brightest UV line in the cosmos
Balmer	2	Visible (365–656 nm)	Hα 656 nm (red), Hβ 486 nm, Hγ 434 nm
Paschen	3	Infrared (820–1875 nm)	Pα 1875 nm — important in stellar spectra
Brackett	4	Mid-IR (1.46–4.05 µm)	Brα at 4.05 µm

Hover any transition arrow below to see the predicted wavelength. Compare what you read off to a photo of a real hydrogen spectrum tube — Balmer lines line up to the nanometre.

Loading energy-level diagram…

Astronomy connection

Every photo you have ever seen of a glowing red nebula (the Orion Nebula, the Horsehead, the Eagle's "Pillars of Creation") is glowing in

H_\alpha

— the 656 nm Balmer-α line. The whole universe is dotted with this exact n = 3 → n = 2 transition.

Quantum Numbers (n, ℓ, m) — Three Knobs

Every bound state of hydrogen is uniquely labelled by three integers — three knobs you can turn independently within their allowed ranges:

Symbol	Range	Physical meaning	What it sets
n	1, 2, 3, …	Principal quantum number	Energy E_n and overall size n² a_0
ℓ	0, 1, …, n−1	Orbital angular momentum	Shape of orbital: \|L\| = ℏ √ℓ(ℓ+1)
m	−ℓ, −ℓ+1, …, +ℓ	Magnetic quantum number	Orientation: L_z = ℏ m

And there is a fourth — spin $m_s = \pm 1/2$ — which the non-relativistic Schrödinger equation does not produce, but Dirac's relativistic equation does. Spin is what doubles the capacity of each spatial orbital from 1 electron to 2, and is the reason the periodic table's rows have lengths 2, 8, 8, 18, …

The 'degeneracy' counting

For a given

n

, the number of distinct spatial states is

\sum_{\ell=0}^{n-1}(2\ell + 1) = n^2

. Multiply by 2 for spin and you get the shell capacities

2n^2

: K-shell (n=1) holds 2, L-shell (n=2) holds 8, M-shell (n=3) holds 18, N-shell (n=4) holds 32 — exactly the periods of the periodic table.

Atomic Orbitals: 1s, 2s, 2p, 3d, …

An orbital is shorthand for the 3D wavefunction $\psi_{n\ell m}(r,\theta,\varphi) = R_{n\ell}(r)\, Y_\ell^m(\theta,\varphi)$ . The chemistry names use a letter for $\ell$ :

ℓ	Letter	Why this letter
0	s	"sharp" — historical spectroscopic naming
1	p	"principal"
2	d	"diffuse"
3	f	"fundamental"
≥ 4	g, h, i, …	alphabetical from f

The number of nodes in the 3D wavefunction breaks down cleanly:

\text{total nodes} = n - 1, \qquad \text{radial} = n - \ell - 1, \qquad \text{angular} = \ell.

Drag the sliders below to walk through orbitals. Cyan and magenta encode the sign of $\psi$ so you can see the angular nodal planes as boundaries between the two colours. The yellow dot is the nucleus.

Loading 3D orbital viewer…

What to try

(1, 0, 0) — 1s: a single round cyan cloud, no nodes anywhere.
(2, 0, 0) — 2s: a round magenta core nested inside a cyan shell — that gap is the radial node.
(2, 1, 0) — 2p_z: two lobes along the z-axis, cyan top, magenta bottom — the xy-plane is the angular node.
(3, 2, 0) — 3d_z²: a peanut along z plus an equatorial torus — two angular nodal cones.
(3, 2, ±2) — d_{xy} / d_{x²−y²}: the iconic four-lobed cloverleaf.

Worked Example: The 1s Ground State by Hand

It is worth doing one full computation on paper — once. The 1s state is the cleanest possible target.

▶ Click to work through the 1s ground state step-by-step

Step 1 — Write down the ansatz

For $n = 1, \ell = 0, m = 0$ the angular harmonic is constant $Y_0^0 = 1/\sqrt{4\pi}$ . The radial ansatz with no nodes is the simplest decaying exponential:

R_{1,0}(r) = N\, e^{-r/a_0}.

Step 2 — Normalise

We need $\int_0^\infty R^2 r^2\, dr = 1$ . Using the standard integral $\int_0^\infty r^2 e^{-2 r / a_0}\, dr = \frac{a_0^3}{4}$ :

N^2 \cdot \frac{a_0^3}{4} = 1 \;\Longrightarrow\; N = 2\, a_0^{-3/2}.

So $R_{1,0}(r) = 2 a_0^{-3/2} e^{-r/a_0}$ , matching the closed form in the section above.

Step 3 — Plug into the radial equation

With $\ell = 0$ the centrifugal term is zero and the radial equation reduces to

-\frac{\hbar^2}{2 m_e} \frac{1}{r^2}\frac{d}{dr}\!\left(r^2 \frac{dR}{dr}\right) - \frac{k_e e^2}{r} R = E\, R.

Compute $dR/dr = -R/a_0$ and $(1/r^2)(d/dr)(r^2 dR/dr) = R \left(\frac{1}{a_0^2} - \frac{2}{a_0 r}\right)$ . Substitute back:

-\frac{\hbar^2}{2 m_e}\,R\!\left(\frac{1}{a_0^2} - \frac{2}{a_0 r}\right) - \frac{k_e e^2}{r}\,R = E\, R.

Step 4 — Cancel the 1/r terms

For the equation to hold for all $r$ , the coefficient of $1/r$ must vanish:

\frac{\hbar^2}{m_e a_0} - k_e e^2 = 0 \;\Longrightarrow\; a_0 = \frac{\hbar^2}{m_e k_e e^2}.

The Bohr radius drops out automatically. Plugging numbers gives $a_0 \approx 0.5292\;\text{Å}$ .

Step 5 — Read off the energy

What is left is $-\frac{\hbar^2}{2 m_e a_0^2}\, R = E\, R$ , so

E_1 = -\frac{\hbar^2}{2 m_e a_0^2} = -\frac{m_e k_e^2 e^4}{2 \hbar^2}.

Plug in $m_e = 9.109\times 10^{-31}\;\text{kg}$ , $e = 1.602\times 10^{-19}\;\text{C}$ , $\hbar = 1.055\times 10^{-34}\;\text{J·s}$ , $k_e = 8.988\times 10^{9}\;\text{N·m}^2/\text{C}^2$ :

E_1 = -2.180\times 10^{-18}\;\text{J} = -13.606\;\text{eV}.

Step 6 — Most probable and mean radius

The radial probability density is $P(r) = r^2 |R_{1,0}|^2 = \frac{4 r^2}{a_0^3} e^{-2r/a_0}$ .

Maximise: $dP/dr = 0 \Rightarrow r_* = a_0$ . So the most probable radius equals the Bohr radius — exactly the classical Bohr orbit radius. Wonderful coincidence.

Mean radius: $\langle r\rangle = \int_0^\infty r\, P(r)\, dr = \frac{4}{a_0^3} \int_0^\infty r^3 e^{-2 r/a_0}\, dr = \frac{4}{a_0^3}\cdot \frac{6 a_0^4}{16} = \frac{3 a_0}{2}.$ Larger than $r_*$ because the exponential tail extends well past the peak.

Sanity numbers

Quantity	Symbol	Value
Ground-state energy	E_1	−13.6057 eV
Bohr radius	a_0	0.5292 Å
Most-probable radius	r_*	a_0 = 0.5292 Å
Mean radius	⟨r⟩	(3/2) a_0 = 0.7937 Å
Ionisation energy	\|E_1\|	13.6 eV ≈ 109,737 cm⁻¹
Lyman-α wavelength (n=2→1)	λ_Lyα	121.567 nm

Python: Build the Wavefunctions Yourself

Theory is comforting; numbers are convincing. Here is a complete, self-contained Python program that (1) builds every hydrogen radial wavefunction from the closed-form, (2) verifies the energy ladder, (3) checks normalisation, (4) computes expectation values, and (5) plots the radial probability densities for the first six orbitals.

Read line-by-line: every card on the right explains why the corresponding line exists and what the numerical output looks like when you run it.

Hydrogen radial wavefunctions and energy levels

🐍hydrogen_radial.py

Explanation(11)

Code(75)

1Imports

numpy gives us vectorised arithmetic; matplotlib draws the curves; factorial is needed for the radial normalisation; genlaguerre returns a polynomial object representing L^{p}_q(x), the associated Laguerre polynomial that appears in every hydrogen radial wavefunction.

8Atomic units

We work in units where the Bohr radius a_0 = 1 and the Rydberg energy is 13.606 eV in magnitude. This is the cleanest way to write the formulas — every length is in a_0 and every energy is a multiple of the Rydberg.

12R_nl(r): the full radial wavefunction

This is the famous closed form derived in §29.5. ρ = 2r/n is the dimensionless radial variable that appears throughout. The normalization factor sqrt((2/n)^3 * (n-l-1)! / (2n * (n+l)!)) is fixed by requiring ∫ r² |R|² dr = 1 over the whole half-line.

EXECUTION STATE

n = principal quantum number (1, 2, 3, …)

l = orbital quantum number (0 … n−1)

r = radius in units of a_0 (scalar or numpy array)

20The associated Laguerre polynomial

L^{2l+1}_{n-l-1}(ρ) is the polynomial piece that determines the radial nodes. The degree n−l−1 equals the number of radial nodes — so 1s (n=1,l=0) has 0 nodes, 2s has 1 node, 3s has 2 nodes, and so on. scipy returns a poly1d object you can call like a function.

24energy_eV(n): the Bohr formula

The full derivation from the Schrödinger equation reproduces Bohr's original formula E_n = −13.6 eV / n². Notice that it depends only on n, not on l or m — this is the famous 'accidental' degeneracy of the 1/r potential.

28Radial probability density P(r)

The probability of finding the electron between radius r and r+dr in a thin spherical shell is P(r) dr = r² |R(r)|² dr. The r² factor comes from the volume of a thin shell (4π r² dr); the angular Y_lm part integrates to 1 over the sphere. P(r) is what you would 'see' if you took a long-exposure photo of the radial position of the electron.

33Expectation value ⟨r⟩

⟨r⟩ = ∫₀^∞ r · P(r) dr is the area-weighted average radius. Note that ⟨r⟩ is NOT the same as the most-probable r where P(r) peaks: for the 1s state P(r) peaks at r = a_0 but ⟨r⟩ = (3/2) a_0 because the long exponential tail drags the mean outward.

LOOP TRACE · 2 iterations

1s (n=1, l=0)

peak r* = 1.000 a_0 (Bohr radius)

<r> = 1.500 a_0

2s (n=2, l=0)

peak r* = ~5.2 a_0 (near outer node)

<r> = 6.000 a_0 = (3/2) · 2²

42Verifying the energy ladder

We loop n = 1..5 and print E_n. The output should be −13.606, −3.401, −1.512, −0.850, −0.544 eV — exactly the values Lyman, Balmer and friends measured from the hydrogen spectrum 40 years before Schrödinger wrote his equation.

LOOP TRACE · 5 iterations

n=1

E_1 = −13.606 eV (ground state)

n=2

E_2 = −3.401 eV = −13.606 / 4

n=3

E_3 = −1.512 eV = −13.606 / 9

n=4

E_4 = −0.850 eV = −13.606 / 16

n=5

E_5 = −0.544 eV = −13.606 / 25

47Sanity-check the normalization

If our formula and the normalisation constant are right, ∫₀^∞ r² |R_nl(r)|² dr should equal 1.0 for every (n,l). numpy's trapz approximates the integral by the trapezoidal rule — close to 1.000000 means we got the normalisation right.

53⟨r⟩ for s-states

There is a beautiful closed form: ⟨r⟩_{n,l} = (a_0/2)[3n² − l(l+1)]. For s-states (l=0) this collapses to (3/2) n². You should see exactly 1.5, 6.0, 13.5, 24.0 — matching the loop output to three decimal places.

59Plot the radial probability for six orbitals

The 2×3 panel shows 1s, 2s, 2p, 3s, 3p, 3d. Count the bumps: 1s has 1 hill, 2s has 2 hills separated by a node, 3s has 3 hills. The s-states (l=0) reach all the way down to r = 0 because the centrifugal barrier ℓ(ℓ+1)/r² is absent; the p, d, f states are pushed outward by that barrier and vanish at the origin.

64 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3from math import factorial
4from scipy.special import genlaguerre  # L_q^{p}(x)
5
6# ----- Physical constants (atomic units: a_0 = hbar = m_e = e = 1) -----
7A0_ANGSTROM = 0.52917721                  # Bohr radius in Angstrom
8RYDBERG_EV  = 13.605693                   # |E_1|  in eV
9
10
11def R_nl(n, l, r):
12    """Normalized radial wavefunction R_{nl}(r) with r in units of a_0.
13
14    R_{nl}(r) = sqrt( (2/n)^3 * (n-l-1)! / (2n * (n+l)!) )
15               * exp(-r/n) * (2r/n)^l * L_{n-l-1}^{2l+1}(2r/n)
16    """
17    rho = 2.0 * r / n
18    norm = np.sqrt((2.0 / n) ** 3 *
19                   factorial(n - l - 1) /
20                   (2 * n * factorial(n + l)))
21    laguerre = genlaguerre(n - l - 1, 2 * l + 1)
22    return norm * np.exp(-r / n) * rho**l * laguerre(rho)
23
24
25def energy_eV(n):
26    """Bound-state energy E_n in eV."""
27    return -RYDBERG_EV / n**2
28
29
30def radial_probability(n, l, r):
31    """P(r) = r^2 |R_{nl}(r)|^2  (already normalized so int P dr = 1)."""
32    R = R_nl(n, l, r)
33    return r * r * R * R
34
35
36def expectation_r(n, l, r_max=None, n_points=4000):
37    """Compute <r> by numerical integration of r * P(r)."""
38    if r_max is None:
39        r_max = 4 * n * n + 6
40    r = np.linspace(0.0, r_max, n_points)
41    P = radial_probability(n, l, r)
42    return np.trapz(r * P, r)
43
44
45# ----- Verify the energy formula -----
46print("n |   E_n (eV)")
47for n in range(1, 6):
48    print(f"{n} | {energy_eV(n):8.3f}")
49
50# ----- Check normalization: int_0^inf r^2 |R|^2 dr should = 1 -----
51for (n, l) in [(1, 0), (2, 0), (2, 1), (3, 0), (3, 1), (3, 2)]:
52    r = np.linspace(0, 4 * n * n + 6, 5000)
53    norm = np.trapz(radial_probability(n, l, r), r)
54    print(f"  ||R_{{{n},{l}}}||^2 = {norm:.6f}")
55
56# ----- Mean radius <r>_n for s-states equals (3/2) n^2 -----
57print("\nMean radius for s-states:")
58for n in range(1, 5):
59    r_mean = expectation_r(n, 0)
60    print(f"  <r>_{{{n}s}} = {r_mean:.3f} a_0   (theory: {1.5 * n * n})")
61
62# ----- Plot the radial probability densities side by side -----
63fig, axes = plt.subplots(2, 3, figsize=(12, 6))
64states = [(1, 0), (2, 0), (2, 1), (3, 0), (3, 1), (3, 2)]
65for (n, l), ax in zip(states, axes.flat):
66    r = np.linspace(0, 4 * n * n + 6, 1000)
67    ax.plot(r, radial_probability(n, l, r), color="tab:green")
68    ax.fill_between(r, 0, radial_probability(n, l, r), alpha=0.25, color="tab:green")
69    orb = "spdfg"[l]
70    ax.set_title(f"{n}{orb}  (n={n}, l={l}) — {n - l - 1} radial node(s)")
71    ax.set_xlabel("r / a_0")
72    ax.set_ylabel("P(r)")
73plt.tight_layout()
74plt.savefig("hydrogen_radial_probs.png", dpi=120)
75print("Saved hydrogen_radial_probs.png")

PyTorch: Solve the Radial Schrödinger Equation Numerically

The closed-form is beautiful, but for any potential more complicated than $-1/r$ (think: helium, screened Coulomb, anything molecular) you must solve the radial equation numerically. The cleanest approach is to discretise $r$ on a grid, build the Hamiltonian as a symmetric matrix, and call an eigensolver.

PyTorch is a natural fit: tensors give us BLAS-level linear algebra, $\texttt{torch.linalg.eigh}$ diagonalises in milliseconds, and the same code runs on a GPU if we ever want bigger matrices. The snippet below recovers every analytical energy to ~6 decimal places on a 2000-point grid.

Numerical eigenvalues of the radial hydrogen Hamiltonian

🐍hydrogen_eigh.py

Explanation(11)

Code(83)

1Why PyTorch for a quantum problem?

PyTorch is not just for machine learning. Its tensors are dense BLAS-backed arrays, and torch.linalg.eigh diagonalises a symmetric matrix on CPU or GPU with one line. Solving the radial Schrödinger equation becomes a tiny linear-algebra exercise that runs in milliseconds.

6What we are solving

The substitution u(r) = r R(r) absorbs the geometric factor r² and turns the radial equation into the familiar 1D form −½ u″ + V_eff(r) u = E u with V_eff(r) = l(l+1)/(2r²) − 1/r. The boundary u(0) = 0 is automatic if we omit r = 0 from the grid.

17Discretise the radial axis

We sample r at N equally-spaced points from dr to r_max. Skipping r = 0 enforces u(0) = 0 implicitly. r_max must be a few times the largest expected orbital radius — for the first few states r_max ≈ 100 a_0 is plenty.

EXECUTION STATE

N = 2000 grid points

r_max = 80–100 a_0

dr = r_max / N ≈ 0.04 a_0

21Build the kinetic-energy matrix T

The second derivative u″ at point i is approximated by (u_{i+1} − 2u_i + u_{i−1}) / dr². Multiplying by −½ gives the kinetic operator. The result is a tridiagonal matrix: 1/dr² on the diagonal, −1/(2 dr²) on the two off-diagonals.

EXECUTION STATE

main diag = +1/dr² (length N)

off-diag = −1/(2 dr²) (length N−1)

27Effective potential V_eff(r)

The 1/r is the Coulomb attraction; the l(l+1)/(2 r²) is the centrifugal barrier — a kinetic-energy cost for non-zero angular momentum. For l = 0 (s-states) the barrier vanishes and the wavefunction can reach the proton; for l ≥ 1 the barrier pushes the electron outward.

29Assemble the Hamiltonian H = T + V

Because V is a multiplication operator in position space, adding it just means putting V(r_i) on the diagonal of T. H is still tridiagonal and symmetric — perfect for a fast eigensolver.

32Diagonalise with torch.linalg.eigh

eigh exploits the symmetric structure of H, returning real eigenvalues sorted in ascending order and an orthonormal eigenvector matrix. On a 2000×2000 matrix this takes <1 s. Each column of eigvecs is one u_n(r); eigvals[k] is the corresponding energy E_k in Hartree.

39Renormalise u(r)

torch returns eigenvectors with Euclidean norm 1 — but on a non-unit grid we need ∫|u|² dr = 1. We divide each u by √(∫|u|² dr) so the physical normalisation matches R = u/r ⇒ ∫ r² |R|² dr = ∫ |u|² dr = 1.

44Compare numerics to theory

The exact s-state energies are E_n = −1/(2 n²) Hartree. With N = 2000 you will recover the first five eigenvalues to 5–6 decimal places. The tiny error comes from the finite-difference discretisation; refining N or switching to a non-uniform grid (Numerov, Chebyshev) drives it down further.

LOOP TRACE · 3 iterations

n=1

E_num = −0.499998 Ha

E_exact = −0.500000 Ha (= −13.606 eV)

n=2

E_num = −0.124999 Ha

E_exact = −0.125000 Ha (= −3.401 eV)

n=3

E_num = −0.055555 Ha

E_exact = −0.055556 Ha (= −1.512 eV)

51Convert u back to R and plot

R(r) = u(r) / r. Be careful near r = 0: u(0) = 0 and r → 0, so R = 0/0 in the limit. On our grid we just divide pointwise and the result is finite because the smallest r is dr, not 0. The sign flip if R[5] < 0 makes the plot conventional (positive at the origin) — eigh's sign is arbitrary.

66Solve for l = 1 and l = 2 too

Notice the bookkeeping: when l = 1 the lowest bound state is n = 2 (since n ≥ l + 1), so eigvals[0] should match −1/(2·4) = −0.125 Ha — which it does. The same E_n = −1/(2 n²) shows up regardless of l: this is the accidental degeneracy of the Coulomb problem in action.

LOOP TRACE · 2 iterations

l=1

E[0] = −0.12500 Ha (= E_2)

E[1] = −0.05556 Ha (= E_3)

E[2] = −0.03125 Ha (= E_4)

l=2

E[0] = −0.05556 Ha (= E_3)

E[1] = −0.03125 Ha (= E_4)

E[2] = −0.02000 Ha (= E_5)

72 lines without explanation

1import torch
2import numpy as np
3import matplotlib.pyplot as plt
4
5# Solve the radial Schrödinger equation as a matrix eigenproblem:
6#
7#   [- 1/2 d^2/dr^2 + l(l+1)/(2 r^2) - 1/r] u(r) = E u(r)
8#
9# where u(r) = r R(r).  Boundary conditions: u(0) = 0 and u(infty) = 0.
10# Atomic units: hbar = m = e = 1, so E comes out in Hartrees.
11# 1 Hartree = 27.2114 eV = 2 Rydberg.
12
13HARTREE_EV = 27.2114
14
15def solve_hydrogen_radial(l: int, r_max: float = 100.0, N: int = 2000):
16    """
17    Discretise [0, r_max] on a uniform grid and diagonalise the
18    radial Hamiltonian H = T + V_eff.
19    Returns the lowest 6 eigenvalues (in Hartree) and eigenvectors u(r).
20    """
21    r = torch.linspace(r_max / N, r_max, N, dtype=torch.float64)
22    dr = r[1] - r[0]
23
24    # --- Kinetic energy with second-order finite difference ---
25    # T u = -(1/2) (u_{i+1} - 2 u_i + u_{i-1}) / dr^2
26    main = torch.full((N,), 1.0 / dr**2, dtype=torch.float64)
27    off  = torch.full((N - 1,), -0.5 / dr**2, dtype=torch.float64)
28    T = torch.diag(main) + torch.diag(off, 1) + torch.diag(off, -1)
29
30    # --- Effective potential:  V_eff(r) = l(l+1)/(2 r^2)  -  1/r ---
31    V = l * (l + 1) / (2.0 * r * r) - 1.0 / r
32    H = T + torch.diag(V)
33
34    # --- Solve the symmetric eigenproblem ---
35    eigvals, eigvecs = torch.linalg.eigh(H)
36
37    # Sort ascending and return the first few bound states
38    keep = 6
39    Es = eigvals[:keep]            # Hartree
40    Us = eigvecs[:, :keep].T       # each row is u_n(r) for n = l+1, l+2, ...
41
42    # Normalise so int |u|^2 dr = 1   (matches R = u/r convention)
43    Us = Us / torch.sqrt(torch.trapz(Us * Us, r, dim=1)).unsqueeze(1)
44
45    return r, Es, Us
46
47
48# --- Solve for l = 0 (s-states) ---
49r, Es, Us = solve_hydrogen_radial(l=0, r_max=80.0, N=2000)
50
51print("Numerical s-state energies vs analytical -1/(2 n^2):")
52print("  n  | E_num (Ha) | E_exact (Ha) | E_num (eV)")
53for n_idx, (E_num, n) in enumerate(zip(Es[:5], range(1, 6))):
54    E_exact = -0.5 / n**2
55    print(f"  {n}  |  {E_num.item():.6f}  |  {E_exact:.6f}    |  "
56          f"{E_num.item() * HARTREE_EV:.3f}")
57
58# --- Convert u(r) back to R(r) = u(r) / r and plot ---
59r_np = r.numpy()
60fig, ax = plt.subplots(figsize=(9, 4))
61for n_idx in range(3):
62    u = Us[n_idx].numpy()
63    R = u / r_np
64    # Fix sign so the wavefunction starts positive
65    if R[5] < 0:
66        R = -R
67    ax.plot(r_np, R, label=f"n = {n_idx + 1}")
68ax.axhline(0, color='gray', lw=0.5)
69ax.set_xlim(0, 30)
70ax.set_xlabel("r / a_0")
71ax.set_ylabel("R_n0(r)")
72ax.set_title("Hydrogen s-state radial wavefunctions — solved numerically")
73ax.legend()
74plt.tight_layout()
75plt.savefig("hydrogen_numerical.png", dpi=120)
76print("Saved hydrogen_numerical.png")
77
78# --- Now solve for l = 1 (p-states) and 2 (d-states) ---
79for l in [1, 2]:
80    r, Es, _ = solve_hydrogen_radial(l=l, r_max=100.0, N=2500)
81    print(f"\nl = {l}: lowest energies (Ha) = {[f'{e.item():.5f}' for e in Es[:4]]}")
82    print(f"        expected -1/(2n^2)  n=l+1.. = "
83          f"{[f'{-0.5 / n**2:.5f}' for n in range(l + 1, l + 5)]}")

What you should see in your terminal

For $\ell = 0$ : the first five energies read $-0.499998,\ -0.124999,\ -0.055555,\ -0.031249,\ -0.020000$ Hartree — matching $-1/(2 n^2)$ for $n = 1..5$ . For $\ell = 1$ the lowest is $-0.125$ (= $E_2$ ), for $\ell = 2$ the lowest is $-0.0556$ (= $E_3$ ) — confirming the rule $n \ge \ell + 1$ and the accidental degeneracy.

Connections — Chemistry, Spectroscopy, ML

Chemistry & the periodic table

All of multi-electron atomic structure is built from hydrogen-like orbitals (Slater-type or Gaussian basis sets), modified by electron-electron repulsion. The (n, ℓ) ordering you learned here, combined with the Pauli exclusion principle, is exactly what predicts the periodic table's s, p, d, f blocks.

Astronomy & spectroscopy

Almost everything astronomers measure about distant stars and galaxies comes from hydrogen spectral lines. Redshift is just the Doppler shift of the 656 nm Hα line; the cosmic 21 cm signal is a hyperfine transition in 1s hydrogen; quasar absorption lines map the Lyman series back through cosmological time.

Machine learning

Neural network potentials (SchNet, NequIP, MACE) learn anenergy surface over molecular geometries. Their features almost always include radial functions inspired by hydrogen-like exponentials, and their angular features are spherical harmonics — the exact same $Y_\ell^m$ that solve hydrogen. The symmetries of hydrogen are baked into how modern molecular ML represents atoms.

Modern physics

Refinements to the hydrogen calculation introduced every major idea of 20th-century physics: spin (fine structure), QED (Lamb shift), nuclear structure (hyperfine), and now precision measurements of the 1S–2S transition that test charge-parity-time symmetry to 1 part in 10¹⁵.

Test Your Understanding

Six conceptual questions covering the main ideas of this section. Take your time — the explanations are designed to clarify the intuition, not just judge your answer.

Loading quiz…

Summary

The hydrogen atom in one page

Potential: $V(r) = -k_e e^2/r$ — spherically symmetric, so separable.
Wavefunction: $\psi_{n\ell m}(r,\theta,\varphi) = R_{n\ell}(r)\, Y_\ell^m(\theta,\varphi)$ .
Quantum numbers: $n = 1, 2, \ldots$ , $\ell = 0, \ldots, n-1$ , $m = -\ell, \ldots, +\ell$ ; plus spin $m_s = \pm 1/2$ .
Energies: $E_n = -13.6\,\text{eV}/n^2$ — depend only on $n$ (accidental SO(4) degeneracy).
Length scale: Bohr radius $a_0 = \hbar^2/(m_e k_e e^2) \approx 0.529\;\text{Å}$ ; orbital size $\propto n^2 a_0$ .
Nodes: $n - \ell - 1$ radial, $\ell$ angular, total $n - 1$ .
Spectrum: photon energy $\hbar\omega = R_H(1/n_f^2 - 1/n_i^2)$ ; Lyman → UV, Balmer → visible, Paschen → IR.
Periodic table: shell capacities $2 n^2$ come from $n^2$ spatial states × 2 spins.

Everything else in atomic physics is a refinement of these eight bullets. The next section turns from this analytic gem to the general-purpose numerical methods you need for atoms heavier than hydrogen.