Higher-Order WKB Mathematica Code

• Higher-order WKB Mathematica code is a set of implementations for computing quasinormal modes and greybody factors in black-hole perturbation theory.
• It employs methods like Padé resummation, local Taylor expansion at the potential peak, and Bender–Wu recursion to generate higher-order corrections.
• Optimized strategies improve performance and scalability by efficiently handling derivatives, enhancing precision in semiclassical and scattering calculations.

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 (Konoplya et al., 12 Mar 2026). The modern literature includes a public 13th-order package with Padé resummation (Konoplya et al., 2019), an updated and optimized version extended to 16th order (Konoplya et al., 12 Mar 2026), and a distinct very-high-order automatic implementation based on Bender–Wu recursion, Padé approximants, and Cesàro means (Konoplya et al., 25 May 2026). Related Mathematica infrastructures also include the general perturbative package BenderWu (Sulejmanpasic et al., 2016).

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 (Konoplya et al., 12 Mar 2026). 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}$ (Konoplya et al., 25 May 2026).

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 (Konoplya et al., 12 Mar 2026). The second is the scattering data, namely the reflection and transmission coefficients $x$0, from which the greybody factor

$x$1

is obtained (Konoplya et al., 12 Mar 2026).

The method is organized around the potential maximum. The WKB formulas depend on

$x$2

and on derivatives at the peak,

$x$3

or, in the $x$4 convention, on $x$5 evaluated at $x$6 (Konoplya et al., 12 Mar 2026, Konoplya et al., 25 May 2026). A recurrent practical point is that an $x$7-th order WKB approximation requires derivatives of the effective potential up to order $x$8 (Konoplya et al., 25 May 2026).

2. Quantization, scattering, and resummation structure

The core higher-order WKB condition implemented in the optimized black-hole package is

$x$9

where $dx = \frac{dr}{f(r)},$0 are higher-order WKB corrections built from derivatives $dx = \frac{dr}{f(r)},$1 (Konoplya et al., 12 Mar 2026). In the alternative notation of the very-high-order implementation, the standard rule is written as

$dx = \frac{dr}{f(r)},$2

which yields the explicit eigenvalue relation

$dx = \frac{dr}{f(r)},$3

(Konoplya et al., 25 May 2026).

At the eikonal level, the WKB parameter is

$dx = \frac{dr}{f(r)},$4

For QNMs, $dx = \frac{dr}{f(r)},$5 is quantized as

$dx = \frac{dr}{f(r)},$6

so the code solves the WKB equation with this quantized value (Konoplya et al., 12 Mar 2026).

For scattering with real potentials, the transmission probability is given by

$dx = \frac{dr}{f(r)},$7

while

$dx = \frac{dr}{f(r)},$8

(Konoplya et al., 12 Mar 2026). 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 (Konoplya et al., 2019, Konoplya et al., 12 Mar 2026).

A persistent feature of the literature is Padé resummation. The WKB series is treated as a formal expansion in a bookkeeping parameter $dx = \frac{dr}{f(r)},$9,

$U(x,\omega)$0

or, in the $U(x,\omega)$1 notation,

$U(x,\omega)$2

(Konoplya et al., 12 Mar 2026, Konoplya et al., 25 May 2026). The polynomial is then replaced by a rational approximant,

$U(x,\omega)$3

or

$U(x,\omega)$4

with $U(x,\omega)$5 (Konoplya et al., 12 Mar 2026, Konoplya et al., 25 May 2026). 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 (Konoplya et al., 12 Mar 2026, Konoplya et al., 25 May 2026).

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 (Konoplya et al., 2019). That framework implemented corrections up to 13th order and became widely used in black-hole spectroscopy and scattering problems (Konoplya et al., 12 Mar 2026).

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 (Konoplya et al., 12 Mar 2026). This updated code extends the earlier public implementation from 13th order to 16th order, using higher-order formulas derived in Matyjasek et al. (Konoplya et al., 12 Mar 2026).

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 (Konoplya et al., 25 May 2026). 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 (Konoplya et al., 25 May 2026).

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 (Sulejmanpasic et al., 2016). 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 (Konoplya et al., 25 May 2026). 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 (Konoplya et al., 12 Mar 2026). In the direct QNM workflow, the input consists of the black-hole metric function $U(x,\omega)$6, the effective potential $U(x,\omega)$7 or $U(x,\omega)$8, spin or field parameters, multipole number $U(x,\omega)$9, overtone number $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$0, and the WKB order (Konoplya et al., 12 Mar 2026). The workflow is then:

1. locate the potential maximum $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$1,
2. compute $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$2, $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$3, and higher derivatives at $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$4,
3. assemble the WKB equation,
4. solve for $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$5 from the QNM quantization condition

