First-Order Perturbation Theory

First-order perturbation theory is the invisible workhorse behind an enormous share of modern spectroscopy and quantum chemistry software. Every NMR chemical shift calculation treats the external magnetic field and electron shielding as perturbations on the bare nuclear spin Hamiltonian. Every prediction of a molecule's fine and hyperfine structure — spin-orbit coupling in transition-metal complexes, zero-field splitting in EPR spectroscopy of radicals and metalloenzymes — starts from an unperturbed electronic Hamiltonian and adds the relativistic spin-orbit term as $\hat H'$. Vibrational spectroscopists routinely use first-order (and second-order) perturbation theory to convert harmonic force fields from density functional theory into anharmonic correction factors that match observed infrared and Raman overtone frequencies, a step built into essentially every modern quantum chemistry package (Gaussian, ORCA, Q-Chem, Psi4). In drug design, perturbation-theory-based free energy methods (thermodynamic integration and free energy perturbation, FEP) use the same $\lambda$-scaling logic introduced here — gradually turning on a "perturbing" difference between two ligands' force fields — to predict relative binding affinities computationally, guiding lead optimization campaigns in the pharmaceutical industry. In materials science, crystal field theory and ligand field theory treat the electrostatic environment around a transition-metal ion as a perturbation on the free-ion atomic states, directly explaining d-orbital splitting patterns, color, and magnetism in coordination compounds and pigments. Møller-Plesset perturbation theory (MP2, MP3) remains one of the most widely used post-Hartree-Fock methods for adding electron correlation energy to molecular calculations, built entirely on the Rayleigh-Schrödinger expansion derived above, with the Hartree-Fock Hamiltonian playing the role of $\hat H_0$ and the residual electron-electron interaction (beyond the mean-field average) playing the role of $\hat H'$.

Ask the tutor: How do I know, before doing any calculation, whether a perturbation is "small enough" for first-order theory to be trustworthy for my specific system?

Imagine an electron trapped in a box whose right wall has been nudged slightly inward, or a hydrogen atom placed in a weak external electric field, or two electrons in helium that mutually repel. In each case, there exists a nearby "reference" problem — the perfect box, the free hydrogen atom, the helium atom without electron-electron repulsion — that we can solve in closed form. The real system's Hamiltonian is

where $\hat H_0$ is the reference (unperturbed) Hamiltonian and $\hat H'$ is everything else. Perturbation theory is the statement that, provided $\hat H'$ doesn't change the physics too violently, the true eigenstates and eigenvalues of $\hat H$ are close to those of $\hat H_0$, and we can compute "close to how much" systematically. The key conceptual move is to imagine turning the perturbation on gradually, governed by a parameter $\lambda$ that runs from 0 (no perturbation, just $\hat H_0$) to 1 (full perturbation, the real system). We then Taylor-expand the true energy and wavefunction in powers of $\lambda$:

At $\lambda = 0$ we recover the known, solved problem. The coefficient $E_n^{(1)}$ — the subject of this article — is the slope of the energy with respect to how strongly the perturbation is switched on, evaluated right at the unperturbed state. Physically, this captures the fact that if the perturbation shifts the potential energy landscape a little, the energy of a state shifts by the average amount that landscape changed, weighted by where the original (unperturbed) electron density sits. The wavefunction itself also deforms, but that deformation is a higher-order effect on the energy — this is a manifestation of the same stationarity that underlies the variational principle: at an energy extremum, first-order changes in the wavefunction produce only second-order changes in the energy. The diagram below shows the conceptual flow: an exactly solvable reference system, a small perturbing operator, and the resulting shifted (but not re-solved) energy level.

A second useful picture is geometric: think of the unperturbed wavefunction as a fixed probability density $|\psi_n^{(0)}|^2$ draped over the perturbing potential $H'(x)$. The first-order energy shift is exactly the weighted average of $H'(x)$ under that density curve — literally an area-under-the-curve calculation, no new Schrödinger equation required.

The full non-degenerate Rayleigh-Schrödinger perturbation expansion writes both the energy and the wavefunction of the $n$th state as power series in $\lambda$:

where:

• $\hat H_0$ has known, exactly solvable eigenstates $\hat H_0 \psi_n^{(0)} = E_n^{(0)}\psi_n^{(0)}$,
• $\lambda$ is a bookkeeping parameter set to 1 at the end (the physical perturbation strength),
• $\psi_n^{(1)}, \psi_n^{(2)},\dots$ are wavefunction corrections expressed as linear combinations of the other unperturbed eigenstates $\psi_k^{(0)}$ ($k\neq n$).

Substituting into $\hat H\psi_n = E_n\psi_n$ and collecting terms order by order in $\lambda$ yields, at zeroth order, simply $\hat H_0\psi_n^{(0)} = E_n^{(0)}\psi_n^{(0)}$ (trivially true), and at first order the central working equation of this article:

