Hydrogen Orbitals: Spherical Harmonics

The hydrogen Hamiltonian in spherical coordinates

The hydrogen atom is a single electron of mass m bound to a proton by the Coulomb potential V(r) = −e²/(4πε₀r). Because the potential depends only on the distance r, not on direction, the problem has spherical symmetry, and the TISE is best written in spherical coordinates (r, θ, φ):

The angular part of this Laplacian is, up to a factor, the operator −L̂²/ℏ² for the total orbital angular momentum. This is the crucial observation that makes separation of variables possible: the spherical symmetry of V(r) means Ĥ commutes with L̂² and L̂z, so the three operators share a common eigenbasis.

Separation of variables: R(r)·Yₗᵐ(θ,φ)

Because Ĥ, L̂² and L̂z commute, we look for solutions of the product form ψ(r, θ, φ) = R(r) · Y(θ, φ). Substituting into the TISE splits it into two independent equations — one purely angular, one purely radial, coupled only through the separation constant l(l+1):

The term l(l+1)ℏ²/(2mr²) in the radial equation is the centrifugal barrier — the same effective repulsion that appears in classical orbital mechanics, here quantised. It pushes states of higher angular momentum away from the origin, which is exactly why p, d, f orbitals have a node (zero probability) at r = 0 while s orbitals do not.

The angular part: spherical harmonics

The angular equation is solved by the spherical harmonics Yₗᵐ(θ, φ) — the same functions that describe multipole fields in electromagnetism and the modes of a vibrating sphere. They factor further into a φ-dependent phase and a θ-dependent associated Legendre polynomial:

The quantum number l = 0, 1, 2, 3 … is conventionally labelled s, p, d, f (a historical spectroscopic notation), and m ranges over 2l+1 integer values from −l to +l for each l — the origin of the familiar "one s orbital, three p orbitals, five d orbitals" rule.

The radial part: Laguerre polynomials

Solving the radial equation with the Coulomb potential requires the substitution R(r) = e^−r/na₀ · (r/a₀)^l · L(r), where a₀ = 4πε₀ℏ²/(me²) is the Bohr radius. The remaining function L(r) must be a polynomial (for the wave function to stay finite as r → ∞) — this quantisation requirement forces the associated Laguerre polynomials L_n−l−1^2l+1, and, crucially, fixes the energy:

Remarkably, Eₙ depends only on the principal quantum number n, not on l or m — a special "accidental" degeneracy of the pure 1/r Coulomb potential (it disappears for any other central potential, and is lifted in real multi-electron atoms by electron-electron screening).

Quantum numbers n, l, m and their rules

For a given n there are n² distinct orbitals (n choices of l times 2l+1 choices of m, summed = n²) — before accounting for the two spin states of the electron, which double this to 2n² and correctly reproduce the lengths of the periods of the periodic table.

Orbital shapes: s, p, d

• s orbitals (l=0): spherically symmetric, no angular dependence, non-zero probability density at the nucleus. The number of radial nodes is n−1.
• p orbitals (l=1): two lobes of opposite sign separated by a nodal plane through the nucleus; three orientations (px, py, pz) point along the Cartesian axes.
• d orbitals (l=2): four-lobed "cloverleaf" shapes (dxy, dxz, dyz, dx²−y²) or a two-lobe-plus-ring shape (dz²), with two angular nodal surfaces.

Rendering electron probability clouds

To visualise |ψ(r,θ,φ)|² in WebGL, sample points via rejection sampling in the radial and angular probability densities, then plot them as a point cloud — far cheaper than volumetric ray marching for the sparse, lobed hydrogen densities:

// Rejection-sample a hydrogen orbital point cloud (n, l, m)
function sampleOrbital(n, l, m, N, a0 = 1) {
  const pts = [];
  const rMax = a0 * n * n * 6;   // generous cutoff radius
  while (pts.length < N) {
    const r = Math.random() * rMax;
    const theta = Math.acos(2 * Math.random() - 1);
    const phi = 2 * Math.PI * Math.random();
    const psi2 = radialR(n, l, r, a0) ** 2 * sphericalHarmonicReal(l, m, theta, phi) ** 2;
    if (Math.random() < psi2 / PSI2_MAX) {
      const x = r * Math.sin(theta) * Math.cos(phi);
      const y = r * Math.sin(theta) * Math.sin(phi);
      const z = r * Math.cos(theta);
      pts.push(x, y, z);
    }
  }
  return new Float32Array(pts);
}