University Quantum Mechanics II · The Variational Principle · 10.8

Linear variation: Rayleigh-Ritz as a matrix problem

This is where variational theory becomes the method everyone actually runs: pick a basis, build two matrices, diagonalise. You will derive H c = E S c, see why every root is an honest bound that can only fall, and watch ten oscillator states pin the x⁴ ground state at a coupling where perturbation theory gives nonsense.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

Every trial function in this unit so far has varied a shape: a width, an effective charge, two exponents. Rayleigh-Ritz varies something simpler, the coefficients of a fixed set of N functions, and that restriction pays for itself completely. Because ⟨H⟩ is then a ratio of two quadratic forms in the coefficients, its stationary points are the roots of a matrix problem, H c = E S c, which a computer solves exactly with N in the thousands.

There are N roots, not one, and each is an upper bound on the exact level of the same rank and symmetry; adding a function can only lower them, the old and new roots interleaving. Run on Unit 6's x⁴ oscillator at g = 1, where first-order perturbation theory gives 1.25 ħω and second order a meaningless −1.375 ħω, ten even oscillator states fix E₀ = 0.80377 ħω to five figures, by the same diagonalisation Unit 6 used as its benchmark. The same algebra in a basis built on r₁, r₂ and r₁₂ carries helium to within the physics its Hamiltonian leaves out.

It stops where the basis stops: it proves an inequality but never a lower bound, it converges slowly when the functions have the wrong scale or cannot form a cusp or a tail, it cannot represent a continuum state, and a nearly linearly dependent basis hands the last digits to round-off.

Simple definition: The Rayleigh-Ritz, or linear variation, method takes the trial state to be an arbitrary linear combination of N fixed functions, so that minimising ⟨H⟩ over the coefficients becomes the N × N matrix problem H c = E S c.
Example: For the x⁴ oscillator at g = 1, the even states |0⟩ and |2⟩ give a 2 × 2 matrix with elements 1.25, 2.121 and 12.25 ħω. Its lower root, 0.855 ħω, bounds E₀ = 0.8038 ħω from above; |0⟩ alone gave 1.25 ħω.

Energy of a linear trial stateE(c) = c†Hc / c†Sc, with Hᵢⱼ = ⟨χᵢ|H|χⱼ⟩ and Sᵢⱼ = ⟨χᵢ|χⱼ⟩

Every trial state in the span of the χᵢ is one vector c, so the search over functions becomes a search over N numbers, and ⟨H⟩ a ratio of two quadratic forms. Derived.

c: column of N complex coefficients; Hᵢⱼ in J (ħω here); Sᵢⱼ, the overlap matrix, dimensionless for normalised χᵢ

Stationarity: the generalised eigenvalue problemH c = E S c, solvable only where det(H − E S) = 0

Setting ∂E/∂cᵢ* = 0 gives it in one line, derived. Drop S for a non-orthogonal basis and an eigenvalue is no longer ⟨H⟩ in the state its eigenvector describes, so it bounds nothing.

N real roots ε₀(N) ≤ ε₁(N) ≤ …, in J; S Hermitian and positive definite; S = 1 for an orthonormal basis

Each root bounds the level of the same rankεₖ(N) ≥ Eₖ, for k = 0, 1, …, N − 1

Derived: the span of the lowest k + 1 eigenvectors contains a state orthogonal to the k exact levels below, where ⟨H⟩ ≥ Eₖ, and εₖ(N) is the largest ⟨H⟩ on that span.

Eₖ in J: the kth exact level, from k = 0, of the basis's symmetry block; in the even x⁴ block that is the oscillator's E₂ₖ

Interleaving when one function is addedεₖ(N + 1) ≤ εₖ(N) ≤ εₖ₊₁(N + 1)

No root can rise as the basis grows, and each old root is trapped between two new ones: adding |4⟩ at g = 1 turns 0.855 and 12.645 ħω into 0.808, 7.383 and 40.56 ħω.

Quoted: Hylleraas and Undheim (1930), MacDonald (1933); energies in J; true whatever function is added

