Part I, Topic 1 | Lectures 1–4

The Hydrogen Atom

The only exactly solvable atomic system — cornerstone of quantum chemistry and atomic physics

1.1 Introduction

The hydrogen atom occupies a unique place in the history of physics and chemistry. As the simplest atom — a single electron bound to a single proton — it is the only atomic system for which the Schrödinger equation can be solved exactly, without any approximations. Every aspect of multi-electron quantum chemistry ultimately traces back to the hydrogen atom solution.

The Hamiltonian for the hydrogen atom consists of the kinetic energy of the electron plus the Coulomb attraction between the electron (charge $-e$) and the proton (charge $+e$):

$$\hat{H} = -\frac{\hbar^2}{2m_e}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 r}$$

Here $r$ is the distance between the electron and the proton, $m_e$ is the electron mass, and $e$ is the elementary charge. In Gaussian units (commonly used in atomic physics), the potential simplifies to $V(r) = -e^2/r$.

Because the potential depends only on the radial distance $r$, and not on the angular coordinates $\theta$ and $\varphi$, the problem possesses full spherical symmetry. This is the key insight that makes separation of variables possible: we can decompose the three-dimensional partial differential equation into three ordinary differential equations, each involving a single coordinate.

The exact solution yields three quantum numbers — $n$, $l$, and $m_l$ — that completely characterize the spatial part of the wavefunction. These quantum numbers determine the energy, the total angular momentum magnitude, and the z-component of angular momentum, respectively. The resulting wavefunctions, the hydrogen atomic orbitals, are products of radial functions and spherical harmonics.

Why Hydrogen Matters

It is the only atom with an exact analytical solution to the Schrödinger equation
The wavefunctions serve as the basis for understanding all multi-electron atoms
Hydrogen spectral lines provided the first quantitative evidence for energy quantization
Hydrogen-like ions (He$^+$, Li$^{2+}$, etc.) are solved by simple scaling
The solution introduces concepts (quantum numbers, angular momentum, degeneracy) used throughout chemistry and physics

1.2 Separation of Variables in Spherical Coordinates

We begin with the full three-dimensional time-independent Schrödinger equation in spherical coordinates $(r, \theta, \varphi)$. The Laplacian in spherical coordinates is:

$$\nabla^2 = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta\frac{\partial}{\partial \theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2}{\partial \varphi^2}$$

The Schrödinger equation $\hat{H}\psi = E\psi$ becomes:

$$-\frac{\hbar^2}{2m_e}\left[\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial \psi}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta\frac{\partial \psi}{\partial \theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 \psi}{\partial \varphi^2}\right] - \frac{e^2}{4\pi\epsilon_0 r}\psi = E\psi$$

We assume the wavefunction separates as a product of functions of individual coordinates:

$$\psi(r, \theta, \varphi) = R(r) \cdot \Theta(\theta) \cdot \Phi(\varphi)$$

Substituting this product into the Schrödinger equation, multiplying through by$-2m_e r^2 \sin^2\theta / (\hbar^2 R\Theta\Phi)$, and rearranging, we isolate terms that depend on different coordinates. The standard separation procedure yields three ordinary differential equations.

Step 1: The Azimuthal Equation ($\Phi$)

The $\varphi$-dependent part separates cleanly because it appears only in the $\partial^2/\partial\varphi^2$ term:

$$\frac{d^2\Phi}{d\varphi^2} = -m_l^2 \Phi$$

where $m_l^2$ is the separation constant. The solution is:

$$\Phi(\varphi) = \frac{1}{\sqrt{2\pi}} e^{i m_l \varphi}$$

The single-valuedness condition $\Phi(\varphi + 2\pi) = \Phi(\varphi)$ requires that $m_l$ be an integer: $m_l = 0, \pm 1, \pm 2, \ldots$. This is the magnetic quantum number.

Step 2: The Polar Equation ($\Theta$)

After separating out the azimuthal part, the $\theta$-dependent equation becomes (with the substitution $x = \cos\theta$):

$$\frac{d}{dx}\left[(1-x^2)\frac{d\Theta}{dx}\right] + \left[l(l+1) - \frac{m_l^2}{1-x^2}\right]\Theta = 0$$

This is the associated Legendre equation. For the solutions to be well-behaved (finite and single-valued) over the range $0 \le \theta \le \pi$, the separation constant must take the form $l(l+1)$ where:

$$l = 0, 1, 2, \ldots \quad \text{and} \quad |m_l| \le l$$

