WKB Versus Exact Spectrum Notebook

This calculation compares the leading WKB spectrum of a quadratic-plus-quartic well with independently converged numerical eigenvalues. A harmonic well checks the phase convention; a quartic contribution reveals the approximation error. The complete Python example requires only NumPy.

Required background. Bohr–Sommerfeld quantization supplies the action condition, and turning-point matching supplies its half-unit shift. The pure quartic worked example evaluates the corresponding action analytically.

Action Quantization in a Smooth Well

Use dimensionless oscillator units and

h=-\frac12\frac{d^2}{dq^2}+V(q),\qquad V(q)=\frac12q^2+gq^4,\qquad g\ge0.

For energy $E>0$ , let $a$ be the positive turning point. A cancellation-free expression is

a^2=\frac{4E}{1+\sqrt{1+16gE}}.

It gives $a^2=2E$ at $g=0$ . The half-orbit action must satisfy

I(E)=2\int_0^a\sqrt{2[E-V(q)]}\,dq=\pi(n+1/2).

The action is monotone in $E$ , so each level can be found by bisection once it is bracketed. The leading rule is exact for harmonic energies, but it is generally approximate for $g>0$ .

Smooth the Turning-Point Integral

Direct quadrature encounters a square-root endpoint. Set $q=a\sin\theta$ and use

E-V(a\sin\theta)=a^2\cos^2\theta \left[\frac12+ga^2(1+\sin^2\theta)\right].

The integrand becomes smooth on $0\le\theta\le\pi/2$ :

I(E)=2a^2\int_0^{\pi/2}\cos^2\theta \sqrt{1+2ga^2(1+\sin^2\theta)}\,d\theta.

At $g=0$ , this gives $I(E)=\pi E$ directly.

Run the Comparison

The numerical reference projects the original Hamiltonian onto oscillator states. Padding before multiplying $q$ preserves the quartic matrix elements; see Perturbation Theory Benchmarks for the reason. The code doubles both the quadrature order and basis size independently.

import sys
import numpy as np
from numpy.polynomial.legendre import leggauss

def spectrum(N, g):
    a = np.diag(np.sqrt(np.arange(1, N + 4)), 1)
    q = (a + a.T) / np.sqrt(2)
    H = np.diag(np.arange(N) + 0.5) + g*np.linalg.matrix_power(q, 4)[:N, :N]
    return np.linalg.eigvalsh(H)[:8]

def wkb_levels(g, order):
    x, weight = leggauss(order)
    theta, weight = (x + 1)*np.pi/4, weight*np.pi/4
    def action(E):
        a2 = 4*E / (1 + np.sqrt(1 + 16*g*E))
        f = np.cos(theta)**2 * np.sqrt(1 + 2*g*a2*(1 + np.sin(theta)**2))
        return 2*a2*np.dot(weight, f)
    result = []
    for n in range(8):
        target = np.pi*(n + 0.5)
        lo, hi = 0.0, n + 1.0
        while action(hi) < target:
            hi *= 2
        for _ in range(70):
            mid = (lo + hi)/2
            if action(mid) < target:
                lo = mid
            else:
                hi = mid
        result.append((lo + hi)/2)
    return np.array(result)

print('Python', sys.version.split()[0], 'NumPy', np.__version__)
for g in (0, 0.01, 0.1, 1):
    exact_small, numerical = spectrum(160, g), spectrum(320, g)
    rough, wkb = wkb_levels(g, 32), wkb_levels(g, 64)
    assert np.max(abs(numerical - exact_small)) < 1e-8
    assert np.max(abs(wkb - rough)) < 1e-10
    if g == 0:
        assert np.allclose(wkb, np.arange(8) + 0.5, atol=1e-12, rtol=0)
    print('\ng =', g)
    print('n   numerical      WKB            signed error (%)')
    for n, (E, W) in enumerate(zip(numerical, wkb)):
        print(f'{n}   {E:.9f}    {W:.9f}    {100*(W-E)/E:+.6f}')

Here “numerical” means a converged approximation to the differential eigenproblem. It is independent of the action rule, but finite-basis stability remains a diagnostic rather than an exact error certificate.

Interpret the Difference

For $g=0$ , disagreement would indicate a numerical or normalization error. For positive $g$ , a stable difference measures the limitation of leading WKB. Its relative error generally falls as $n$ grows because the orbit contains more phase oscillations. Keep the sign: WKB is not a variational bound, even if every displayed error for this family has the same sign.

Do not enlarge the level range without rechecking the basis. The highest retained numerical states are the least reliable. For a separate real-space check, solve on $[-L,L]$ with Dirichlet boundaries and vary both $L$ and the grid spacing; changing only the grid cannot reveal a box that is too small.

Exercises

Two different numerical errors

The WKB energies stop changing when the quadrature order doubles, but the numerical spectrum changes when the basis grows. Which comparison is ready to interpret?

Solution

Only the action integral has converged. The reference spectrum still has unresolved basis error, so its difference from WKB cannot yet be assigned entirely to semiclassics. Increase the basis or choose a better basis scale, then repeat both checks.

Recover the quartic scaling

When the quartic term dominates, why should $E_n$ grow as $g^{1/3}(n+1/2)^{4/3}$ ?

Solution

The turning point scales as $(E/g)^{1/4}$ , while momentum scales as $E^{1/2}$ . Thus $I(E)\propto g^{-1/4}E^{3/4}$ . Setting this proportional to $n+1/2$ gives the stated powers. The pure quartic calculation determines the coefficient and its error.

References

C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory, Springer, 1999.
M. V. Berry and K. E. Mount, “Semiclassical approximations in wave mechanics,” Reports on Progress in Physics 35, 315–397 (1972). doi:10.1088/0034-4885/35/1/306.
L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory, 3rd ed., Pergamon, 1977.
J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.