Second-Order Perturbation Theory

The Møller–Plesset second-order method (MP2), the workhorse "cheap correlation" method in computational chemistry, is nothing but this formula applied with $\hat H^{(0)}$ taken as the sum of Fock operators and $\hat H'$ as the difference between the true electron–electron repulsion and the mean-field (Hartree–Fock) approximation to it — the excited "states" $k$ are Slater determinants with electrons promoted into virtual orbitals, and the sum over $k$ becomes a sum over double excitations (Brillouin's theorem kills singles). This single second-order correction typically recovers 80–90% of the correlation energy missing from Hartree–Fock, making MP2 geometry optimizations and single-point energies a standard first step in drug-discovery pipelines, reaction-mechanism screening, and materials modeling, before more expensive coupled-cluster methods are applied for final accuracy. Second-order theory is also the direct theoretical basis of static and dynamic polarizabilities, computed routinely for molecules to predict Raman intensities, refractive indices, and nonlinear optical response (hyperpolarizabilities extend the same expansion to third order); of London dispersion energies, the attractive $1/R^6$ tail responsible for graphite's interlayer binding, the stability of DNA base-stacking, and the packing of drug molecules in a protein binding pocket — all traceable to the second-order coupling between instantaneous charge fluctuations on two separated fragments. In vibrational spectroscopy, second-order (and higher) vibrational perturbation theory (VPT2) corrects harmonic vibrational frequencies for anharmonicity, giving IR and Raman band positions accurate enough to assign real experimental spectra rather than idealized harmonic ones. In magnetic resonance, the second-order hyperfine and quadrupolar shifts in EPR and NMR spectroscopy, and indirect spin–spin coupling constants ($J$-couplings) computed by response theory, are sum-over-states second-order expressions in disguise. In solid-state and materials chemistry, the second-order Jahn–Teller effect explains distortions in molecules and crystals (like the pyramidalization of certain main-group compounds, or ferroelectric off-centering in perovskites such as BaTiO₃) that have no orbital degeneracy to invoke the (first-order) Jahn–Teller theorem, but do have a low-lying excited state strongly coupled by a vibronic operator — precisely the mechanism captured by the $E_n^{(2)}$ sum.

Ask the tutor: Can you show me how the second-order formula changes when two nearly-degenerate excited states are both included, and how that connects to setting up degenerate perturbation theory?

Picture the unperturbed Hamiltonian $\hat H^{(0)}$ as furnishing a ladder of exact energy levels $E_n^{(0)}$ with wavefunctions $|n^{(0)}\rangle$ — for instance, the harmonic oscillator ladder, or the hydrogen atom's orbitals. Now switch on a small perturbation $\hat H'$: an anharmonic term in a vibrational potential, an external electric field, or the true electron–electron repulsion left out of the mean-field starting point. First-order perturbation theory tells you the diagonal matrix element $\langle n^{(0)}|\hat H'|n^{(0)}\rangle$ shifts the energy — this is the average of the perturbation over the unperturbed density, a purely classical-feeling correction with no change in the shape of the state. But $\hat H'$ generally also has nonzero off-diagonal matrix elements $\langle k^{(0)}|\hat H'|n^{(0)}\rangle$ connecting state $n$ to every other state $k$. These off-diagonal couplings do two things simultaneously: they distort the wavefunction (mixing in a bit of $|k^{(0)}\rangle$), and — as a byproduct of that distortion — they lower the energy further, at second order in $\hat H'$. The size of the mixing coefficient is controlled by two competing factors: how strong the coupling matrix element is, and how far away in energy the other state sits. States that are strongly coupled and nearly degenerate mix in heavily and dominate the sum; states that are weakly coupled or energetically remote contribute almost nothing.

The energy gaps in the denominator explain a great deal of chemistry. When a molecule has a low-lying excited state that couples strongly to the ground state via some operator (the dipole operator for an electric field, the vibronic coupling operator for a nuclear distortion), the second-order energy lowering is large — this is exactly the second-order Jahn–Teller effect, which drives many non-symmetric distortions in molecules that are not Jahn–Teller active in the strict (degenerate) sense. Conversely, a large HOMO–LUMO gap suppresses second-order mixing, which is part of why aromatic, closed-shell molecules with big gaps are comparatively rigid and hard to polarize.

We start, as always in Rayleigh–Schrödinger perturbation theory, from an exactly solvable Hamiltonian $\hat H^{(0)}$ with known orthonormal eigenstates,

where: $\hat H^{(0)}$ is the zeroth-order (solvable) Hamiltonian, $E_n^{(0)}$ and $|n^{(0)}\rangle$ its exact eigenvalues/eigenvectors, and $\delta_{mn}$ the Kronecker delta expressing orthonormality.

