Third-Order Phase-Integral Method

• The third-order phase-integral method is an asymptotic approach that removes dominant oscillations via a second-order WKB phase, enabling efficient computation in non-turning regimes.
• It employs a one-step third-order accurate discretization using Picard expansions and specialized quadratures to approximate nested oscillatory integrals with controlled error bounds.
• The method is versatile, applying to both numerical Schrödinger problems and weighted stationary-phase evaluations by incorporating explicit boundary corrections and uniform remainder estimates.

Searching arXiv for the cited papers and closely related foundational work. arXiv search query: "WKB-based third order method for the highly oscillatory 1D stationary Schrödinger equation" arXiv search query: "Weighted stationary phase of higher orders" The third-order phase-integral method denotes a class of asymptotic and numerical constructions that retain phase-integral, WKB, or stationary-phase information through third order. In the semiclassical ordinary-differential-equation setting, it is a WKB-based one-step scheme for the highly oscillatory $1$D stationary Schrödinger equation that analytically removes the dominant oscillations and then integrates a smoother first-order system with third-order accuracy in the step size (Arnold et al., 2024). In the oscillatory-integral setting, it is a third-order weighted stationary-phase expansion on a finite interval with a single non-degenerate stationary point, explicit boundary terms, and a quantified remainder (McKee et al., 2016). The two usages are closely related through the classical phase-integral/JWKB idea that oscillatory structure should be encoded analytically in the phase before one approximates amplitudes or residual terms.

1. Semiclassical equation and classical phase-integral ansatz

In the numerical Schrödinger setting, the governing model is the $1$D stationary Schrödinger equation in semiclassical scaling

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$

with $0<\varepsilon\ll 1$. In quantum mechanical applications one has $a(x)=E-V(x)$, so that the standard physical form

$-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$

is recovered by identifying

$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$

For small $\varepsilon$, solutions oscillate with local wavelength

$\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$

The regime treated in the third-order WKB-based method is specified by Hypothesis A: $a\in C^7(I)$, $1$0 real-valued, and $1$1. This is a non-turning, oscillatory regime; turning points with $1$2 are excluded (Arnold et al., 2024).

The initial-value formulation is

$1$3

which matches semiclassical scaling. The associated classical phase-integral/WKB ansatz is

$1$4

Matching powers of $1$5 yields the first terms

$1$6

$1$7

and

$1$8

Hence the second-order WKB approximation used for normalization is

$1$9

with corrected phase

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$0

This choice already incorporates the $$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$1 phase correction associated with the standard Langer/Fröman correction term $$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$2.

2. Removal of fast oscillations and reduction to a smooth first-order system

The numerical third-order phase-integral method is based on an analytic change of variables that removes the leading oscillation before discretization. The second-order scalar equation is first written as a $$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$3 first-order system by defining

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$4

Then

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$5

where

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$6

The fast oscillations are then diagonalized by the unitary transformation

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$7

with

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$8

The transformed problem becomes

$$\varepsilon^2\,\varphi''(x)+a(x)\,\varphi(x)=0,\qquad x\in I=[x_0,x_{\mathrm{end}}],$$9

where $0<\varepsilon\ll 1$0 is Hermitian and purely off-diagonal: $0<\varepsilon\ll 1$1

The significance of this transformation is structural. The leading $0<\varepsilon\ll 1$2 oscillatory part has been removed analytically, and the residual dynamics is of order $0<\varepsilon\ll 1$3 with bounded coefficients. The solution $0<\varepsilon\ll 1$4 oscillates around $0<\varepsilon\ll 1$5 with amplitude $0<\varepsilon\ll 1$6 uniformly in $0<\varepsilon\ll 1$7, which is the mechanism that permits coarse spatial grids in a regime where the original solution $0<\varepsilon\ll 1$8 is highly oscillatory (Arnold et al., 2024). Recovery of the physical variable is accomplished through

$0<\varepsilon\ll 1$9

A common misunderstanding is that the third-order designation refers to a third-order WKB phase. In the Arnold–Körner scheme it does not: the preprocessing uses a second-order WKB phase $a(x)=E-V(x)$0, while third order refers to the step-size accuracy of the discretization applied to the smoother transformed system.

3. Picard expansion, oscillatory quadrature, and the one-step third-order update

On a single step $a(x)=E-V(x)$1, the transformed system is expanded by Picard iteration: $a(x)=E-V(x)$2 where

$a(x)=E-V(x)$3

equivalently

$a(x)=E-V(x)$4