(Konoplya et al., 12 Mar 2026). The output is the complex quasinormal frequency, optionally improved by Padé resummation (Konoplya et al., 12 Mar 2026).

The greybody-factor workflow uses the same potential data and a real frequency range. It computes $$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$6 from the WKB condition and then obtains

$$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$7

(Konoplya et al., 12 Mar 2026). The code can also output the reflection coefficient if needed (Konoplya et al., 12 Mar 2026).

The analytic large-$$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$8 workflow expands in inverse powers of

$$\frac{d^2\Psi}{dr_*^2}+Q(r_*)\Psi=0,\qquad Q(r_*)=\omega^2-V(r_*),$$9

The peak position is written as

$r_*$0

and the WKB expression is expanded consistently in $r_*$1, producing analytic approximations for QNMs and greybody factors (Konoplya et al., 12 Mar 2026).

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 (Konoplya et al., 25 May 2026). 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 $r_*$2,
5. use Padé or Wynn-based routines,
6. optionally apply CesaroMean,
7. compare to accurate reference frequencies obtained via the Leaver continued-fraction method

(Konoplya et al., 25 May 2026).

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

$r_*$3

reduces perturbations to

$r_*$4

constructs the scalar potential

$r_*$5

finds the barrier maximum, evaluates derivatives with respect to $r_*$6, inserts the standard 13th-order $r_*$7 corrections, and solves for $r_*$8 (Avalos et al., 2023). 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 (Konoplya et al., 12 Mar 2026). The earlier package computed higher derivatives symbolically by differentiating the full analytic expression for $r_*$9 with respect to the tortoise coordinate and then substituting the maximum (Konoplya et al., 12 Mar 2026). This becomes expensive when the potential contains complicated metric functions, nested radicals, exponentials, non-rational terms, or generally messy analytic structures (Konoplya et al., 12 Mar 2026).

The updated implementation instead finds the maximum numerically, expands the potential locally as

$V$0

with $V$1 at the maximum, and evaluates the derivatives numerically there (Konoplya et al., 12 Mar 2026). 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 (Konoplya et al., 12 Mar 2026). 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 (Konoplya et al., 12 Mar 2026).

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 (Konoplya et al., 12 Mar 2026). 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$2 for complicated potentials (Konoplya et al., 12 Mar 2026).

The very-high-order Bender–Wu implementation pursues a different form of scalability. Rather than optimize a fixed set of $V$3, it automatically generates coefficients to very high order and is described as limited only by memory and computational time (Konoplya et al., 25 May 2026). 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 (Sulejmanpasic et al., 2016). 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 (Konoplya et al., 2019, Konoplya et al., 25 May 2026). It works best when $V$4, or more conservatively for low overtones and dominant modes where the barrier approximation is clean (Konoplya et al., 2019). The optimized 2026 package likewise emphasizes black-hole scattering and ringing problems in which only the behavior near the potential maximum matters (Konoplya et al., 12 Mar 2026).

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 (Konoplya et al., 2019). 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 $V$5, and the possibility that fourth-order WKB is not uniformly better than third order in some higher-derivative-gravity modes (Tang et al., 19 Dec 2025).

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 (Konoplya et al., 25 May 2026). A specific example reported in the paper shows that for $V$6, even 800th-order WKB does not converge to the correct frequency despite apparent stabilization (Konoplya et al., 25 May 2026). 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 (Konoplya et al., 25 May 2026).

Cesàro averaging is introduced precisely as an internal convergence diagnostic. The running mean

$V$7

and especially the Cesàro mean

$V$8

often become monotonically convergent once sufficiently high order is reached, even when the raw Padé sequence oscillates non-monotonically (Konoplya et al., 25 May 2026). 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 (Konoplya et al., 25 May 2026).

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 (Konoplya et al., 25 May 2026).

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 (Sulejmanpasic et al., 2016); 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 (Kreshchuk et al., 2018). In that setting, the second-order correction takes the universal form

$V$9

and is then reduced to derivatives of the classical action by symbolic basis reduction (Kreshchuk et al., 2018). A further development computes quantum periods for higher-genus hyperelliptic curves through Picard–Fuchs operators and quantum differential operators $V_0$0 acting on classical periods (Fischbach et al., 2018).

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

$V_0$1

into a smoother first-order system for a transformed variable $V_0$2, then apply Picard-based one-step schemes of order $V_0$3 built from asymptotic quadratures of oscillatory integrals (Arnold et al., 2023, Arnold et al., 2024). 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 (Konoplya et al., 2019, Konoplya et al., 12 Mar 2026, Konoplya et al., 25 May 2026). 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 (Sulejmanpasic et al., 2016, Kreshchuk et al., 2018, Fischbach et al., 2018, Arnold et al., 2023, Arnold et al., 2024). 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.