The x⁴ oscillator as a banded matrix⟨m|H|n⟩ = ħω[(n + ½)δₘₙ + (g/4)⟨m|(a + a†)⁴|n⟩], zero unless m − n = 0, ±2 or ±4

Derived by counting ladder paths. Parity splits the matrix: the even states n = 0 to 18 form a 10 × 10 block with only five nonzero diagonals.

g = λħ/(m²ω³), dimensionless; ħω in J; diagonal term ¾g(2n² + 2n + 1), next ½g(2n + 3)√((n + 1)(n + 2))

The benchmark at g = 1ε₀(10) = 0.803774 ħω ≥ E₀ = 0.8037707 ħω

Five figures from a 10 × 10 matrix at a coupling where the series fails in the first figure. A basis of frequency 2ω, converged, agrees to nine figures.

Ten even states against eighty, computed here; ħω in J. First order, 1.25 ħω, is ε₀(1); second order gives −1.375 ħω

Fixed functions turn minimisation into a matrix problem

Choose N functions |χ₁⟩ to |χ_N⟩ and let the trial state be |ψ⟩ = Σ cᵢ|χᵢ⟩, with only the coefficients free. Then ⟨ψ|H|ψ⟩ = Σ cᵢ* Hᵢⱼ cⱼ and ⟨ψ|ψ⟩ = Σ cᵢ* Sᵢⱼ cⱼ, so E(c) = c†Hc / c†Sc is a ratio of two quadratic forms. Setting ∂E/∂cᵢ* = 0 for each i gives Σⱼ (Hᵢⱼ − E Sᵢⱼ) cⱼ = 0: the generalised eigenvalue problem H c = E S c, with solutions only where det(H − ES) = 0. That is derived, not assumed. H is Hermitian, and S is Hermitian and positive definite whenever the χᵢ are linearly independent, so the N roots are real; for an orthonormal basis S = 1 and the problem is ordinary diagonalisation. Each root equals E(c) at its own eigenvector, so the lowest root is ⟨H⟩ in a genuine trial state and the variational theorem already makes it at least E₀. Dimensions check themselves: Hᵢⱼ and E Sᵢⱼ carry the same units whatever normalisation the χᵢ have.

Why every root is a bound, and why none can rise

Adding a function enlarges the space the minimum is taken over, so the lowest root cannot rise. The upper roots need one more step, derived here. Rank the roots from k = 0. The root εₖ(N) is the largest value of ⟨H⟩ on the span of the lowest k + 1 eigenvectors, because on that span ⟨H⟩ is a weighted average of ε₀(N) to εₖ(N). A (k + 1)-dimensional space always contains a state orthogonal to the k exact eigenstates below level k, since that is only k conditions, and in that state ⟨H⟩ ≥ Eₖ. So εₖ(N) ≥ Eₖ for every k. Topic 10.03's warning still stands: a single state orthogonal to an approximate ground state proves nothing. The argument works because the root is a maximum over a subspace. The stronger statement, that adding one function makes old and new roots interleave, εₖ(N + 1) ≤ εₖ(N) ≤ εₖ₊₁(N + 1), is the theorem of Hylleraas and Undheim (1930) and MacDonald (1933), quoted here. It promises bounds that only fall; it gives no lower bound and no measure of the gap.

The quartic oscillator becomes a banded matrix

Take Unit 6's H = p²/2m + ½mω²x² + λx⁴ and, as the basis, the eigenstates |n⟩ of its oscillator part, so S = 1. Quoting x = x₀(a + a†)/√2 from Quantum Mechanics I, with x₀ = √(ħ/mω), gives λx⁴ = (g/4)ħω(a + a†)⁴ with g = λħ/(m²ω³); dimensions, (J m⁻⁴)(J s)/(kg² s⁻³) = 1, confirm g is a pure number. Each a or a† shifts n by one, so four of them reach only n, n ± 2 and n ± 4. Counting the paths gives ⟨n|X⁴|n⟩ = ¾(2n² + 2n + 1), ⟨n + 2|X⁴|n⟩ = ½(2n + 3)√((n + 1)(n + 2)) and ⟨n + 4|X⁴|n⟩ = ¼√((n + 1)(n + 2)(n + 3)(n + 4)), with X = x/x₀. Parity kills every element between even and odd n, so n = 0, 2, …, 18 form their own 10 × 10 block with five nonzero diagonals. At g = 1 its first row reads 1.25, 2.121 and 1.225 ħω and its second starts 2.121, 12.25 and 12.12 ħω. Limiting case: at g = 0 the matrix is diagonal and its roots are n + ½ exactly.