The quantum number $l$ is the angular momentum quantum number(also called the azimuthal or orbital quantum number). The solutions are the associated Legendre polynomials $P_l^{m_l}(\cos\theta)$:

$$\Theta_{l,m_l}(\theta) = N_{l,m_l} P_l^{|m_l|}(\cos\theta)$$

where $N_{l,m_l}$ is a normalization constant. The constraint $|m_l| \le l$means that for a given $l$, there are $2l + 1$ allowed values of $m_l$.

Step 3: The Radial Equation ($R$)

The remaining radial equation, after separating the angular parts, is:

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) + \left[\frac{2m_e}{\hbar^2}\left(E + \frac{e^2}{4\pi\epsilon_0 r}\right) - \frac{l(l+1)}{r^2}\right]R = 0$$

This equation contains the Coulomb potential and the centrifugal barrier term $l(l+1)/r^2$. It is convenient to define $u(r) = rR(r)$, which transforms the equation into a form resembling a one-dimensional Schrödinger equation with an effective potential:

$$-\frac{\hbar^2}{2m_e}\frac{d^2 u}{dr^2} + \left[-\frac{e^2}{4\pi\epsilon_0 r} + \frac{\hbar^2 l(l+1)}{2m_e r^2}\right]u = Eu$$

The solution of this equation, demanding that the wavefunction be normalizable (bound states with $E < 0$), introduces the principal quantum number$n = 1, 2, 3, \ldots$ with the constraint $l \le n - 1$. The full derivation is presented in the next section.

1.3 The Radial Equation and Energy Eigenvalues

To solve the radial equation, we introduce dimensionless variables. Define the Bohr radius:

$$a_0 = \frac{4\pi\epsilon_0\hbar^2}{m_e e^2} = \frac{\hbar^2}{m_e e^2} \approx 0.529\,\text{\AA}$$

(The second expression uses Gaussian units.) The Bohr radius sets the natural length scale for the hydrogen atom. Now introduce the dimensionless radial variable:

$$\rho = \frac{2r}{na_0}$$

where $n$ will turn out to be the principal quantum number. With this substitution, and defining $\kappa^2 = -2m_e E/\hbar^2$ (since $E < 0$ for bound states,$\kappa$ is real), the radial equation transforms to:

$$\frac{d^2 u}{d\rho^2} + \left[-\frac{1}{4} + \frac{n}{\rho} - \frac{l(l+1)}{\rho^2}\right]u = 0$$

Asymptotic Analysis

We analyze the behavior of $u(\rho)$ in the two limits:

As $\rho \to \infty$: the equation becomes $d^2u/d\rho^2 \approx u/4$, giving $u \sim e^{-\rho/2}$ (rejecting the growing exponential)
As $\rho \to 0$: the centrifugal term dominates, giving $u \sim \rho^{l+1}$ (rejecting the singular solution $\rho^{-l}$)

This motivates the ansatz:

$$u(\rho) = \rho^{l+1} e^{-\rho/2} v(\rho)$$

Series Solution and Quantization

Substituting this ansatz into the radial equation yields a differential equation for $v(\rho)$:

$$\rho \frac{d^2v}{d\rho^2} + (2l + 2 - \rho)\frac{dv}{d\rho} + (n - l - 1)v = 0$$

This is the associated Laguerre equation. We expand $v(\rho)$ as a power series $v(\rho) = \sum_{j=0}^{\infty} c_j \rho^j$ and substitute to obtain the recurrence relation:

$$c_{j+1} = \frac{j + l + 1 - n}{(j+1)(j + 2l + 2)} c_j$$

If this series does not terminate, the large-$j$ behavior of the coefficients gives $c_{j+1}/c_j \to 1/j$, which means $v(\rho) \sim e^{\rho}$, and so $u(\rho) \sim \rho^{l+1}e^{\rho/2}$, which diverges as $\rho \to \infty$. For a normalizable wavefunction, the series must terminate. This requires:

$$n - l - 1 = j_{\max} \ge 0 \quad \Rightarrow \quad n = l + 1 + j_{\max}$$

Since $j_{\max}$ is a non-negative integer, $n$ must be a positive integer with $n \ge l + 1$, i.e., $l = 0, 1, \ldots, n-1$. The terminated series gives the associated Laguerre polynomials$L_{n-l-1}^{2l+1}(\rho)$.

The Energy Eigenvalues

The quantization condition, combined with the relation between $\kappa$ and $E$, gives the celebrated hydrogen energy levels:

$$E_n = -\frac{m_e e^4}{2\hbar^2(4\pi\epsilon_0)^2} \cdot \frac{1}{n^2} = -\frac{13.6\,\text{eV}}{n^2}, \quad n = 1, 2, 3, \ldots$$

Key features of this result:

The energy depends only on $n$, not on $l$ or $m_l$. This gives rise to a characteristic degeneracy: for each $n$, the number of degenerate states is $\sum_{l=0}^{n-1}(2l+1) = n^2$.
The ground state energy is $E_1 = -13.6$ eV, matching the ionization energy of hydrogen.
Energy levels converge as $n \to \infty$, with the ionization threshold at $E = 0$.
This $1/n^2$ dependence exactly reproduces the empirical Rydberg formula discovered in the 19th century.

The Complete Radial Wavefunction

Combining the asymptotic factors and the Laguerre polynomial, the normalized radial wavefunction is:

$$R_{nl}(r) = -\sqrt{\left(\frac{2}{na_0}\right)^3 \frac{(n-l-1)!}{2n[(n+l)!]^3}} \; e^{-r/(na_0)} \left(\frac{2r}{na_0}\right)^l L_{n-l-1}^{2l+1}\!\left(\frac{2r}{na_0}\right)$$

The number of radial nodes (zeros of $R_{nl}$ for $r > 0$) is $n - l - 1$. The most probable radius — the maximum of $r^2|R_{nl}|^2$ — for the ground state ($n = 1, l = 0$) occurs at $r = a_0$, the Bohr radius.

1.4 Angular Momentum and Spherical Harmonics

The angular part of the hydrogen wavefunction is universally important in quantum mechanics. The product $\Theta(\theta)\Phi(\varphi)$ forms the spherical harmonics, which are eigenfunctions of both the total angular momentum squared operator $\hat{L}^2$and its z-component $\hat{L}_z$.

Definition of Spherical Harmonics

The spherical harmonics are defined as:

$$Y_l^{m_l}(\theta, \varphi) = (-1)^{m_l} \sqrt{\frac{2l+1}{4\pi}\frac{(l-|m_l|)!}{(l+|m_l|)!}} \; P_l^{|m_l|}(\cos\theta) \; e^{im_l\varphi}$$

where $P_l^{|m_l|}$ is the associated Legendre polynomial and the prefactor$(-1)^{m_l}$ is the Condon-Shortley phase convention. These functions form a complete orthonormal set on the unit sphere:

$$\int_0^{2\pi}\int_0^{\pi} Y_{l'}^{m'*}(\theta,\varphi) \, Y_l^m(\theta,\varphi) \sin\theta \, d\theta \, d\varphi = \delta_{ll'}\delta_{mm'}$$

Eigenvalue Equations

The angular momentum operators in spherical coordinates are:

$$\hat{L}_z = -i\hbar\frac{\partial}{\partial \varphi}$$

$$\hat{L}^2 = -\hbar^2\left[\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right) + \frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\varphi^2}\right]$$

Acting on the spherical harmonics:

$$\hat{L}^2 Y_l^{m_l}(\theta,\varphi) = l(l+1)\hbar^2 \, Y_l^{m_l}(\theta,\varphi)$$

$$\hat{L}_z Y_l^{m_l}(\theta,\varphi) = m_l\hbar \, Y_l^{m_l}(\theta,\varphi)$$

These eigenvalue equations tell us that:

The magnitude of the orbital angular momentum is $|\mathbf{L}| = \hbar\sqrt{l(l+1)}$
The z-component of angular momentum is $L_z = m_l\hbar$
Since $\sqrt{l(l+1)} > l$ for $l \ge 1$, the angular momentum vector can never be perfectly aligned with the z-axis — a consequence of the uncertainty principle
The x and y components of $\mathbf{L}$ are not simultaneously determined (they do not commute with $\hat{L}_z$)

Explicit Forms of Low-Order Spherical Harmonics

$Y_0^0 = \dfrac{1}{\sqrt{4\pi}}$ — spherically symmetric (s orbital)

$Y_1^0 = \sqrt{\dfrac{3}{4\pi}}\cos\theta$ — p$_z$ orbital

$Y_1^{\pm 1} = \mp\sqrt{\dfrac{3}{8\pi}}\sin\theta \, e^{\pm i\varphi}$ — p$_{\pm 1}$ orbitals

$Y_2^0 = \sqrt{\dfrac{5}{16\pi}}(3\cos^2\theta - 1)$ — d$_{z^2}$ orbital

The Complete Hydrogen Wavefunction

The full hydrogen atom wavefunction is the product of the radial and angular parts:

$$\psi_{nlm_l}(r,\theta,\varphi) = R_{nl}(r) \cdot Y_l^{m_l}(\theta,\varphi)$$

1.5 Selection Rules for Hydrogen

Not all transitions between hydrogen energy levels are allowed. The probability of an electric dipole transition between states $|n, l, m_l\rangle$ and $|n', l', m_l'\rangle$ is proportional to the square of the transition dipole matrix element:

$$\boldsymbol{\mu}_{fi} = \langle n', l', m_l' | \, e\hat{\mathbf{r}} \, | n, l, m_l \rangle = e\int \psi_{n'l'm_l'}^{*} \, \mathbf{r} \, \psi_{nlm_l} \, d^3r$$

A transition is "allowed" if this integral is nonzero. Since $\mathbf{r}$ can be expressed in terms of spherical harmonics of order 1 (i.e., $Y_1^m$), the angular integral factorizes as:

$$\int Y_{l'}^{m_l'*}(\theta,\varphi) \, Y_1^q(\theta,\varphi) \, Y_l^{m_l}(\theta,\varphi) \, \sin\theta \, d\theta \, d\varphi$$

where $q = 0, \pm 1$ corresponds to the z, x, and y components of $\mathbf{r}$. This integral of three spherical harmonics can be evaluated using the properties of the Clebsch-Gordan coefficients (or equivalently, the Wigner 3j symbols). The integral vanishes unless the following conditions are satisfied:

Electric Dipole Selection Rules

$$\Delta l = l' - l = \pm 1$$

$$\Delta m_l = m_l' - m_l = 0, \pm 1$$

There is no restriction on $\Delta n$ (transitions between any principal quantum numbers are allowed, provided the angular momentum rules are satisfied).

Physical Interpretation

The selection rule $\Delta l = \pm 1$ has a deep physical origin: the emitted or absorbed photon carries one unit of angular momentum ($s = 1$ for a photon). Conservation of angular momentum therefore requires the atom's orbital angular momentum to change by exactly one unit.

This means, for example:

$2p \to 1s$: allowed ($\Delta l = -1$) — this is the Lyman-alpha transition
$2s \to 1s$: forbidden ($\Delta l = 0$) — the 2s state is metastable
$3d \to 2p$: allowed ($\Delta l = -1$)
$3d \to 1s$: forbidden ($\Delta l = -2$)
$3s \to 2p$: allowed ($\Delta l = +1$)

The metastability of the $2s$ state is physically significant: it can only decay via two-photon emission (a much slower process) or by magnetic dipole / electric quadrupole transitions, which have much smaller rates than electric dipole transitions.

1.6 Applications

Atomic Spectroscopy: The Hydrogen Spectral Series

The Rydberg formula gives the wavelength of light emitted when an electron transitions from level $n_2$ to level $n_1$:

$$\frac{1}{\lambda} = R_\infty\left(\frac{1}{n_1^2} - \frac{1}{n_2^2}\right), \quad n_2 > n_1$$

where $R_\infty = 1.097 \times 10^7$ m$^{-1}$ is the Rydberg constant. The named series are:

Lyman series ($n_1 = 1$): ultraviolet, 91.2–121.6 nm
Balmer series ($n_1 = 2$): visible/near-UV, 364.6–656.3 nm
Paschen series ($n_1 = 3$): near-infrared, 820.4–1875.1 nm
Brackett series ($n_1 = 4$): infrared
Pfund series ($n_1 = 5$): far-infrared

Hydrogen in Astrophysics

Hydrogen is the most abundant element in the universe, comprising about 75% of all baryonic matter by mass. Its spectral lines are among the most important tools in astrophysics:

The 21 cm line (hyperfine transition of the ground state) is used to map the distribution of neutral hydrogen in galaxies and the intergalactic medium
The Lyman-alpha forest in quasar spectra reveals the distribution of intergalactic hydrogen and is a probe of large-scale cosmic structure
The Balmer series (particularly H-alpha at 656.3 nm) is used to identify star-forming regions and measure the temperatures of stellar atmospheres
Redshifted hydrogen lines are the primary tool for measuring the expansion rate of the universe

Hydrogen-like Ions

Any one-electron ion with nuclear charge $Z$ (such as He$^+$, Li$^{2+}$, Be$^{3+}$, etc.) has exactly the same mathematical structure as hydrogen, with the replacement $e^2 \to Ze^2$. The energy levels become:

$$E_n = -\frac{Z^2 \times 13.6\,\text{eV}}{n^2}$$

