---
title: Higher-Order WKB Mathematica Code
url: https://www.emergentmind.com/topics/higher-order-wkb-mathematica-code
type: topic
---

# Higher-Order WKB Mathematica Code

Searching arXiv for the cited higher-order WKB Mathematica papers and closely related implementations.
Higher-order WKB Mathematica code denotes a family of Mathematica implementations of the higher-order Wentzel–Kramers–Brillouin approximation used to compute quasinormal modes (QNMs), greybody factors (GBFs), and related semiclassical quantities from one-dimensional barrier problems, especially in black-hole perturbation theory. In the black-hole setting, these codes solve the standard wave equation in tortoise coordinate, extract the effective potential near its maximum, assemble higher-order WKB corrections, and determine either the complex quasinormal frequency or the reflection and transmission amplitudes [2603.12466]. The modern literature includes a public 13th-order package with Padé resummation [1904.10333], an updated and optimized version extended to 16th order [2603.12466], and a distinct very-high-order automatic implementation based on Bender–Wu recursion, Padé approximants, and Cesàro means [2605.25705]. Related Mathematica infrastructures also include the general perturbative package **BenderWu** [1608.08256].

## 1. Scope and mathematical setting

The principal target of higher-order WKB Mathematica code is the one-dimensional wave equation for perturbations,
\[
\frac{d^2\Psi}{dx^2} = U(x,\omega)\,\Psi,
\]
where \(x\) is the tortoise coordinate defined by
\[
dx = \frac{dr}{f(r)},
\]
and \(U(x,\omega)\) is a single-barrier effective potential [2603.12466]. An equivalent form used in black-hole perturbation theory is
\[
\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),
\]
with \(r_*\) the tortoise coordinate and \(V\) a potential with a single peak \(V_0\) at \(r_{*0}\) [2605.25705].

From this barrier data, the code computes two classes of quantities. The first is the set of quasinormal frequencies \(\omega\), defined by outgoing or ingoing asymptotic boundary conditions at infinity and the horizon [2603.12466]. The second is the scattering data, namely the reflection and transmission coefficients \(R,T\), from which the greybody factor
\[
\Gamma(\omega)=|T|^2
\]
is obtained [2603.12466].

The method is organized around the potential maximum. The WKB formulas depend on
\[
U_0(\omega)\equiv U(x_m,\omega),
\]
and on derivatives at the peak,
\[
U_1=0,\qquad U_2=\left.\frac{d^2U}{dx^2}\right|_{x_m},\qquad U_3=\left.\frac{d^3U}{dx^3}\right|_{x_m},\ \ldots
\]
or, in the \(Q=\omega^2-V\) convention, on \(Q_0,Q_2,\dots,Q_{2N}\) evaluated at \(r_{*0}\) [2603.12466; 2605.25705]. A recurrent practical point is that an \(N\)-th order WKB approximation requires derivatives of the effective potential up to order \(2N\) [2605.25705].

## 2. Quantization, scattering, and resummation structure