Ten states fix five figures where the series fixed none

At g = 1 Unit 6's series gives 1.25 ħω at first order and −1.375 ħω at second, below zero for a Hamiltonian that can never go below ħω/2, since λx⁴ ≥ 0. First order is not a lucky guess: E⁽⁰⁾ + E⁽¹⁾ = ⟨0|H|0⟩ is the N = 1 root, so it is a bound; second order is no bound at all. Diagonalising the even block, as computed here, the lowest root runs 1.25, 0.8551, 0.8082, 0.8059 and 0.8056 ħω for N = 1 to 5, then 0.804698, 0.804076, 0.803838 and 0.803781 for N = 6 to 9, and 0.803774 ħω at N = 10, against 0.8037707 ħω from eighty even states. Ten states fix 0.80377 to five figures, and every entry lies above the converged value, as the theorem demands. The second root falls from 12.64 through 7.38 and 5.86 to 5.1815 ħω at N = 10, bounding the n = 2 level, 5.1793 ħω; higher levels need more of the basis, so their roots converge more slowly. The odd block, n = 1, 3, …, gives E₁ = 2.7379 ħω the same way.

The basis scale is a free choice, and a good one saves states

Nothing forces the basis to share the Hamiltonian's frequency. Eigenstates of an oscillator of frequency Ω = sω still give S = 1 and a banded matrix; the diagonal becomes (2n + 1)(s + 1/s)/4 plus the quartic term scaled by 1/s², and each n, n + 2 element gains (1/s − s)√((n + 1)(n + 2))/4. One state then gives (s + 1/s)/4 + 3g/(4s²), smallest where s³ − s = 6g: at g = 1 that is s = 2 and 0.8125 ħω, the best bound a single Gaussian of Topic 10.02's kind can give. At that same s, ⟨2|H|0⟩ = √2[(1/s − s)/4 + 3g/(2s²)] vanishes, since setting it to zero is the same equation; stretching a Gaussian changes it, to first order, only by |2⟩. So N = 2 gains nothing at s = 2, yet N = 4 gives 0.803795 and N = 7 gives 0.803771 ħω, five figures with three fewer states than at s = 1. A basis whose functions already have the right size converges fast and one that does not converges slowly, but both approach from above.

What a converged number certifies, and where the method stops

Rayleigh-Ritz certifies an inequality, not an error bar: 0.803774 ħω ≥ E₀ is proved, while E₀ = 0.80377 ħω to five figures is an inference from roots that have stopped moving, confirmed here by a second basis, s = 2, settling on the same digits. The same machinery with a few thousand functions of r₁, r₂ and r₁₂ gives helium's fixed-nucleus energy, −79.014 eV (quoted), against the measured −79.005 eV, the sum of its two ionisation energies. The 0.009 eV left, about one part in 10⁴, is nuclear recoil, of relative size mₑ/M ≈ 1.4 × 10⁻⁴, plus smaller relativistic and QED corrections of relative order α² ≈ 5 × 10⁻⁵ and beyond, all quoted: physics missing from H, which no basis can supply. The method stops where the basis cannot follow the state. A continuum state is not a normalisable sum of bound functions, a cusp or a long tail needs many smooth functions, and near-duplicate functions make S nearly singular, so round-off rather than physics sets the last digits.

Change one variable at a time

Make the relationship visible.

Interactive model

Coupling g = λħ/(m²ω³)1.00