Reading: this is the entire input to the theory — no perturbation appears yet; it is simply the exactly solved reference problem (particle in a box, harmonic oscillator, hydrogen atom, or a Hartree–Fock determinant basis). We write the true Hamiltonian as $\hat H = \hat H^{(0)} + \lambda \hat H'$ with a bookkeeping parameter $\lambda$ eventually set to 1, and expand both the energy and wavefunction as power series,

where: $E_n^{(1)}, E_n^{(2)}$ are the first- and second-order energy corrections, and $|n^{(1)}\rangle, |n^{(2)}\rangle$ the corresponding wavefunction corrections, each order-by-order smaller than the last.

Reading: this is a Taylor expansion in the "strength" of the perturbation; each successive term is assumed much smaller than the one before it, which is precisely the validity condition for the whole approach. Substituting into the full Schrödinger equation $\hat H|n\rangle = E_n|n\rangle$ and collecting powers of $\lambda$ gives, at first order,

where: $\hat H'$ is the perturbing operator, and the equation must be solved order by order because it holds for every power of $\lambda$ independently.

Reading: projecting this equation onto $\langle n^{(0)}|$ gives the familiar first-order energy formula $E_n^{(1)} = \langle n^{(0)}|\hat H'|n^{(0)}\rangle$; projecting instead onto $\langle k^{(0)}|$ for $k\ne n$ is what yields the wavefunction correction, derived below. The first-order wavefunction correction, expanded in the complete unperturbed basis, is

where: the sum runs over all states $k$ other than $n$, $\langle k^{(0)}|\hat H'|n^{(0)}\rangle$ is the coupling matrix element between state $n$ and state $k$, and $E_n^{(0)} - E_k^{(0)}$ is the unperturbed energy gap.

Reading: the perturbed state is the original state plus a linear combination of every other unperturbed state, with mixing coefficients that grow with coupling strength and shrink with energy separation — exactly the "borrowing" picture above. Finally, the second-order energy correction, obtained by projecting the second-order equation onto $\langle n^{(0)}|$, is

where: $|\langle k^{(0)}|\hat H'|n^{(0)}\rangle|^2 \ge 0$ is always non-negative, and $E_n^{(0)} - E_k^{(0)}$ carries the sign.

Reading: for the ground state ($n=0$), every denominator $E_0^{(0)} - E_k^{(0)}$ is negative (since all other states lie higher), so $E_0^{(2)}$ is a sum of non-positive terms — the ground state can only be pushed down by second-order mixing, never up.

Example 1 — Foundational: the quartic-perturbed harmonic oscillator, second-order shift of the ground state. Consider a harmonic oscillator with perturbation $\hat H' = c\,\hat x^4$ where $c$ has units of energy/length⁴, a standard toy problem for anharmonicity. Using ladder operators, $\hat x = \sqrt{\hbar/(2m\omega)}(\hat a + \hat a^\dagger)$. The first-order energy shift of the ground state is well known from the first-order theory,

but here we want the second-order piece, which requires the off-diagonal matrix elements of $\hat x^4$. Expanding $\hat x^4$ in ladder operators, the only nonzero off-diagonal matrix elements out of $|0\rangle$ connect to $|2\rangle$ and $|4\rangle$:

Take asymmetric, concrete numbers: let $\hbar\omega = 0.15\ \text{eV}$ (a C–H bending-like mode) and let $\beta \equiv \hbar/(2m\omega) = 0.020\ \text{Å}^2$ (not a round number, chosen to avoid masking algebra), and $c = 0.05\ \text{eV/Å}^4$. Then:

so the coupling matrix element is $\langle 2|\hat H'|0\rangle = c \times 3.394\times10^{-3} = 0.05 \times 3.394\times10^{-3} = 1.697\times10^{-4}\ \text{eV}$.

so $\langle 4|\hat H'|0\rangle = 0.05 \times 1.1758\times10^{-2} = 5.879\times10^{-4}\ \text{eV}$. The unperturbed energy gaps are $E_0^{(0)} - E_2^{(0)} = -2\hbar\omega = -0.30\ \text{eV}$ and $E_0^{(0)} - E_4^{(0)} = -4\hbar\omega = -0.60\ \text{eV}$. Applying the second-order formula:

This is small (as expected, since $\beta \ll 1\ \text{Å}^2$ is a genuinely weak-perturbation regime) and — exactly as the sign rule of the second derivation guarantees — negative. Limiting-case check: as $c \to 0$ (perturbation switched off), both matrix elements vanish quadratically in $c$ so $E_0^{(2)} \propto c^2 \to 0$ faster than $E_0^{(1)}\propto c$, confirming that second order is a genuinely smaller correction that must be checked for convergence before being trusted.

