“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Linear Response Theory

Perturbation Methods

⛓ What you need to know first 7 concepts, 4 layers

The requisite-knowledge inventory for this page, bottom-up: the primitives at the base, combined upward until you reach what this page assumes. Skim the layers you already own; start wherever the ground gets unfamiliar.

base
- Linear algebraconcept
- Quantum mechanics (states & operators)
L1
- Matrix eigenvalue problems
- Perturbation theory
- The eigenvalue problem (Ax = λx)concept
L2
- Expanding in the eigenbasis
L3
- The resolvent & spectral function
↳you are here

2 of these are concepts without a dedicated page yet — the grey chips. Following the linked ones first makes the rest land.

Kick a system gently and watch what comes back. The response is not mush — it rings, and it rings only at the system's own transition frequencies. Linear response theory is the bookkeeping for that fact: one function $χ (ω)$ whose poles are the excitation energies and whose residues are how brightly each one couples to your probe. Every spectroscopy is a rediscovery of this sentence, and Casida's equation is what it turns into for molecules. Five ways in; read the one that lands.

read it as

Strike a bell once and it does not hum at the frequency of your mallet — it rings at its own frequencies, the ones fixed by its geometry before you ever walked up to it. Record the ring, Fourier transform it, and you have measured the bell's spectrum without ever solving for a single vibrational mode.

Quantum systems play the same game. Give the ground state a sudden weak kick (in the code below, literally multiply by $e^{iκ \overset{x}{^}}$ and let it evolve) and the observable $⟨ \overset{x}{^} (t)⟩$ oscillates as a superposition of pure tones — one per excited state the kick can reach, each at exactly $ω_{n} = E_{n} - E_{0}$ . The spectrum was sitting in the dynamics the whole time. Linear response theory is the guarantee that the ring you record and the eigenvalues you would have computed are the same information.

One formula, three altitudes

The sum-over-states formula, with real orbitals and the two time-reversal terms combined over a common denominator:

χ (ω) = n \sum \frac{2 ω _{n} s _{n}}{ω _{n}^{2} - ( ω + i η ) ^{2}}, s_{n} = ∣ ⟨ n ∣ \overset{x}{^} ∣0 ⟩ ∣^{2} .

Now bring it down to earth. For a two-level system — one ground state, one bright excited state at gap $ω_{10}$ — the sum has a single term:

χ (ω) = \frac{2 ω _{10} ∣ x _{01} ∣ ^{2}}{ω _{10}^{2} - ( ω + i η ) ^{2}} .

Set $ω = 0$ and you get the static polarizability $χ (0) = 2∣ x_{01} ∣^{2} / ω_{10}$ — push gently and constantly, and the system leans by that much per unit force. Sweep $ω$ up toward $ω_{10}$ and the denominator heads for zero: resonance. Past it, the sign flips and the response lags the drive. One pole, and already the whole phenomenology of a driven oscillator.

Finally, real numbers from a real Hamiltonian. Take the quartic anharmonic oscillator $H = \frac{p ^{2}}{2} + \frac{x ^{2}}{2} + 0.1 x^{4}$ — same conventions as the resummation pages — and diagonalize it in a 120-function basis (the script below; levels converged to $1 0^{- 13}$ ). The bright poles of $χ$ come out

(ω_{1}, s_{1}) = (1.210356, 0.412282), (ω_{3}, s_{3}) = (4.069736, 2.43 \times 1 0^{- 4}),

and the single-pole approximation of the formula above becomes $χ (ω) \approx 0.998/ (1.465 - ω^{2})$ , giving a static polarizability of $0.6813$ — within $1 0^{- 4}$ of the full answer $0.68138$ , because the higher poles are 1,700 and 2,000,000 times dimmer. At $λ = 0$ the same machinery returns exactly one pole at $ω = 1$ with $s = 1/2$ , so $χ (ω) = 1/ (1 - ω^{2})$ — the classical Lorentz oscillator, and $χ (0) = 1$ is just Hooke's law $x = F / k$ with $k = 1$ . The quantum formula knows classical mechanics when the potential is quadratic.

Watch the poles, on real numbers