Basis frequency Ω/ω1.0

At s = 1 and g = 1, one state gives 1.25 ħω, first order exactly, three give 0.808 ħω, and second order has fallen to −1.375 ħω. Now raise s to 2: the N = 2 curve touches N = 1 at g = 1, where ⟨2|H|0⟩ vanishes, and all three sit within 0.009 ħω of the circle at 0.80377 ħω.

ONE STATE, N = 11.2500 ħω

TWO STATES, N = 20.8551 ħω

SECOND ORDER, SERIES-1.3750 ħω

COUPLING ⟨2|H|0⟩2.1213 ħω

Live interpretationONE STATE, N = 1: 1.2500 ħω. TWO STATES, N = 2: 0.8551 ħω. SECOND ORDER, SERIES: −1.3750 ħω. COUPLING ⟨2|H|0⟩: 2.1213 ħω

Catch the common trap

Explain before calculating.

For the x⁴ oscillator at g = 1, the even states |0⟩ and |2⟩ give roots 0.8551 and 12.645 ħω; adding |4⟩ gives 0.8082, 7.383 and 40.56 ħω. What does the Rayleigh-Ritz method guarantee?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyFor the x⁴ oscillator at g = λħ/(m²ω³) = 0.5, use the even oscillator states |0⟩ and |2⟩ as the basis. Build the 2 × 2 matrix of H in units of ħω, find both roots and say what each bounds, then compare with first- and second-order perturbation theory and with the converged E₀ = 0.69618 ħω.

Diagonal: ⟨n|H|n⟩/ħω = n + ½ + ¾g(2n² + 2n + 1). For n = 0 this is 0.5 + 0.375 = 0.875; for n = 2 it is 2.5 + ¾ × 0.5 × 13 = 2.5 + 4.875 = 7.375.
Off-diagonal: ⟨2|H|0⟩/ħω = ½g(2n + 3)√((n + 1)(n + 2)) at n = 0, which is 0.5 × 1.5 × √2 = 1.0607. The oscillator part is diagonal and S = 1, so this is the only coupling.
Roots: E = 4.125 ∓ √(3.25² + 1.0607²) = 4.125 ∓ √11.6875 = 4.125 ∓ 3.4187, giving 0.7063 and 7.5437 ħω.
Bounds: 0.7063 ħω ≥ E₀, and 7.5437 ħω bounds the second even level, n = 2, converged at 4.3275 ħω. The lower root is 0.0101 ħω, or 1.5%, above the converged 0.69618 ħω.
Series: first order is the one-state root, 0.875 ħω, 26% high. Second order gives 0.875 − (21/8)(0.5)² = 0.2188 ħω, 69% low and below ħω/2, which H ≥ H₀ forbids since λx⁴ ≥ 0.

AnswerE₀ ≤ 0.7063 ħω and the n = 2 level ≤ 7.544 ħω. The lower root is 1.5% above the converged 0.6962 ħω, against 26% high at first order and 69% low at second.

MediumAt g = 1 a student builds the x⁴ oscillator's matrix in the basis |0⟩, |2⟩ by squaring the 2 × 2 matrix of X² = x²/x₀² instead of computing ⟨m|X⁴|n⟩. Find the root this gives, compare it with the converged E₀ = 0.80377 ħω, and find which element is wrong and why the bound is lost.