where:

• $\psi_n^{(0)}$ is the unperturbed, already-known wavefunction of the state of interest,
• $\hat H'$ is the full (physical, $\lambda=1$) perturbing operator,
• the integral runs over all coordinates of the system ($d\tau$ = volume element, possibly including spin sums).

Reading: the first-order correction to the energy of state $n$ is the expectation value of the perturbing operator, computed using the old, unperturbed wavefunction — you are not required to find the new wavefunction at all to get this leading correction to the energy. The total energy through first order is then

The companion result — needed if you want the wavefunction correction, or if you want to go on to second order — expands $\psi_n^{(1)}$ in the complete basis of unperturbed states:

where:

• the sum excludes $k=n$ (that piece is absorbed into normalization convention),
• $E_n^{(0)} - E_k^{(0)}$ in the denominator is the unperturbed energy gap between state $n$ and state $k$,
• this correction is used to build $E_n^{(2)}$, the leading topic of second-order perturbation theory.

Reading: the wavefunction "borrows" character from every other unperturbed state, with a weight that grows the stronger the coupling matrix element $\langle k|\hat H'|n\rangle$ is, and grows the closer in energy state $k$ lies to state $n$. This is exactly why near-degeneracies ($E_n^{(0)} \approx E_k^{(0)}$) blow up the expansion and demand a different treatment.

Example 1 — Foundational: asymmetric tilt in an infinite square well. Consider a particle of mass $m$ in an infinite square well of width $L$, with unperturbed eigenstates

Now perturb it with a linear "tilt" potential representing, say, a weak uniform field pulling the particle toward $x=L$:

We want the first-order shift in the ground state ($n=1$) and the second excited state ($n=3$) — deliberately asymmetric choices, not a "clean" case. The first-order correction is

Using the standard integral $\int_0^L x\sin^2(n\pi x/L)\,dx = L^2/4$ (true for every integer $n$, since the sine-squared term integrates to zero by symmetry about $x=L/2$, leaving only the linear-ramp average), we get

Remarkably, every level shifts by exactly $V_0/2$ at first order — the average value of a linear ramp over the box is always its midpoint value, regardless of quantum number, because $\sin^2(n\pi x/L)$ is symmetric about $x=L/2$ for all $n$. So for $n=1$: $E_1 \approx E_1^{(0)} + V_0/2$; for $n=3$: $E_3 \approx E_3^{(0)} + V_0/2$.

Numerical check. Take $L = 0.50\ \text{nm}$, $m = m_e$, and $V_0 = 0.10\ \text{eV}$ (chosen deliberately small compared to $E_1^{(0)}$, and non-round to stress-test the arithmetic). First,

Evaluating: $\hbar^2\pi^2 = (1.0546\times10^{-34})^2 \times 9.8696 = 1.0972\times10^{-67}$. Denominator: $2\times9.109\times10^{-31}\times(0.5\times10^{-9})^2 = 2\times9.109\times10^{-31}\times2.5\times10^{-19} = 4.5545\times10^{-49}$. So

Then $E_3^{(0)} = 9 \times 1.504\ \text{eV} = 13.54\ \text{eV}$. Both shift by $V_0/2 = 0.050\ \text{eV}$:

Limiting check. As $V_0 \to 0$ the correction vanishes and $E_n \to E_n^{(0)}$ exactly, as it must — with no perturbation there is nothing to correct. The correction is also independent of $n$ here only because $\hat H'$ is linear; a quadratic perturbation would not give an $n$-independent shift, since $\langle x^2\rangle$ in a box does depend on $n$.