The core higher-order WKB condition implemented in the optimized black-hole package is
\[
0 = U_0(\omega) + A_2(\mathcal{K}^2) + A_4(\mathcal{K}^2) + A_6(\mathcal{K}^2) + \cdots
 - i \mathcal{K}\sqrt{-2U_2(\omega)} \left(1 + A_3(\mathcal{K}^2) + A_5(\mathcal{K}^2) + A_7(\mathcal{K}^2)+\cdots\right),
\]
where \(A_k(\mathcal{K}^2)\) are higher-order WKB corrections built from derivatives \(U_2,\dots,U_{2k}\) [2603.12466]. In the alternative notation of the very-high-order implementation, the standard rule is written as
\[
\frac{iQ_0}{\sqrt{2Q_0''}-\sum_{j=2}^{N}\Lambda_j=\mathcal K,\qquad \mathcal K\equiv n+\frac{1}{2},
\]
which yields the explicit eigenvalue relation
\[
\omega^2=V_0-i\mathcal K\sqrt{-2V_0''}-i\sqrt{-2V_0''}\sum_{j=2}^{N}\Lambda_j
\]
[2605.25705].

At the eikonal level, the WKB parameter is
\[
\mathcal{K}=-i\,\frac{U_0(\omega)}{\sqrt{-2U_2(\omega)}}.
\]
For QNMs, \(\mathcal{K}\) is quantized as
\[
\mathcal{K}= \begin{cases} +n+\tfrac12, & \Re(\omega)>0,\\[4pt] -n-\tfrac12, & \Re(\omega)<0, \end{cases}
\qquad n=0,1,2,\ldots
\]
so the code solves the WKB equation with this quantized value [2603.12466].

For scattering with real potentials, the transmission probability is given by
\[
\Gamma(\omega)=|T|^2=\frac{1}{1+e^{2\pi i\mathcal{K}}},
\]
while
\[
|R|^2=\frac{1}{1+e^{-2\pi i\mathcal{K}}},\qquad |T|^2 = \frac{1}{1+e^{2\pi i\mathcal{K}}}=1-|R|^2
\]
[2603.12466]. The same barrier-top quantity that determines QNM poles also controls the greybody factor, which is why the same code base can serve both ringing and scattering computations [1904.10333; 2603.12466].

A persistent feature of the literature is Padé resummation. The WKB series is treated as a formal expansion in a bookkeeping parameter \(\epsilon\),
\[
\mathcal{K}(\epsilon)=\sum_{n=0}^{N} a_n \epsilon^n
\]
or, in the \(P_N(\lambda)\) notation,
\[
P_N(\lambda)=V_0-i\mathcal K\sqrt{-2V_0''}\,\lambda -i\sqrt{-2V_0''}\sum_{j=2}^{N}\lambda^j\Lambda_j,\qquad \omega^2\approx P_N(1)
\]
[2603.12466; 2605.25705]. The polynomial is then replaced by a rational approximant,
\[
P_{m/n}(\epsilon)=\frac{\sum_{i=0}^{m} p_i \epsilon^i}{1+\sum_{j=1}^{n} q_j \epsilon^j},
\]
or
\[
\mathcal P^{\tilde m}_{\tilde n}(\lambda)= \frac{\sum_{j=0}^{\tilde m} A_j\lambda^j}{1+\sum_{j=1}^{\tilde n} B_j\lambda^j},
\qquad \tilde m+\tilde n=N,
\]
with \(\omega^2\approx \mathcal P^{\tilde m}_{\tilde n}(1)\) [2603.12466; 2605.25705]. The updated package emphasizes balanced Padé approximants, while the very-high-order code highlights diagonal and near-diagonal Padé approximants as the most useful choices [2603.12466; 2605.25705].

## 3. Evolution of Mathematica implementations

The public code history reflected in the literature begins with a 13th-order WKB package and Padé resummation for QNMs and greybody factors [1904.10333]. That framework implemented corrections up to 13th order and became widely used in black-hole spectroscopy and scattering problems [2603.12466].

A later optimized version retained the same method but changed the internal strategy for obtaining peak derivatives. Its principal improvement is to expand the effective potential in a Taylor series around its maximum rather than evaluate the full analytic expression of the WKB formula for each specific potential [2603.12466]. This updated code extends the earlier public implementation from 13th order to 16th order, using higher-order formulas derived in Matyjasek et al. [2603.12466].

A separate line of development replaced fixed-order hand-coded correction formulas with automatic generation at very high order. The fully automatic Mathematica implementation described in 2026 is built around the Bender–Wu algorithm, produces WKB coefficients recursively to very high order, resums them with Padé or Wynn-type procedures, and monitors convergence through Cesàro means [2605.25705]. This approach is intended to go well beyond the usual 3rd- or 6th-order treatments and is described as fully automatic and, in principle, limited only by memory and CPU time [2605.25705].

The Bender–Wu-based black-hole code rests on the more general **BenderWu** package for locally harmonic one-dimensional quantum-mechanical potentials. That package automates exact recursion relations for perturbative wave-function coefficients and energy corrections to arbitrarily high order, and exposes three principal functions: `BenderWu`, `BWProcess`, and `BWLevelPolynomial` [1608.08256]. In the black-hole WKB setting, the local barrier-top problem is reformulated as a perturbative anharmonic-oscillator problem, which allows the same recursion machinery to be used [2605.25705]. This suggests a close methodological continuity between high-order perturbation theory and very-high-order WKB generation.

## 4. Algorithmic organization and Mathematica workflows

The optimized black-hole package is organized around three workflows: direct WKB/QNM computation, greybody-factor computation, and analytic expansion beyond eikonal order [2603.12466]. In the direct QNM workflow, the input consists of the black-hole metric function \(f(r)\), the effective potential \(V(r)\) or \(U(x,\omega)\), spin or field parameters, multipole number \(\ell\), overtone number \(n\), and the WKB order [2603.12466]. The workflow is then:

1. locate the potential maximum \(x_m\),
2. compute \(U_0\), \(U_2\), and higher derivatives at \(x_m\),
3. assemble the WKB equation,
4. solve for \(\omega\) from the QNM quantization condition

[2603.12466]. The output is the complex quasinormal frequency, optionally improved by Padé resummation [2603.12466].

The greybody-factor workflow uses the same potential data and a real frequency range. It computes \(\mathcal K\) from the WKB condition and then obtains
\[
\Gamma(\omega)=\frac{1}{1+e^{2\pi i\mathcal K}}
\]
[2603.12466]. The code can also output the reflection coefficient if needed [2603.12466].

The analytic large-\(\ell\) workflow expands in inverse powers of
\[
\kappa\equiv \ell+\tfrac12.
\]
The peak position is written as
\[
x_m=x_0+x_1\kappa^{-1}+x_2\kappa^{-2}+\cdots+x_k\kappa^{-k}+\mathcal O(\kappa^{-k-1}),
\]
and the WKB expression is expanded consistently in \(\kappa^{-1}\), producing analytic approximations for QNMs and greybody factors [2603.12466].

The very-high-order automatic implementation has a more explicit Mathematica routine structure. It uses `PrecisePotentialPeak` to locate the relevant maximum of the potential, `WKBOrdersList` to pass the peak data to `BenderWu` and generate the correction list, `WKBOrdersPadefication` to construct Padé approximants, Wynn-based routines to generate Padé-type sequences, and `CesaroMean` to compute Cesàro averages [2605.25705]. Its practical workflow is:

1. specify the metric and field potential,
2. use `PrecisePotentialPeak`,
3. feed the peak and derivatives into `WKBOrdersList`,
4. obtain high-order WKB coefficients and form \(P_N(\lambda)\),
5. use Padé or Wynn-based routines,
6. optionally apply `CesaroMean`,
7. compare to accurate reference frequencies obtained via the Leaver continued-fraction method

[2605.25705].

A distinct but related Mathematica coding pattern appears in applications that implement fixed-order WKB directly from a chosen metric function. For example, the 13th-order treatment of hairy black holes starts from the metric
\[
ds^2=f(r)dt^2-\frac{dr^2}{f(r)}-r^2(d\theta^2+\sin^2\theta\,d\phi^2),
\]
reduces perturbations to
\[
\left(\frac{d^2}{dr_*^2}+\omega^2-V(r_*)\right)\chi(r_*)=0,
\qquad \frac{dr_*}{dr}=\frac1{f(r)},
\]
constructs the scalar potential
\[
V_L(r)=f(r)\left(\frac{L(L+1)}{r^2}+\frac{f'(r)}{r}\right),
\]
finds the barrier maximum, evaluates derivatives with respect to \(r_*\), inserts the standard 13th-order \(\Lambda_j\) corrections, and solves for \(\omega\) [2302.09148]. This application-specific pattern is representative of how the more general public packages are typically used.

## 5. Optimization strategies and performance characteristics

The defining optimization in the 2026 updated package is the replacement of symbolic differentiation of the full potential by local numerical Taylor expansion around the peak [2603.12466]. The earlier package computed higher derivatives symbolically by differentiating the full analytic expression for \(U(x,\omega)\) with respect to the tortoise coordinate and then substituting the maximum [2603.12466]. This becomes expensive when the potential contains complicated metric functions, nested radicals, exponentials, non-rational terms, or generally messy analytic structures [2603.12466].

The updated implementation instead finds the maximum numerically, expands the potential locally as
\[
U(x,\omega)=U_0+\frac{U_2}{2!}(x-x_m)^2+\frac{U_3}{3!}(x-x_m)^3+\cdots,
\]
with \(U_1=0\) at the maximum, and evaluates the derivatives numerically there [2603.12466]. Because higher-order WKB formulas depend only on derivatives at the peak and not on the full global form of the potential, the symbolic algebra bottleneck is bypassed [2603.12466]. The reported consequence is that, for complicated models, run times can be reduced from minutes or hours to fractions of a second while preserving the accuracy of the method [2603.12466].

Benchmarking in that work uses a Dirac field in Schwarzschild–de Sitter and a scalar field in a regular black-hole background from proper-time flow in quantum gravity [2603.12466]. The results show much smaller runtimes, less floating-point precision loss, especially dramatic improvement at higher WKB orders, and near-insensitivity to the analytic complexity of \(V(r)\) for complicated potentials [2603.12466].

The very-high-order Bender–Wu implementation pursues a different form of scalability. Rather than optimize a fixed set of \(\Lambda_j\), it automatically generates coefficients to very high order and is described as limited only by memory and computational time [2605.25705]. In the underlying general-purpose BenderWu package, the published performance claims are about 100 orders in 10–30 seconds and 250 orders in 1–2 hours on a modern home computer, with symbolic computation by default and optional numerical evaluation through `Evaluation -> "Numerical"` and `WorkingPrecision` [1608.08256]. A plausible implication is that the black-hole adaptation inherits the same emphasis on recursive generation, symbolic structure, and controllable resource use.

## 6. Domains of applicability, diagnostics, and failure modes

The standard black-hole higher-order WKB method assumes a single barrier with two turning points [1904.10333; 2605.25705]. It works best when \(\ell\gg n\), or more conservatively for low overtones and dominant modes where the barrier approximation is clean [1904.10333]. The optimized 2026 package likewise emphasizes black-hole scattering and ringing problems in which only the behavior near the potential maximum matters [2603.12466].

The literature repeatedly warns against extending the method beyond this regime. The 2019 review lists superradiance, instability, quasiresonances, high overtones, asymptotically non-constant potentials, and multiple turning points as situations in which the standard WKB barrier matching can become unreliable or misleading [1904.10333]. The rotating-black-hole extension similarly stresses the potential barrier requirement, the deterioration for higher overtones, the complexity introduced by complex-valued Kerr potentials, self-consistent dependence of angular separation constants on \(\omega\), and the possibility that fourth-order WKB is not uniformly better than third order in some higher-derivative-gravity modes [2512.17786].

The very-high-order automatic code refines this caution with the notion of **apparent convergence**. In strongly nonmoderate black-hole metrics, especially when higher near-horizon coefficients become large, the Padé sequence may stabilize numerically around a value that is nevertheless far from the accurate quasinormal frequency [2605.25705]. A specific example reported in the paper shows that for \(a_2=10\), even 800th-order WKB does not converge to the correct frequency despite apparent stabilization [2605.25705]. The practical guidance is therefore not to trust mere stabilization at modest order, to inspect diagonal Padé subsequences, compare with Cesàro means, and benchmark against an independent method such as Leaver or Frobenius when possible [2605.25705].

Cesàro averaging is introduced precisely as an internal convergence diagnostic. The running mean
\[
\overline{\omega}_{N,\Delta}=\frac{1}{\Delta+1}\sum_{s=0}^{\Delta}\omega_{N-s}
\]
and especially the Cesàro mean
\[
\sigma_N=\frac{1}{N+1}\sum_{k=0}^{N}\omega_k
\]
often become monotonically convergent once sufficiently high order is reached, even when the raw Padé sequence oscillates non-monotonically [2605.25705]. The paper stresses, however, that Cesàro means usually converge more slowly than the Padé sequence itself and should be used as a diagnostic rather than automatically as the best estimator [2605.25705].

Another subtlety concerns massive-field potentials with multiple turning points. The 2026 paper notes that massive fields may develop three turning points or even lose the local maximum. As long as the potential retains a local maximum, high-order Padé-resummed WKB often still works well, but this is presented as an empirical success rather than a strict derivation for the three-turning-point case [2605.25705].

## 7. Related Mathematica traditions beyond black-hole barrier codes

Although the phrase “higher-order WKB Mathematica code” most commonly refers to black-hole QNM and greybody-factor packages, Mathematica implementations of higher-order WKB also appear in broader semiclassical and quantum-mechanical contexts. One line uses the Bender–Wu recursion for perturbation theory around a locally harmonic point [1608.08256]; another uses Picard–Fuchs equations to express higher WKB corrections as linear combinations of classical actions and their energy derivatives for sextic and Lamé potentials [1803.07566]. In that setting, the second-order correction takes the universal form
\[
\mathcal S_2(E)=\oint_{\mathcal C_R}\frac{-mV''(x)}{24[2m(E-V(x))]^{3/2}}\,dx
\]
and is then reduced to derivatives of the classical action by symbolic basis reduction [1803.07566]. A further development computes quantum periods for higher-genus hyperelliptic curves through Picard–Fuchs operators and quantum differential operators \(\mathcal D_{2n}\) acting on classical periods [1803.11222].

A different usage appears in numerical solvers for the highly oscillatory one-dimensional stationary Schrödinger equation. There, Mathematica-oriented WKB codes implement analytic preprocessing that transforms
\[
\varepsilon^2\varphi''(x)+a(x)\varphi(x)=0
\]
into a smoother first-order system for a transformed variable \(Z\), then apply Picard-based one-step schemes of order \(\mathcal O(h^3)\) built from asymptotic quadratures of oscillatory integrals [2310.00963; 2402.18406]. These are not black-hole QNM packages, but they share the same core idea: remove the dominant oscillation analytically and compute higher-order corrections from local structure.

This broader landscape clarifies that higher-order WKB Mathematica code is not a single package but a technical ecosystem. In black-hole spectroscopy, it refers mainly to automated barrier-top codes for QNMs and greybody factors [1904.10333; 2603.12466; 2605.25705]. In semiclassical quantum mechanics, it also denotes symbolic or asymptotic frameworks that compute quantum periods, perturbative coefficients, or oscillatory Schrödinger solutions at high order [1608.08256; 1803.07566; 1803.11222; 2310.00963; 2402.18406]. The common denominator is the conversion of local analytic information—usually derivatives at a turning point, maximum, or harmonic extremum—into systematically improvable WKB data within Mathematica.

Source: https://www.emergentmind.com/topics/higher-order-wkb-mathematica-code