Borel-Padé Resummation

Perturbation Methods

Borel-Padé is the workhorse resummation technique for factorially divergent asymptotic series in quantum mechanics and field theory. The strategy combines two ideas. (1) BOREL: divide each coefficient by $n!$ to tame the factorial growth — the formal series of $a_{n} / n!$ has a finite radius of convergence even when the original series diverges everywhere. (2) PADÉ: rather than try to extend the Borel transform analytically by hand, fit a rational function to its first few coefficients. Together, the pair turns a wildly divergent perturbation expansion into a convergent numerical procedure whose accuracy improves order by order.

The setup. Consider an asymptotic series of the form

f (g) \sim n = 0 \sum \infty a_{n} g^{n}, a_{n} \sim C (- 1)^{n} n! / S^{n} for large n .

The factorial growth means the radius of convergence is zero — the partial sums diverge for every $g \neq = 0$ . This is the universal large-order behavior of perturbation theory in Hamiltonians with non-perturbative tunneling (the quartic anharmonic oscillator, the Stark effect, double-wells, gauge instantons in QFT). The constant $S$ in the denominator is precisely the action of the leading instanton — connecting perturbative coefficients to non-perturbative physics. See the resurgence page for that story.

The Borel transform

Define the Borel transform of the series by

B (t) = n = 0 \sum \infty \frac{a _{n}}{n !} t^{n} .

With $a_{n} / n! \sim C (- 1)^{n} / S^{n}$ , the coefficients of $B (t)$ are bounded — the series has finite radius of convergence $∣ t ∣ < S$ . The original function is recovered by the Laplace-type integral

f (g) = \int_{0}^{\infty} e^{- t} B (g t) d t .

The integral converges when $B (t)$ can be evaluated along the positive real axis and grows slower than $e^{t}$ there. Two conditions for this to make sense as a valid resummation: (i) the integral converges, (ii) the result is the same function the asymptotic series was approximating. Watson's theorem and its generalizations (Nevanlinna, Sokal) give the technical conditions; for Bender-Wu type series with definite-sign Stieltjes structure, both conditions are met.

But: we don't HAVE the Borel transform analytically. We have its first few coefficients $a_{0}, a_{1} /1!, a_{2} /2!, \dots, a_{N} / N!$ . Truncating $B (t)$ to a polynomial and integrating gives back exactly the original partial sum — no improvement. The polynomial truncation cannot extrapolate $B (t)$ for $t > ∣ S ∣$ or capture its singularity structure.

Padé on the Borel transform

The cure is to fit a RATIONAL function to the Borel coefficients instead of integrating their polynomial truncation. Build the $[m / n]$ Padé approximant of $B (t)$ , then integrate that rational function against $e^{- t}$ :

f_{m, n} (g) = \int_{0}^{\infty} e^{- t} [m / n]_{B} (g t) d t .

The rational approximant captures the leading SINGULARITIES of $B (t)$ as poles of its denominator. For the AHO this is structurally important: the Bender-Wu coefficients $a_{n} \sim (- 3)^{n} n!$ imply $B (t)$ has a singularity near $t = - 1/3$ (off the integration contour, fortunately for resummation), and Padé recovers that singularity from the first few coefficients. The integrand decays exponentially, so the integral is dominated by $g t ≲ 1/ g$ ; for not-too-large $g$ , the rational interpolation is excellent on the relevant range.

Code

# Borel-Pade resummation of a factorially divergent asymptotic series.
#
# Given coefficients a_n with a_n ~ (-1)^n n! / S^n for large n
# (Bender-Wu growth), the formal sum f(g) = sum_n a_n g^n is divergent
# for every g != 0. The Borel transform
#     B(t) = sum_n (a_n / n!) t^n
# has finite radius of convergence (= S), and the original function can
# be reconstructed via the Laplace integral
#     f(g) = int_0^infty exp(-t) B(g t) dt.
# Borel-Pade: replace B by a Pade approximant [m/n] of its truncated
# series, then integrate.

import numpy as np
from scipy.linalg import eigh
from scipy.integrate import quad
from math import factorial