and the Bohr radius scales as $a_0/Z$, meaning the orbitals contract for higher nuclear charges. The wavefunctions are obtained by the simple substitution $a_0 \to a_0/Z$in all formulas.

Hydrogen-like ions are important in plasma physics, X-ray astronomy, and as test systems for high-precision QED calculations. For example, measurements on hydrogen-like uranium (U$^{91+}$) test quantum electrodynamics in extreme electric fields.

Quantum Defect Theory for Alkali Metals

The alkali metals (Li, Na, K, Rb, Cs) each have a single valence electron, making them approximately hydrogen-like. Their energy levels are well-described by the quantum defect formula:

$$E_{nl} = -\frac{13.6\,\text{eV}}{(n - \delta_l)^2}$$

where $\delta_l$ is the quantum defect, which depends on $l$ but is approximately independent of $n$. The quantum defect is largest for s orbitals ($l = 0$), which penetrate most deeply into the core electron cloud, and decreases rapidly with increasing $l$. For sodium, $\delta_0 \approx 1.37$,$\delta_1 \approx 0.88$, $\delta_2 \approx 0.01$. This explains why the hydrogen degeneracy ($E$ independent of $l$) is broken in multi-electron atoms.

1.7 Historical Context

The Bohr Model (1913)

Niels Bohr proposed the first successful model of the hydrogen atom by postulating that the electron moves in circular orbits around the nucleus, with angular momentum quantized in integer multiples of $\hbar$:

$$m_e v r = n\hbar, \quad n = 1, 2, 3, \ldots$$

Combined with the classical force balance (Coulomb force = centripetal force), this yields exactly the same energy levels $E_n = -13.6/n^2$ eV as the full quantum mechanical treatment. The Bohr model correctly predicted the hydrogen spectrum and introduced the concept of quantized energy levels, but it could not explain the fine structure, the Zeeman effect, or the spectra of multi-electron atoms. It also failed to provide wavefunctions or any notion of orbital angular momentum beyond $l = n - 1$.

Schrödinger's Solution (1926)

Erwin Schrödinger published his wave equation in January 1926 and almost immediately applied it to the hydrogen atom. His solution confirmed the Bohr energy levels but also yielded:

The full three-dimensional wavefunctions $\psi_{nlm}$
The correct angular momentum quantum numbers and their relationships
Probability distributions that replaced the classical notion of electron orbits
The $l = 0$ states (s orbitals), which have no classical analogue in the Bohr model

Schrödinger's solution is the one presented in this chapter and remains the foundation of all atomic structure calculations.

Dirac Equation and Fine Structure

Paul Dirac's relativistic wave equation (1928) provided a deeper treatment of the hydrogen atom. The Dirac solution naturally incorporates:

Electron spin: the spin quantum number $s = 1/2$ emerges naturally, without ad hoc postulates
Spin-orbit coupling: the interaction between the electron's spin and orbital angular momentum splits energy levels
Fine structure: the energy depends on the total angular momentum quantum number $j$, breaking the $l$-degeneracy

The Dirac fine structure formula is:

$$E_{nj} = -\frac{13.6\,\text{eV}}{n^2}\left[1 + \frac{\alpha^2}{n^2}\left(\frac{n}{j + 1/2} - \frac{3}{4}\right)\right]$$

where $\alpha \approx 1/137$ is the fine-structure constant.

The Lamb Shift (1947)

Willis Lamb and Robert Retherford discovered experimentally that the $2s_{1/2}$and $2p_{1/2}$ states of hydrogen, which are degenerate according to the Dirac equation (both have $j = 1/2$), are in fact split by about 1057 MHz (approximately$4.4 \times 10^{-6}$ eV).

This Lamb shift was the first direct evidence of quantum electrodynamic (QED) effects. It arises from the electron's interaction with the quantum vacuum — specifically, the emission and reabsorption of virtual photons, which causes the electron to undergo rapid fluctuations in position ("Zitterbewegung"). This smears out the electron's charge distribution, slightly modifying the Coulomb potential felt by s-orbital electrons (which have nonzero probability density at the nucleus).

Hans Bethe's 1947 calculation of the Lamb shift, in excellent agreement with experiment, was a triumph for quantum electrodynamics and helped establish QED as the most precisely tested theory in all of physics.

1.8 Python Simulation: Hydrogen Radial Wavefunctions

The following simulation computes and plots the hydrogen radial wavefunctions $R_{nl}(r)$and radial probability densities $r^2|R_{nl}(r)|^2$ for the first several orbitals (1s through 3d). The associated Laguerre polynomials are computed using a recurrence relation implemented in pure numpy (no scipy dependency). The normalization is verified by numerical integration using $\texttt{np.trapezoid}$.