Every curve below is assembled from the exact poles and residues printed by the validation script — nothing is synthesized at runtime except the Lorentzian broadening, which is the same $η$ that damps the time signal in the code. Switch the anharmonicity, then drag $η$ down and watch the peaks sharpen toward delta functions parked at the excitation energies.

linewidth η = 0.050

Exact poles for λ = 0.1: ω = 1.2104 (strength 0.4123) · ω = 4.0697 (strength 2.4e-4). Each peak height is strength/η — shrink η and the peaks march toward delta functions parked at the poles.

Three things to read off. First, the harmonic oscillator shows one line — out of infinitely many excited states, only $n = 1$ is bright, because $\overset{x}{^}$ only reaches one rung up the ladder. Second, turning on anharmonicity does two separate things: it shifts the main line (from $ω = 1$ to $1.21$ at $λ = 0.1$ , to $1.63$ at $λ = 0.5$ ) and it switches on a new line at $ω_{3}$ — a transition that was exactly dark by the harmonic selection rule, now allowed at three orders of magnitude down. Third, the linewidth is not a property of the system at all: it is $η$ , your resolution. The poles underneath are infinitely sharp; a finite observation window (or a finite excited-state lifetime, in a real molecule) is what smears them.

Two ways to compute it — and why they must agree

The claim "poles of the response ARE the excitation energies" deserves a numerical cross-examination, so the script computes $χ (ω)$ twice with no shared machinery. Way one is the Lehmann formula: diagonalize, read off $ω_{n}$ and $s_{n}$ , sum the fractions.

import numpy as np

# H = p^2/2 + x^2/2 + lam x^4 in a harmonic-oscillator basis (hbar = m = 1)
def aho_eigen(lam, nbasis=120):
    n = np.arange(nbasis)
    x = np.zeros((nbasis, nbasis))
    off = np.sqrt((n[:-1] + 1) / 2.0)      # <n|x|n+1> = sqrt((n+1)/2)
    x[n[:-1], n[:-1] + 1] = off
    x[n[:-1] + 1, n[:-1]] = off
    H = np.diag(n + 0.5) + lam * np.linalg.matrix_power(x, 4)
    E, U = np.linalg.eigh(H)
    return E, U.T @ x @ U                  # energies, x in the eigenbasis

E, x_eig = aho_eigen(lam=0.1)
w = E[1:] - E[0]                           # poles:    excitation energies
s = x_eig[0, 1:] ** 2                      # residues: |<n|x|0>|^2

def chi_sos(omega, eta=0.05):              # the sum-over-states formula
    return np.sum(s * 2 * w / (w**2 - (omega + 1j * eta)**2))

Way two never sees an eigensolver. Cool a Gaussian to the ground state in imaginary time, kick it with $e^{iκ x}$ , propagate on a grid with the split-operator method, record $⟨ x (t)⟩$ , and Fourier transform the transient. First-order perturbation theory in $κ$ says the kicked signal is the response kernel: $δ ⟨ x (t)⟩ / κ = 2 \sum_{n} s_{n} sin (ω_{n} t) = χ (t)$ — so its transform must land on the same poles.

# Independent check: no eigensolver anywhere. Ground state by imaginary
# time, then a delta-kick exp(i kappa x), then split-operator propagation.
xs = np.linspace(-L, L, N, endpoint=False); dx = xs[1] - xs[0]
k = 2 * np.pi * np.fft.fftfreq(N, d=dx)
V = 0.5 * xs**2 + lam * xs**4

psi = np.exp(-xs**2 / 2).astype(complex)   # imaginary time: cool to |0>
eV, eT = np.exp(-V * dt / 2), np.exp(-0.5 * k**2 * dt)
for _ in range(20000):
    psi = eV * np.fft.ifft(eT * np.fft.fft(eV * psi))
    psi /= np.sqrt(np.sum(np.abs(psi)**2) * dx)

psi = np.exp(1j * kappa * xs) * psi        # the kick, kappa = 1e-3
expV = np.exp(-1j * V * dt / 2)            # real time: listen
expT = np.exp(-1j * 0.5 * k**2 * dt)
sig = []
for j in range(int(T / dt)):               # T = 400, dt = 0.002
    psi = expV * np.fft.ifft(expT * np.fft.fft(expV * psi))
    sig.append(np.sum(xs * np.abs(psi)**2).real * dx)