def pade(coeffs, m, n):
    """[m/n] Pade approximant from coefficients a_0, ..., a_{m+n}."""
    M = m + n
    c = np.asarray(coeffs[:M+1], dtype=float)
    if n > 0:
        A = np.array([[c[m + i - k] if 0 <= m + i - k < len(c) else 0
                       for k in range(1, n + 1)]
                      for i in range(1, n + 1)])
        b_tail = np.linalg.solve(A, -c[m+1 : m+n+1])
        b = np.concatenate([[1.0], b_tail])
    else:
        b = np.array([1.0])
    a = np.zeros(m + 1)
    for i in range(m + 1):
        a[i] = sum(b[k] * c[i - k] for k in range(min(i, n) + 1))
    return a, b

def borel_pade_sum(coeffs, g, m, n):
    """Borel-Pade resummation of sum_k a_k g^k."""
    # Borel coefficients B_k = a_k / k!
    B = [coeffs[k] / factorial(k) for k in range(m + n + 1)]
    a_p, b_p = pade(B, m, n)
    def integrand(t):
        num = sum(a_p[i] * (g * t)**i for i in range(len(a_p)))
        den = sum(b_p[i] * (g * t)**i for i in range(len(b_p)))
        return np.exp(-t) * num / den
    val, _ = quad(integrand, 0, np.inf, limit=300,
                  epsabs=1e-14, epsrel=1e-14)
    return val

# ─── Apply to the AHO ground state ──────────────────────────────────────
# Convention: H = -1/2 d²/dx² + 1/2 x² + g x⁴.
# 'True' E_0(g) by SHO-basis diagonalization; coefficients by RS PT.

Nbasis = 80
X = np.zeros((Nbasis, Nbasis))
for n in range(Nbasis - 1):
    X[n, n+1] = np.sqrt((n + 1) / 2)
    X[n+1, n] = np.sqrt((n + 1) / 2)
X4 = X @ X @ X @ X
H0 = np.diag(np.arange(Nbasis) + 0.5)

def aho_exact(g):
    return float(eigh(H0 + g * X4, eigvals_only=True, subset_by_index=[0, 0])[0])

def bender_wu(Nmax):
    V = X4
    E = [0.5]
    psi = [np.zeros(Nbasis) for _ in range(Nmax + 1)]
    psi[0][0] = 1.0
    R = np.array([0.0] + [1.0 / n for n in range(1, Nbasis)])
    for k in range(1, Nmax + 1):
        Vpsi = V @ psi[k-1]
        Ek = psi[0] @ Vpsi
        E.append(Ek)
        rhs = -Vpsi + sum(E[j] * psi[k-j] for j in range(1, k))
        psi[k] = R * rhs; psi[k][0] = 0.0
    return E

A = bender_wu(11)

# ─── Compare partial sums, Pade, and Borel-Pade ─────────────────────────
def pade_eval(a, b, g):
    return (sum(a[i] * g**i for i in range(len(a))) /
            sum(b[i] * g**i for i in range(len(b))))

print(f"{'g':>5s}  {'true E_0':>14s}  {'[5/5] Pade err':>18s}  {'[5/5] B-Pade err':>20s}")
for g in [0.02, 0.05, 0.1, 0.2]:
    true = aho_exact(g)
    a_p, b_p = pade(A, 5, 5)
    pade_v   = pade_eval(a_p, b_p, g)
    bp_v     = borel_pade_sum(A, g, 5, 5)
    print(f"{g:5.2f}  {true:14.10f}  "
          f"{abs(pade_v - true):>18.2e}  "
          f"{abs(bp_v   - true):>20.2e}")

Output, comparing diagonal Padé and Borel-Padé on the AHO ground-state series:

    g        true E_0      [5/5] Pade err      [5/5] B-Pade err
 0.02    0.5140864273            2.18e-11              9.98e-14
 0.05    0.5326427548            2.71e-08              1.75e-10
 0.10    0.5591463272            2.12e-06              1.82e-08
 0.20    0.6024051637            6.92e-05              8.85e-07