Hydrogen Radial Wavefunctions and Probability Densities

Python

hydrogen_radial.py124 lines

import numpy as np
import matplotlib.pyplot as plt

# === Hydrogen Radial Wavefunctions (numpy only, no scipy) ===

def factorial(n):
    """Compute factorial of n."""
    result = 1
    for i in range(2, n + 1):
        result *= i
    return result

def associated_laguerre(p, q, x):
    """Compute associated Laguerre polynomial L_p^q(x) using the recurrence relation.
    L_0^q(x) = 1
    L_1^q(x) = 1 + q - x
    (p+1) L_{p+1}^q(x) = (2p + 1 + q - x) L_p^q(x) - (p + q) L_{p-1}^q(x)
    """
    if p == 0:
        return np.ones_like(x, dtype=float)
    elif p == 1:
        return 1.0 + q - x
    else:
        L_prev2 = np.ones_like(x, dtype=float)
        L_prev1 = 1.0 + q - x
        for k in range(1, p):
            L_curr = ((2*k + 1 + q - x) * L_prev1 - (k + q) * L_prev2) / (k + 1)
            L_prev2 = L_prev1
            L_prev1 = L_curr
        return L_prev1

def R_nl(n, l, r):
    """Radial wavefunction R_nl(r) for hydrogen atom.
    r is in units of a0 (Bohr radius)."""
    rho = 2.0 * r / n
    norm = np.sqrt((2.0 / n)**3 * factorial(n - l - 1) / (2.0 * n * factorial(n + l)))
    L = associated_laguerre(n - l - 1, 2*l + 1, rho)
    return norm * np.exp(-rho / 2) * rho**l * L

# Set up the grid
r = np.linspace(0, 40, 1000)

# Orbitals to plot: (n, l, label, color)
orbitals = [
    (1, 0, '1s', '#34d399'),
    (2, 0, '2s', '#60a5fa'),
    (2, 1, '2p', '#f472b6'),
    (3, 0, '3s', '#fbbf24'),
    (3, 1, '3p', '#a78bfa'),
    (3, 2, '3d', '#fb923c'),
]

fig, axes = plt.subplots(2, 3, figsize=(16, 10))
fig.patch.set_facecolor('#0f172a')
fig.suptitle('Hydrogen Radial Wavefunctions and Probability Densities',
             color='#34d399', fontsize=16, fontweight='bold', y=0.98)

for idx, (n, l, label, color) in enumerate(orbitals):
    row = idx // 3
    col = idx % 3
    ax = axes[row, col]
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    ax.xaxis.label.set_color('#cbd5e1')
    ax.yaxis.label.set_color('#cbd5e1')
    ax.title.set_color(color)
    for sp in ax.spines.values():
        sp.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

R = R_nl(n, l, r)
    P = r**2 * R**2

# Plot R_nl(r) as thin dashed line
    ax.plot(r, R, color=color, lw=1.2, ls='--', alpha=0.5, label=f'R(r)')

# Plot radial probability density r^2|R|^2
    ax.fill_between(r, P, alpha=0.25, color=color)
    ax.plot(r, P, color=color, lw=2.5, label=r'$r^2|R|^2$')

ax.set_title(f'{label} orbital (n={n}, l={l})', fontsize=13, fontweight='bold')
    ax.set_xlabel('r / a\u2080', fontsize=10)
    ax.set_ylabel('Probability density', fontsize=10)

# Mark the most probable radius
    if np.max(P) > 0:
        i_max = np.argmax(P)
        r_max = r[i_max]
        ax.axvline(r_max, color=color, ls=':', alpha=0.6, lw=1)
        ax.text(r_max + 0.5, P[i_max] * 0.85, f'r_mp={r_max:.1f}a\u2080',
                color=color, fontsize=8)

# Number of radial nodes
    nodes = n - l - 1
    ax.text(0.95, 0.95, f'Radial nodes: {nodes}', transform=ax.transAxes,
            color='#94a3b8', fontsize=9, va='top', ha='right',
            bbox=dict(boxstyle='round,pad=0.3', facecolor='#0f172a', edgecolor='#334155'))

ax.legend(loc='upper right', fontsize=8, facecolor='#1e293b',
              edgecolor='#334155', labelcolor='#cbd5e1',
              bbox_to_anchor=(0.95, 0.82))

plt.tight_layout(pad=2, rect=[0, 0, 1, 0.95])
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