# d<x(t)>/kappa is chi(t); a damped Fourier transform gives chi(omega)
resp = (np.array(sig) - 0.0) / kappa
chi_grid = fourier(resp * np.exp(-eta * t))   # e^{-eta t} <-> Lorentzian width

They agree to six digits across the whole frequency window — the residual is the Trotter error of the propagator, not physics. The Thomas-Reiche-Kuhn sum rule $\sum_{n} ω_{n} s_{n} = 1/2$ holds to ten digits for both values of $λ$ , and the harmonic limit recovers its single exact pole. (The grid run's weight for the weak $ω_{3}$ line reads high — $3.3 \times 1 0^{- 4}$ against the exact $2.4 \times 1 0^{- 4}$ — because at $η = 0.05$ the peak sits on the Lorentzian tail of a line 1,700 times brighter. Shrink $η$ in the demo and watch that contamination die.)

basis convergence (lam=0.1): max |E_80 - E_120| over first 8 levels = 4.37e-13

lam = 0.0: E0 = 0.500000
  poles (sum over states):
    w = 1.000000   |<n|x|0>|^2 = 5.000000e-01
  TRK sum rule  sum w_n s_n = 0.5000000000   (exact: 0.5)
  grid vs sum-over-states: max relative deviation of chi(w) = 4.12e-06
  peaks found in the grid spectrum:
    w = 1.0000   weight = 4.996346e-01

lam = 0.1: E0 = 0.559146
  poles (sum over states):
    w = 1.210356   |<n|x|0>|^2 = 4.122817e-01
    w = 4.069736   |<n|x|0>|^2 = 2.434491e-04
    w = 7.340621   |<n|x|0>|^2 = 1.956422e-07
  TRK sum rule  sum w_n s_n = 0.5000000000   (exact: 0.5)
  grid vs sum-over-states: max relative deviation of chi(w) = 6.20e-06
  peaks found in the grid spectrum:
    w = 1.2104   weight = 4.121054e-01
    w = 4.0694   weight = 3.324795e-04

Full script: scripts/gen_linear_response.py. Everything above — poles, strengths, both spectra, the sum rule — is its verbatim output.

Try First

Each prompt asks a checkable question about the working code or math above — predict an output, derive a sign, state an invariant, find a bug. Commit to an answer before clicking "reveal." That commitment is the whole point: if your answer matched, you understand the piece you were looking at; if it didn't, that's the part worth re-reading.

predict

The harmonic oscillator has infinitely many excited states, evenly spaced at

ω = 1, 2, 3, \dots

above the ground state. Its absorption spectrum has how many lines? Predict, then argue from

s_{n} = ∣ ⟨ n ∣ \overset{x}{^} ∣0 ⟩ ∣^{2}

why does this work

Turning on

λ x^{4}

breaks the harmonic selection rule: the

0 \to 3

line switches on (strength

2.4 \times 1 0^{- 4}

λ = 0.1

). But

0 \to 2

and

0 \to 4

stay exactly dark at every value of

λ

the script tries. Why do some selection rules die and others survive?

invariant

The output shows

\sum_{n} ω_{n} s_{n} = 0.5000000000

for

λ = 0

and for

λ = 0.1

— ten digits, no drift. The poles moved, the strengths moved, their weighted sum did not. What protects it?

Where this goes

Everything downstream of this page is a change of arena, not of idea. Swap $\overset{x}{^}$ for the dipole operator and $∣0 ⟩$ for a Kohn-Sham ground state, and demanding that the density response be consistent with the response of the exchange-correlation potential turns the pole condition into a matrix eigenproblem in the space of occupied→virtual orbital pairs — that is Casida's equation, and its $(X, Y)$ two-block structure is the absorption/emission pairing from the Lehmann representation. When that matrix is too large to build, the resolvent framing takes over: kick, Lanczos, read the spectrum off a tridiagonal — Liouville-Lanczos. And the residues, dressed with a factor $\frac{2}{3} ω_{n}$ and summed against the TRK rule, become the oscillator strengths an experimental spectrum is plotted in.