Truncation at $a(x)=E-V(x)$5 gives

$a(x)=E-V(x)$6

The key difficulty is that the matrices $a(x)=E-V(x)$7 contain highly oscillatory phase factors $a(x)=E-V(x)$8 and nested integrals. The paper develops quadratures for these terms through two mechanisms. The AM (Asymptotic Method) repeatedly integrates by parts in order to increase the $a(x)=E-V(x)$9-order of the remainder. The SAM (Shifted Asymptotic Method) shifts the oscillatory factor so that it has a zero inside the integration interval, specifically at $-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$0, thereby increasing the $-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$1-order of the remainder (Arnold et al., 2024).

The construction uses the step size $-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$2, the phase increment

$-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$3

the recursive coefficients

$-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$4

and the auxiliary functions

$-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$5

For non-oscillatory terms arising in $-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$6, Simpson’s rule is used: $-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$7

The resulting one-step method is organized through three step matrices.

Quantity	Definition	Structural role
$-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$8	$-\frac{\hbar^2}{2m}\,u''(x)+V(x)\,u(x)=E\,u(x)$9	First Picard term
$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$0	$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$1	Second Picard term
$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$2	$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$3	Third Picard term

These matrices are then combined as

$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$4

The physical variables at nodes are recovered by

$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$5

The first quadrature $$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$6 approximates $$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$7 with the mixed estimate

$$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$8

and the final scheme takes $$\varepsilon:=\frac{\hbar}{\sqrt{2m}},\qquad a(x):=E-V(x),\qquad u\equiv\varphi.$$9. The second and third quadratures satisfy

$\varepsilon$0

$\varepsilon$1

The method is therefore genuinely third order in the step size while remaining adapted to the oscillatory structure of the transformed system.

4. Accuracy, stability, numerical behavior, and implementation

Under Hypothesis A, the local discretization error is

$\varepsilon$2

and the global bounds are

$\varepsilon$3

These estimates are uniform in $\varepsilon$4. When $\varepsilon$5, the factor $\varepsilon$6 yields an observed $\varepsilon$7 regime with $\varepsilon$8-independent constant; when $\varepsilon$9, the error behaves like $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$0 (Arnold et al., 2024).

The method is uniformly stable in $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$1, due to the smallness of the step matrices and boundedness of $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$2. No hard CFL-type limit is required. Accuracy improves as $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$3 decreases, but the method is designed to work on coarse grids even as $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$4, including $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$5. The stated practical guidance is that $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$6 should resolve the slow variation of the coefficients and the phase increment $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$7, rather than the original wavelength of $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$8.

The numerical examples illustrate both accuracy and efficiency. For the Airy equation with $\lambda(x)=\frac{2\pi\,\varepsilon}{\sqrt{a(x)}}.$9 on $a\in C^7(I)$0, and for the exponential potential $a\in C^7(I)$1 on $a\in C^7(I)$2 with initial data $a\in C^7(I)$3, $a\in C^7(I)$4, the reported error metric is the max-norm global error in the $a\in C^7(I)$5-variable, $a\in C^7(I)$6. In both examples, WKB3 and the simplified WKB3s show observed $a\in C^7(I)$7 global convergence, consistent with theory. In the regime $a\in C^7(I)$8, WKB3 displays $a\in C^7(I)$9 superconvergence. The paper also reports that work–precision diagrams show WKB3 and WKB3s reach target accuracies with up to $1$00 less CPU time than the second-order WKB2 method, and that for very small $1$01, WKB3 can be up to $1$02 faster at the same accuracy in quad-precision runs (Arnold et al., 2024).

Implementation is $1$03 per step and $1$04 overall. The step formulas require endpoint data and, for the Simpson terms in $1$05, midpoint evaluations. If the phase is available analytically, phase evaluation is $1$06. If it is computed numerically, the paper recommends a high-accuracy spectral procedure: one-time Clenshaw–Curtis $1$07 for the integrals and barycentric interpolation $1$08 per node. This matters because low-order phase computation can dominate the error. If $1$09 is approximated numerically with error $1$10, then the inverse transform contributes a term $1$11 to the $1$12-error, and composite Simpson phase computation can generate an additional $1$13 effect. The reported remedy is spectral phase computation with Clenshaw–Curtis quadrature and barycentric interpolation, which removes the visible degradation in the experiments. The simplified third-order scheme WKB3s uses weaker quadratures and still achieves

$1$14

but without the additional $1$15 gain in the bound.

5. Third-order weighted stationary phase on finite intervals