Example 2 — Applied: quadratic Stark effect and polarizability of a two-state atomic model. A simplified atomic model retains just the ground state $|0\rangle$ ($2s$-like, $E_0^{(0)} = -14.0\ \text{eV}$) and one dipole-coupled excited state $|1\rangle$ ($2p_z$-like, $E_1^{(0)} = -9.4\ \text{eV}$), a caricature of the lowest $s\to p$ transition in an alkali-metal atom. Take the transition dipole matrix element to be $\langle 1|\mu_z|0\rangle = 2.20\ e a_0$ (non-round, representative of a strongly allowed atomic transition), and apply a static field $\mathscr{E} = 5.00\times10^{7}\ \text{V/m}$ (a large but laboratory-attainable field). The perturbation is $\hat H' = -\mu_z \mathscr{E}$, so the coupling matrix element is

Compute step by step: $e a_0 = (1.602\times10^{-19})(5.292\times10^{-11}) = 8.478\times10^{-30}\ \text{C}\,\cdot\,\text{m}$. Then $2.20\, e a_0 = 1.865\times10^{-29}\ \text{C}\,\cdot\,\text{m}$, and

The energy gap is $E_0^{(0)} - E_1^{(0)} = -14.0 - (-9.4) = -4.6\ \text{eV}$. The second-order shift is

Confirming the sign rule: negative, as required for a ground state. Converting to the polarizability using the two-derivation formula $\alpha_{zz} = 2|\langle1|\mu_z|0\rangle|^2/(E_1^{(0)}-E_0^{(0)})$:

which, converted to the conventional atomic-polarizability units ($\alpha/4\pi\varepsilon_0$, in ų), gives roughly $\alpha_{zz}/(4\pi\varepsilon_0) \approx 8.49\ \text{Å}^3$ — squarely in the range of real alkali-atom polarizabilities (sodium is about $24\ \text{Å}^3$; this stripped-down two-level model captures the right order of magnitude despite ignoring all other $np$ states). Limiting-case check: as the gap $E_1^{(0)}-E_0^{(0)} \to \infty$ (a very "stiff," hard-to-excite atom), $\alpha_{zz}\to 0$ — a rigid atom cannot be polarized, exactly as physically expected.

Example 3 — Exam/competition level: three-state system, computing $E_0^{(2)}$ and testing near-degeneracy breakdown. Consider a model three-level system (e.g., a truncated basis for a conjugated π-system's frontier orbitals) with unperturbed energies $E_0^{(0)} = 0$, $E_1^{(0)} = 1.8\ \text{eV}$, $E_2^{(0)} = 2.1\ \text{eV}$ (states 1 and 2 are close together — a near-degeneracy — while state 0 is well separated). The perturbation matrix elements out of the ground state are $\langle 1|\hat H'|0\rangle = 0.35\ \text{eV}$ and $\langle 2|\hat H'|0\rangle = 0.12\ \text{eV}$ (asymmetric, unequal couplings by design). Applying the sum-over-states formula directly to the ground state:

So the ground state is stabilized by about $0.0749\ \text{eV}$, dominated (by a factor of roughly 10) by mixing with state 1 despite state 2 having a smaller gap, because the coupling to state 1 is nearly three times larger — illustrating that coupling strength, not gap size alone, controls which state dominates the sum, exactly as the two-panel figure above shows. Now probe excited-state 1's second-order energy, where the sign rule does not guarantee stabilization. State 1 couples to state 0 (gap $E_1^{(0)}-E_0^{(0)} = +1.8\ \text{eV}$, giving a positive-denominator, destabilizing contribution) and to state 2 (gap $E_1^{(0)}-E_2^{(0)} = 1.8-2.1 = -0.3\ \text{eV}$, a small negative denominator — a near-degeneracy). Suppose $\langle 2|\hat H'|1\rangle = 0.09\ \text{eV}$. Then

Here the first term is positive (destabilizing — borrowing from a lower state pushes an excited state up) while the second, from the near-degenerate partner, is large and negative because the small denominator amplifies it. If the gap $E_1^{(0)}-E_2^{(0)}$ were shrunk further, say to $-0.03\ \text{eV}$ instead of $-0.3\ \text{eV}$ (a ten-fold closer near-degeneracy) with the same coupling, that single term would balloon to $0.0081/(-0.03) = -0.27\ \text{eV}$ — comparable to or larger than the entire unperturbed gap itself, signaling that the perturbative series is no longer trustworthy and that states 1 and 2 must instead be treated by degenerate (quasi-degenerate) perturbation theory, diagonalizing the $2\times2$ block exactly rather than expanding it order by order.