# Verify normalization using numpy trapezoid
print("=== Hydrogen Radial Wavefunctions ===")
print()
for n, l, label, _ in orbitals:
    R = R_nl(n, l, r)
    integrand = r**2 * R**2
    norm = np.trapezoid(integrand, r)
    nodes = n - l - 1
    i_max = np.argmax(r**2 * R**2)
    r_max = r[i_max]
    print(f"  {label}: norm = {norm:.4f}, radial nodes = {nodes}, r_mp = {r_max:.1f} a0")

print()
print("For hydrogen: E_n = -13.6 eV / n^2")
for n in range(1, 5):
    print(f"  n={n}: E = {-13.6/n**2:.3f} eV")
print()
print("Bohr radius a0 = 0.529 Angstroms = 52.9 pm")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

1.9 Fortran Simulation: Hydrogen Energy Levels and Spectral Series

This Fortran program computes the hydrogen energy levels for $n = 1$ through $n = 7$, then generates formatted tables of transition wavelengths for the Lyman, Balmer, and Paschen spectral series. The Rydberg energy ($13.606$ eV) and fundamental physical constants are used to calculate wavelengths in nanometers, photon energies in electron volts, and wavenumbers in cm$^{-1}$.

Hydrogen Energy Levels and Spectral Series

Fortran

hydrogen_spectrum.f90109 lines

program hydrogen_spectrum
  implicit none

! === Hydrogen Energy Levels and Spectral Series ===
  ! Compute energy levels, transition wavelengths, and frequencies
  ! for the Lyman, Balmer, and Paschen series.

integer :: n_upper, n_lower, i
  double precision :: E_n_upper, E_n_lower, delta_E
  double precision :: wavelength_nm, frequency_hz, wavenumber
  double precision :: R_inf, h_eV, c_nm, eV_to_J

! Physical constants
  double precision, parameter :: Rydberg_eV = 13.605693d0  ! eV
  double precision, parameter :: h_planck = 6.62607015d-34  ! J*s
  double precision, parameter :: c_light = 2.99792458d8     ! m/s
  double precision, parameter :: eV_J = 1.602176634d-19     ! J/eV
  double precision, parameter :: R_Rydberg = 1.0973731568d7 ! m^-1

write(*,'(A)') '=================================================================='
  write(*,'(A)') '         HYDROGEN ATOM: ENERGY LEVELS AND SPECTRAL SERIES'
  write(*,'(A)') '=================================================================='
  write(*,*)

! Print energy levels
  write(*,'(A)') '--- Energy Levels ---'
  write(*,'(A6, A15, A15)') 'n', 'E_n (eV)', 'E_n (cm^-1)'
  write(*,'(A)') '--------------------------------------'
  do i = 1, 7
    E_n_upper = -Rydberg_eV / dble(i*i)
    wavenumber = -R_Rydberg / dble(i*i) / 100.0d0  ! convert to cm^-1
    write(*,'(I6, F15.4, F15.1)') i, E_n_upper, wavenumber
  end do
  write(*,*)

! === Lyman Series (n_lower = 1, UV) ===
  write(*,'(A)') '=== LYMAN SERIES (n_lower = 1, Ultraviolet) ==='
  write(*,'(A8, A8, A16, A16, A16)') 'n_upper', 'n_low', &
    'lambda (nm)', 'E_photon (eV)', 'wavenumber(cm-1)'
  write(*,'(A)') '----------------------------------------------------------------'

n_lower = 1
  do n_upper = 2, 7
    E_n_upper = -Rydberg_eV / dble(n_upper * n_upper)
    E_n_lower = -Rydberg_eV / dble(n_lower * n_lower)
    delta_E = E_n_upper - E_n_lower
    wavelength_nm = (h_planck * c_light) / (delta_E * eV_J) * 1.0d9
    wavenumber = 1.0d7 / wavelength_nm  ! cm^-1
    write(*,'(I8, I8, F16.2, F16.4, F16.1)') n_upper, n_lower, &
      wavelength_nm, delta_E, wavenumber
  end do
  write(*,'(A,F10.2,A)') '  Series limit: ', &
    (h_planck * c_light) / (Rydberg_eV * eV_J) * 1.0d9, ' nm'
  write(*,*)

! === Balmer Series (n_lower = 2, Visible) ===
  write(*,'(A)') '=== BALMER SERIES (n_lower = 2, Visible/Near-UV) ==='
  write(*,'(A8, A8, A16, A16, A16)') 'n_upper', 'n_low', &
    'lambda (nm)', 'E_photon (eV)', 'wavenumber(cm-1)'
  write(*,'(A)') '----------------------------------------------------------------'