A distinct but related meaning of third-order phase-integral method appears in the theory of oscillatory integrals. The weighted stationary-phase problem is to asymptotically evaluate

$1$16

when the phase $1$17 has a single stationary point $1$18 and the amplitude $1$19 need not vanish at the endpoints. In the formulation of “Weighted stationary phase of higher orders,” $1$20 is real-valued and $1$21, $1$22 is real-valued and $1$23, $1$24, $1$25 changes sign only once, and the stationary point is non-degenerate in the sense that

$1$26

with $1$27 and $1$28. Degenerate stationary points are not treated (McKee et al., 2016).

The finite-interval structure forces boundary contributions to be retained explicitly. The integration-by-parts hierarchy is

$1$29

The paper’s weighted first derivative test yields

$1$30

with an explicit multi-parameter error term $1$31. This extends Huxley’s first derivative test by preserving boundary terms and sharpening the error on finite intervals.

For the stationary contribution itself, the normalization is chosen so that

$1$32

with $1$33. Writing

$1$34

the coefficients of the asymptotic series are expressed through universal combinatorial combinations $1$35, which encode the Taylor data of both the phase and the weight. The main theorem states that, for $1$36 and sufficiently large $1$37 satisfying $1$38,

$1$39

At third order, the local stationary contribution is organized as

$1$40

where

$1$41

Thus the third-order weighted stationary-phase approximation is

$1$42

The paper provides a detailed remainder estimate involving the parameters $1$43, contributions from the outer intervals controlled by the weighted first derivative test, and local Taylor-probability-integral tail bounds. A worked example with constant weight and quadratic phase shows that when $1$44 and $1$45, all higher $1$46 vanish and only the leading stationary term remains, with the finite-interval boundary terms quantifying the truncation. The authors note applications in harmonic analysis and analytic number theory, including mass distribution of cusp forms, spectral square moments in resonance sums, and improved subconvexity bounds for Rankin–Selberg $1$47-functions.

6. Relation to classical JWKB, admissible regimes, and limitations

Both formulations are phase-integral in the classical sense that the dominant oscillation is encoded through a phase before residual corrections are expanded. In the Schrödinger method, the second-order WKB phase

$1$48

plays the central role. The amplitude prefactor $1$49 is the first transport term, and $1$50 is the standard Langer/Fröman correction entering the $1$51-phase. The numerical scheme can therefore be viewed as combining classical PI/WKB preprocessing with a third-order discretization of the residual smoother system (Arnold et al., 2024).

The weighted stationary-phase method is mathematically rigorous for finite-interval integrals $1$52 with explicit weight dependence through the coefficients $1$53, explicit finite-endpoint boundary terms through the hierarchy $1$54, and quantified remainder bounds based on derivative scales $1$55 (McKee et al., 2016). The resemblance to WKB lies in the inverse-power series structure governed by the local curvature of the phase; the main difference is that finite-interval boundary effects are part of the principal asymptotic statement rather than a peripheral correction.

The admissible regimes are correspondingly different. The numerical Schrödinger method is tailored to non-turning regions with $1$56. It is well suited to initial-value propagation and scattering-type problems on intervals where $1$57. Bound states and quantization conditions are not treated directly; a plausible implication is that shooting-type constructions would require segmentwise propagation together with eigenvalue search, matching, and connection formulae near turning points. The weighted stationary-phase method requires a single non-degenerate stationary point; degenerate stationary points are خارج its stated scope.

The principal limitations are explicit. For the numerical method, turning points are excluded because $1$58 may vanish and the WKB preprocessing breaks down. Remedies listed in the source include switching to a standard ODE solver such as Runge–Kutta in a small neighborhood where $1$59 is small or vanishes, adaptive step-size control, and, in principle, a Langer-type transform, though such a transform is not implemented in the paper. Rapidly varying or discontinuous potentials also lie outside the smooth $1$60 framework required for the third-order quadratures. For the weighted stationary-phase method, the non-degeneracy assumption $1$61 is essential, and no modified coefficient structure for degenerate critical points is given.

A second common misunderstanding is terminological. In the numerical Schrödinger context, third order means third-order accuracy with respect to the step size $1$62, obtained after an analytic phase-integral transformation. In the weighted stationary-phase context, third order means truncation of the local stationary expansion after the $1$63 term, together with explicit boundary contributions and a controlled remainder. The shared phrase reflects a common asymptotic philosophy, but the objects being approximated—solutions of highly oscillatory ODEs in one case and finite-interval oscillatory integrals in the other—are different.