The columns tell the story. At $g = 0.02$ , both Padé and Borel-Padé hit single-digit-of-precision accuracy or better, but Borel-Padé wins by three orders of magnitude ( $1 0^{- 13}$ vs $1 0^{- 11}$ ). At $g = 0.1$ , where the partial sum of the bare series has exploded to $\sim 6$ (vs true value $\sim 0.56$ ), Borel-Padé gives error $2 \times 1 0^{- 8}$ — eight digits correct from a series that nominally diverges everywhere. At $g = 0.2$ , Padé loses a couple of orders to Borel-Padé. The gap widens with $g$ : the harder the resummation problem, the more important it is to tame the factorial growth before fitting.

Why Borel-Padé works

Three structural reasons, all visible in the integrand.

(1) Factorial → exponential. The Borel transform converts $a_{n} \sim n!$ growth into the bounded coefficients $a_{n} / n! \sim 1$ . Rational fits to slowly varying coefficients are dramatically better than rational fits to factorially growing ones. The Padé linear system is well-conditioned.

(2) Singularity capture. The dominant non-analytic structure of $f (g)$ — the cut along the negative $g$ -axis encoding tunneling, the branch behavior, the resurgent singularities — gets packaged into a small number of SINGULARITIES OF $B (t)$ in the Borel plane. A low-order Padé denominator can model those singularities accurately. The Laplace integral then transports that singularity information back through $t \to g t$ , giving $f (g)$ nearly exactly.

(3) The integral is forgiving. The $e^{- t}$ kernel kills the contribution from large $t$ where the Padé approximant becomes unreliable. Effectively you only need the Padé to be accurate on the interval where $g t \sim 1$ , which is a SHORT interval at small $g$ and a longer one at large $g$ . Convergence in $g$ and the order of Padé approximation are linked through the support of the integrand.

When it fails

Borel-Padé is robust for series of the AHO type — alternating signs, factorial growth, Borel-summable. Two failure modes worth knowing:

Borel singularities on the positive real axis. If $B (t)$ has a pole or branch point at some $t = t_{+} > 0$ , the Laplace contour HITS the singularity. The integral becomes ambiguous — the contour can pass above or below the singularity, giving complex-conjugate values whose imaginary part is exponentially small in $g$ . This ambiguity is real physics: it is precisely the imaginary part of $f (g)$ from instanton-anti-instanton pairs. Resolving it requires the resurgent transseries machinery rather than naive Borel-Padé.
Branch points near zero. If $B (t)$ has a branch point near the origin of the Borel plane (some QFT series do), Padé struggles and you should switch to CONFORMAL MAPPING in the Borel plane — map the cut Borel plane to the unit disk and apply Padé in the mapped variable. The Le Guillou-Zinn-Justin technique used for high-precision φ⁴ critical exponents.

Practical recipe

For a Bender-Wu type asymptotic series in physics, here is the standard workflow.

Compute as many coefficients $a_{n}$ as you reasonably can. Each additional order makes the Padé more accurate. For the AHO, perturbation theory in the SHO basis gives them recursively; in QFT, diagrammatic computation; in operator-overlap methods, recursive matrix-element evaluation.
Form Borel coefficients $B_{n} = a_{n} / n!$ .
Build the diagonal $[m / m]$ Padé approximant of $B (t)$ for several $m$ . Check stability of the result across different $m$ — divergence between approximants is a red flag.
Integrate against $e^{- t}$ numerically (Gauss-Laguerre quadrature is the natural choice). Sanity-check by comparing to a finite-element / SHO-basis diagonalization at small $g$ .
If the result is unstable or the integral diverges, the series may not be Borel summable in the naïve sense — fall back to resurgent transseries or conformal mapping.

Padé approximants — the rational-fit half of the procedure.
Perturbation theory and resurgence — the broader story of what divergent series ENCODE and how to extract non-perturbative information from them.
Variational perturbation theory — alternative resummation that works at strong coupling where Borel-Padé loses accuracy.
WKB and semiclassical methods — the instanton action $S$ appearing in the Bender-Wu growth $a_{n} \sim n! / S^{n}$ is the same $S$ the WKB tunneling exponential.