n_lower = 2
  do n_upper = 3, 8
    E_n_upper = -Rydberg_eV / dble(n_upper * n_upper)
    E_n_lower = -Rydberg_eV / dble(n_lower * n_lower)
    delta_E = E_n_upper - E_n_lower
    wavelength_nm = (h_planck * c_light) / (delta_E * eV_J) * 1.0d9
    wavenumber = 1.0d7 / wavelength_nm
    write(*,'(I8, I8, F16.2, F16.4, F16.1)') n_upper, n_lower, &
      wavelength_nm, delta_E, wavenumber
  end do
  write(*,'(A,F10.2,A)') '  Series limit: ', &
    (h_planck * c_light) / (Rydberg_eV * 0.25d0 * eV_J) * 1.0d9, ' nm'
  write(*,*)

! === Paschen Series (n_lower = 3, Infrared) ===
  write(*,'(A)') '=== PASCHEN SERIES (n_lower = 3, Infrared) ==='
  write(*,'(A8, A8, A16, A16, A16)') 'n_upper', 'n_low', &
    'lambda (nm)', 'E_photon (eV)', 'wavenumber(cm-1)'
  write(*,'(A)') '----------------------------------------------------------------'

n_lower = 3
  do n_upper = 4, 9
    E_n_upper = -Rydberg_eV / dble(n_upper * n_upper)
    E_n_lower = -Rydberg_eV / dble(n_lower * n_lower)
    delta_E = E_n_upper - E_n_lower
    wavelength_nm = (h_planck * c_light) / (delta_E * eV_J) * 1.0d9
    wavenumber = 1.0d7 / wavelength_nm
    write(*,'(I8, I8, F16.2, F16.4, F16.1)') n_upper, n_lower, &
      wavelength_nm, delta_E, wavenumber
  end do
  write(*,'(A,F10.2,A)') '  Series limit: ', &
    (h_planck * c_light) / (Rydberg_eV * (1.0d0/9.0d0) * eV_J) * 1.0d9, ' nm'
  write(*,*)

! === Summary ===
  write(*,'(A)') '=== SUMMARY ==='
  write(*,'(A)') 'Rydberg formula: 1/lambda = R_inf * (1/n1^2 - 1/n2^2)'
  write(*,'(A,ES15.7,A)') 'R_infinity = ', R_Rydberg, ' m^-1'
  write(*,'(A,F10.6,A)') 'Rydberg energy = ', Rydberg_eV, ' eV'
  write(*,'(A)') ''
  write(*,'(A)') 'Lyman series: UV region (91.2 - 121.6 nm)'
  write(*,'(A)') 'Balmer series: Visible/near-UV (364.6 - 656.3 nm)'
  write(*,'(A)') 'Paschen series: Near-IR (820.4 - 1875.1 nm)'
  write(*,'(A)') ''
  write(*,'(A)') 'The hydrogen spectrum provided key evidence for'
  write(*,'(A)') 'quantized energy levels in atoms.'

end program hydrogen_spectrum

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Summary of Key Results

Quantum Numbers

$n = 1, 2, 3, \ldots$ (principal); $l = 0, 1, \ldots, n-1$ (angular momentum);$m_l = -l, \ldots, +l$ (magnetic). Total degeneracy of level $n$: $n^2$(or $2n^2$ including spin).

Energy Levels

$E_n = -13.6\,\text{eV}/n^2$. Bohr radius: $a_0 = \hbar^2/(m_e e^2) = 0.529$ \u00C5. Most probable radius for 1s: $r_{mp} = a_0$.

Angular Momentum

$|\mathbf{L}| = \hbar\sqrt{l(l+1)}$, $L_z = m_l\hbar$. Spherical harmonics$Y_l^{m_l}(\theta,\varphi)$ are simultaneous eigenfunctions of $\hat{L}^2$and $\hat{L}_z$.

Selection Rules

Electric dipole transitions require $\Delta l = \pm 1$ and $\Delta m_l = 0, \pm 1$. No restriction on $\Delta n$.

Wavefunction

$\psi_{nlm_l}(r,\theta,\varphi) = R_{nl}(r) \cdot Y_l^{m_l}(\theta,\varphi)$ where $R_{nl}$involves associated Laguerre polynomials and $Y_l^{m_l}$ are spherical harmonics. Radial nodes: $n - l - 1$. Angular nodes: $l$. Total nodes: $n - 1$.

← Part I Overview Next: Molecular Orbital Theory →