From X² = (a + a†)²/2 = (a² + a†² + 2n + 1)/2: ⟨0|X²|0⟩ = ½, ⟨2|X²|2⟩ = 5/2 and ⟨2|X²|0⟩ = √2/2.
Squaring the 2 × 2 matrix: the (0, 0) element is ¼ + ½ = ¾ and the (2, 0) element is ½ × √2/2 + √2/2 × 5/2 = 3√2/2 = 2.1213, both equal to the true ⟨0|X⁴|0⟩ and ⟨2|X⁴|0⟩. But the (2, 2) element is ½ + 25/4 = 6.75, where the true ⟨2|X⁴|2⟩ = ¾ × 13 = 9.75.
The missing 3 is |⟨4|X²|2⟩|² = (√12/2)² = 3: the path 2 → 4 → 2 runs through a state outside the basis. In general, with P the projector on the basis and Q = 1 − P, PX⁴P − (PX²P)² = PX²QX²P, which is never negative, so the product always undershoots.
Roots with the wrong element: H₂₂ = 2.5 + 6.75 = 9.25 ħω, so E = 5.25 − √(4.00² + 2.1213²) = 5.25 − √20.5 = 5.25 − 4.5277 = 0.7223 ħω.
Verdict: 0.7223 ħω lies 0.0815 ħω below the converged 0.80377 ħω, which no Rayleigh-Ritz root can do; the true elements give 0.8551 ħω, above it. The squared matrix is not ⟨ψ|H|ψ⟩ for any ψ in the span, so the theorem never applied.

AnswerThe squared-X² matrix gives 0.722 ħω, 0.081 ħω below the converged 0.804 ħω, because its ⟨2|X⁴|2⟩ is 6.75 instead of 9.75. The exact elements give 0.855 ħω, a genuine upper bound.

HardFor the x⁴ oscillator at g = 1, with X = x/x₀ and H/ħω = −½ d²/dX² + ½X² + X⁴, take the non-orthogonal basis χ₁ = e^(−X²) and χ₂ = e^(−2X²). Build the normalised S and H, solve H c = E S c, and compare with the best single Gaussian, 0.8125 ħω, and the converged 0.80377 ħω. Then find the lowest eigenvalue of H alone and say why it is not a bound.

For e^(−aX²) and e^(−bX²), with c = a + b: ∫e^(−cX²)dX = √(π/c), the X² and X⁴ integrals are √(π/c)/(2c) and 3√(π/c)/(4c²), and the kinetic term, after one integration by parts, is ab√(π/c)/c. So S(a, b) = √(π/c) and H(a, b) = √(π/c)[ab/c + 1/(4c) + 3g/(4c²)].
Dividing each element by √(S(a, a)S(b, b)): the diagonal elements are each Gaussian's own ⟨H⟩, 0.8125 ħω for χ₁ and 1.109375 ħω for χ₂. The overlap is S₁₂ = √(π/3)/(π²/8)¹⁄⁴ = (2√2/3)¹⁄² = 0.970984, and, with the bracket at c = 3 equal to 2/3 + 1/12 + 1/12, H₁₂ = 1.023327 × 0.833333/1.053907 = 0.809153 ħω.
det(H − ES) = 0 reads (1 − S₁₂²)E² − (H₁₁ + H₂₂ − 2S₁₂H₁₂)E + (H₁₁H₂₂ − H₁₂²) = 0, that is 0.057191E² − 0.350527E + 0.246639 = 0.
Roots: E = [0.350527 ∓ √(0.122869 − 0.056422)]/0.114382 = (0.350527 ∓ 0.257773)/0.114382, giving 0.81091 and 5.3181 ħω.
Compare: 0.81091 ħω is 0.00159 ħω below the better single Gaussian and 0.00714 ħω above the converged 0.80377 ħω, as a bound must be; the upper root bounds the n = 2 level, 5.1793 ħω. The small 1 − S₁₂² = 0.057 warns that the two functions are nearly dependent.
Dropping S: the eigenvalues of H alone are 0.960938 ∓ √(0.148438² + 0.809153²) = 0.960938 ∓ 0.822656, so 0.1383 ħω, below ħω/2, which H ≥ H₀ forbids. Such an eigenvalue is c†Hc/c†c, but the state c describes has norm c†Sc, so the eigenvalue is not that state's ⟨H⟩; lying below ħω/2, it is ⟨H⟩ of no state at all.

AnswerE₀ ≤ 0.8109 ħω, 0.0016 ħω below one Gaussian and 0.0071 ħω above the converged 0.8038 ħω, and the n = 2 level ≤ 5.318 ħω. Ignoring S gives 0.138 ħω, below the floor ħω/2: not a bound.