Example 2 — Applied: helium's ground-state energy from electron-electron repulsion. Model helium's two electrons as independently occupying hydrogen-like $1s$ orbitals with nuclear charge $Z=2$ (this defines $\hat H_0$, exactly solvable since it's two independent hydrogen-like problems), and treat the electron-electron repulsion as the perturbation:

The zeroth-order ground state is a product of two hydrogen-like $1s$ orbitals with $Z=2$:

and the zeroth-order energy (sum of two hydrogen-like $1s$ energies, each $-Z^2\times 13.606$ eV) is

This is far below the experimental ground-state energy of helium, $-79.0$ eV, because $\hat H_0$ completely neglects the repulsion that pushes the electrons apart and raises the energy. The first-order correction is the standard textbook integral (evaluated exactly using hydrogen-like $1s$ orbitals):

(the standard $5Z/8$ Ry result derivable by direct integration of $1/r_{12}$ over the product of two $1s$ densities, expanded in Legendre polynomials — a calculation that belongs to the theory of two-electron integrals but whose use here is exactly the first-order formula above). Adding this correction:

Comparison to experiment. The true ground-state energy of helium is $-79.0$ eV. First-order perturbation theory gets us from a crude $-108.8$ eV to $-91.8$ eV — a large improvement (closing roughly 68% of the 29.8 eV gap), but not exact, because $\hat H'$ (repulsion $\sim 34$ eV out of $\sim 109$ eV, i.e. roughly 30% of $|\hat H_0|$) is not actually all that small. The residual 12.8 eV discrepancy is exactly the kind of correction that Second-Order Perturbation Theory and, in practice, full variational optimization of an effective nuclear charge, are built to capture.

Limiting check. If we let the repulsion strength be scaled by a fictitious parameter $\lambda$ ($\hat H' \to \lambda e^2/4\pi\epsilon_0 r_{12}$) and take $\lambda \to 0$, we recover exactly the unperturbed value $-108.8$ eV, confirming the correction vanishes continuously as the perturbation is switched off — exactly consistent with the general expansion $E(\lambda) = E^{(0)} + \lambda E^{(1)} + \cdots$.

Example 3 — Exam-level: quadratic Stark effect setup meets anharmonic correction (multi-part). (a) Why the linear Stark shift vanishes for a harmonic oscillator dipole in a field, forcing us to examine parity. Consider a charged harmonic oscillator (charge $q$, force constant $k$, mass $m$) placed in a uniform external electric field $\mathcal E$ along $x$. The perturbation is

The first-order energy correction for the $v$th vibrational level is

Harmonic oscillator eigenstates have definite parity: $\psi_v^{(0)}(-x) = (-1)^v \psi_v^{(0)}(x)$. Since $x$ is an odd operator, $|\psi_v^{(0)}|^2 \, x$ is an odd function of $x$ for any $v$ (because $|\psi_v^{(0)}(-x)|^2 = |\psi_v^{(0)}(x)|^2$ is even, times $x$ odd = odd overall), so the integral over all $x$ vanishes identically:

This is a rigorous, general result — not an approximation — following purely from the symmetry (parity) of harmonic oscillator wavefunctions. It is exactly why the Stark effect for a purely harmonic dipole is a second-order phenomenon (the famous quadratic Stark shift), foreshadowing the machinery of Second-Order Perturbation Theory. (b) Now add a genuine first-order-active perturbation: a cubic anharmonicity term, evaluated for specific asymmetric numbers. Take instead

Because $x^3$ is also odd, by the identical parity argument $\langle v^{(0)}|x^3|v^{(0)}\rangle = 0$ for every $v$ — the cubic anharmonic term likewise gives zero first-order shift, which is the textbook reason cubic anharmonicity is treated at second order in vibrational spectroscopy (it's what generates the anharmonic red-shift of overtones). (c) A first-order-active perturbation: quartic anharmonicity. Now take the physically standard next term in a Morse-like expansion,

with harmonic oscillator parameters representative of a C–H stretch: $m = 1.628\times10^{-27}$ kg (reduced mass of C–H), $k = 560\ \text{N/m}$, giving $\omega = \sqrt{k/m}$. Compute:

The needed expectation value is the standard harmonic oscillator result

For $v=2$ (an asymmetric, non-ground-state choice): $6(4) + 6(2) + 3 = 24+12+3=39$. The prefactor:

Squaring: $(5.521\times10^{-23})^2 = 3.048\times10^{-45}\ \text{m}^4$. So

The first-order energy correction is

which is unphysically large in SI joules for a single quantum (as expected since $c$ was chosen only for illustrative magnitude) — converting to more natural spectroscopic units, $5.94\times10^{-9}$ J corresponds to $\approx 3.7\times10^{10}$ eV, signaling that this particular $c$ is far too large to be a "small" perturbation for this oscillator; a realistic quartic anharmonicity constant would be many orders of magnitude smaller. This is deliberately included as a check of scope: first-order perturbation theory produces a numerical answer regardless of input, but the answer is only trustworthy when $E^{(1)} \ll \hbar\omega$. Here $\hbar\omega = 1.0546\times10^{-34}\times5.865\times10^{14} = 6.186\times10^{-20}$ J, and our computed correction ($5.94\times10^{-9}$ J) is about $10^{11}$ times larger — an immediate red flag that this $c$ value is unphysical for a real molecular bond and the perturbative expansion would not converge.

Limiting check. Setting $c \to 0$ recovers $E_2^{(1)} \to 0$ and the pure harmonic energy $E_2^{(0)} = (2+\tfrac12)\hbar\omega$ survives untouched, exactly as required.

Erwin Schrödinger developed much of the mathematical apparatus of quantum perturbation theory in his celebrated series of 1926 papers, adapting a classical technique — the perturbation methods used for centuries in celestial mechanics to compute how planets nudge each other's orbits — into the new quantum formalism. The specific "Rayleigh-Schrödinger" name honors an even earlier debt: Lord Rayleigh had used essentially the same expansion technique decades earlier to study the vibrations of mechanical and acoustic systems, and Schrödinger recognized that his new wave equation admitted an exactly analogous treatment for perturbed quantum states.