Perturbation Theory Benchmarks

How far can a low-order energy expansion be trusted? This calculation compares the quartic oscillator’s first two perturbative orders with numerical diagonalization and an optimized Gaussian bound. The code runs in Python with NumPy; it prints its library versions and checks the matrix construction before producing a table.

Required background. The worked perturbation calculation derives the coefficients. Rayleigh–Ritz explains why an exactly projected Hamiltonian gives upper bounds.

Quartic Energies in Oscillator Units

Measure position in $\sqrt{\hbar/(m\omega)}$ and energy in $\hbar\omega$ . Then

h=\frac{p^2+q^2}{2}+gq^4,\qquad g=\frac{\lambda\hbar}{m^2\omega^3}\ge0.

The ground-state approximations are

\epsilon_0^{[1]}=\frac12+\frac34g,\qquad \epsilon_0^{[2]}=\frac12+\frac34g-\frac{21}{8}g^2.

For the Gaussian bound, solve $y^3-y-6g=0$ for $y\ge1$ and use $\epsilon_G=(3y+1/y)/8$ . Its derivation is in the variational worked example.

Multiply Operators Before Projecting

Let $P_N$ retain oscillator states $0,\ldots,N-1$ . The Ritz matrix contains $P_Nq^4P_N$ , which is generally different from $(P_NqP_N)^4$ . The latter discards paths that temporarily leave the retained space before returning.

Construct $q=(a+a^\dagger)/\sqrt2$ in $N+4$ states, form its fourth power, and only then take the $N\times N$ block. Four applications change oscillator number by at most four, so this padding includes every path needed. Build $h_0$ directly from its exact diagonal energies.

import sys
import numpy as np

def quartic_matrix(N, g):
    a = np.diag(np.sqrt(np.arange(1, N + 4)), 1)
    q = (a + a.T) / np.sqrt(2)
    q4 = np.linalg.matrix_power(q, 4)[:N, :N]
    n = np.arange(N)
    # Independent oscillator moment, including the highest retained state.
    assert np.allclose(np.diag(q4), 3 * (2*n*n + 2*n + 1) / 4)
    assert np.allclose(q4[(n[:, None] + n) % 2 == 1], 0)
    return np.diag(n + 0.5) + g * q4

def gaussian(g):
    lo, hi = 1.0, 1.0 + (6*g)**(1/3)
    for _ in range(70):
        y = (lo + hi) / 2
        if y**3 - y > 6*g:
            hi = y
        else:
            lo = y
    y = (lo + hi) / 2
    return (3*y + 1/y) / 8

print('Python', sys.version.split()[0], 'NumPy', np.__version__)
print('g       numerical       first          second         Gaussian       N-change')
for g in (0, 0.001, 0.01, 0.1, 1, 10):
    energies = []
    for N in (80, 160, 240):
        H = quartic_matrix(N, g)
        E, states = np.linalg.eigh(H)
        residual = np.linalg.norm(H @ states[:, 0] - E[0]*states[:, 0])
        assert residual < 1e-8
        energies.append(E[0])
    change = abs(energies[-1] - energies[-2])
    # Nested Ritz spaces decrease toward the ground energy, up to roundoff.
    assert np.all(np.diff(energies) <= 1e-9)
    assert change < 1e-8
    first, second = 0.5 + 0.75*g, 0.5 + 0.75*g - 2.625*g*g
    print(f'{g:<7g} {energies[-1]:.9f}    {first: .9f}    '
          f'{second: .9f}    {gaussian(g):.9f}    {change:.1e}')

The matrix residual checks the finite eigenproblem. Basis enlargement checks what it misses. Neither test alone supplies a rigorous bound on the distance from the infinite-basis limit.

Reading the Results

Representative dimensionless ground energies are

$g$	Numerical energy	First order	Second order	Gaussian bound
$0.01$	$0.507256205$	$0.507500$	$0.5072375$	$0.507288$
$0.1$	$0.559146327$	$0.575000$	$0.548750$	$0.560307$
$1$	$0.803770651$	$1.250000$	$-1.375000$	$0.812500$

At small $g$ , the second-order estimate improves on first order. At $g=1$ it is already unphysical: adding a positive potential cannot lower the ground energy below $1/2$ . The minimized Gaussian stays an upper bound, although it does not reproduce the exact second-order coefficient.

Increasing $g$ makes the state narrower in position, not broader. A fixed-frequency basis may nevertheless need more functions to represent that narrower shape; an adjustable basis frequency can improve convergence. Variational Optimization explores that choice.

Exercises

Detect the projection error

Construct $q$ with only one oscillator state and compare $(P_1qP_1)^4$ with $\langle0|q^4|0\rangle$ .

Solution

The one-by-one position matrix is zero, so its fourth power is zero. The true moment is $3/4$ . The missing paths pass through excited oscillator states and return to the ground state. The code’s diagonal-moment assertion detects this error even at the upper edge of a larger basis.

Identify the omitted order

Repeat the calculation for $g=0.001,0.002,0.004$ . What happens to $[\epsilon_0^{\mathrm{num}}-\epsilon_0^{[2]}]/g^3$ as $g$ decreases?

Solution

It approaches the third-order coefficient $333/16=20.8125$ . Choose $g$ large enough that the $g^3$ difference remains above numerical roundoff, and reduce basis error before interpreting the trend. The full weak-coupling series is asymptotic; this fixed-order limit does not imply convergence as the order grows.

References

C. M. Bender and T. T. Wu, “Anharmonic oscillator,” Physical Review 184, 1231–1260 (1969). doi:10.1103/PhysRev.184.1231.
G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
R. Shankar, Principles of Quantum Mechanics, 2nd ed., Springer, 1994.
J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.