Asymptotic Expansions of Special Functions

A physicist almost never wants a special function at a generic point; they want it in a limit — a Bessel function at large radius, a gamma function at large argument, a Legendre polynomial at large degree, a partition function at large energy. In exactly those regimes the convergent series defining these functions become useless: a Taylor or hypergeometric series near $x=\infty$ converges with agonizing slowness and hides the dominant exponential or oscillatory behavior that carries the physics. Asymptotic analysis supplies what convergence cannot — a short, explicit expansion whose error is sharply controlled after finitely many terms, even though the series usually diverges. Two engines produce nearly every result. The first is the analysis of integral representations: Watson's lemma for Laplace-type integrals, and the method of steepest descent (the saddle-point method) for contour integrals $\int_C e^{x\phi(z)}h(z)\,\mathrm{d}z$ carrying a large parameter $x$ in the exponent. One deforms the contour onto a path of constant phase through a saddle $\phi'(z_0)=0$ , where the integrand collapses to a single Gaussian. The second is the WKB / Liouville–Green treatment of the defining differential equation, which yields uniform expansions valid even across turning points, where Airy functions take over. From these come Stirling's series, the Hankel forms of $J_\nu,Y_\nu,I_\nu,K_\nu$ , and the Plancherel–Rotach asymptotics of orthogonal polynomials — and, underneath them all, the Stokes phenomenon, the discontinuous jump of a subdominant exponential's coefficient that stitches the local expansions into one global function.

Key takeaways

An asymptotic expansion $f(x)\sim\sum_n a_n\phi_n(x)$ is a usually-divergent series whose partial sums approximate $f$ with error smaller than the last retained term as $x\to\infty$ — not convergent, yet enormously accurate.
The method of steepest descent puts the large parameter in an exponential and localizes at a saddle $\phi'(z_0)=0$ , reducing the integral to a Gaussian — one move that yields Stirling's formula and the Bessel/Hankel asymptotics.
Watson's lemma is the endpoint-dominated cousin: the large- $x$ behavior of $\int_0^\infty e^{-xt}g(t)\,\mathrm{d}t$ is fixed by the power law of $g$ at $t=0$ .
WKB / Liouville–Green analysis of the underlying ODE gives uniform large-argument and large-order expansions, with Airy functions bridging turning points.
Divergence is a feature: terms shrink to a minimum near $n\approx x$ then grow factorially; truncating there gives an error as small as $e^{-2\pi x}$ , and two functions can share a series yet differ by such amounts — the Stokes phenomenon.
Key idea: put the large parameter into an exponential, $\int e^{x\phi(z)}h(z)\,\mathrm{d}z$ , then localize at the stationary point $\phi'(z_0)=0$ where the exponent is locally quadratic and the integral becomes Gaussian. Stirling, Hankel, and Plancherel–Rotach are all bookkeeping around that one Gaussian collapse.

The core insight

Put the large parameter into an exponential, $\int e^{x\phi(z)}h(z)\,\mathrm{d}z$ , then localize: the integral is dominated by the neighborhood of the stationary point $\phi'(z_0)=0$ , where the exponent is locally quadratic and the integral becomes Gaussian. Everything else — Stirling, Hankel, Plancherel–Rotach — is bookkeeping around this one Gaussian collapse.

Where this lives today

The saddle-point method is a daily instrument in physics. In quantum theory it is the semiclassical limit: stationary-phase evaluation of the Feynman path integral is steepest descent with $x=1/\hbar$ , and the WKB connection formulae — the Airy asymptotics of equation (6) — fix the quantized energies of every bound state through Bohr–Sommerfeld $\oint p\,\mathrm{d}q=2\pi\hbar(n+\tfrac12)$ , the $+\tfrac12$ Maslov index being the accumulated $\pi/4$ turning-point phases. In statistical mechanics the density of states $\rho(E)=\frac{1}{2\pi i}\int e^{\beta E}Z(\beta)\,\mathrm{d}\beta$ is a saddle-point inversion: the saddle $\beta^\ast$ is the inverse temperature and the Gaussian curvature gives the entropy fluctuations — the computation behind the Hardy–Ramanujan partition asymptotic and the Cardy formula in 2D conformal field theory. In scattering, the high-energy partial-wave sum is dominated by large $\ell$ , where Hilb's Legendre asymptotics convert it into the eikonal impact-parameter representation, $\ell+\tfrac12\to kb$ . In random-matrix theory the spectral edge of a large Gaussian ensemble is governed by the Hermite polynomials at their turning point, where Plancherel–Rotach asymptotics produce the Airy kernel and the Tracy–Widom distribution. And resurgence theory uses the divergent tails themselves — via Borel summation — to reconstruct non-perturbative physics in field theory and matrix models, turning divergence into a tool.

The physical picture

Think of the exponent as terrain. Over the complex $z$ -plane, plot the height $\operatorname{Re}[x\phi(z)]$ ; the integrand's magnitude is $e^{\text{height}}$ , so a contour integral is dominated by the highest ground the contour must cross. But the contour is free to deform, and it wants the lowest high point — a mountain pass, a saddle. At a saddle $z_0$ the surface curves up one way and down the perpendicular way; the optimal route crosses the pass and descends into the valleys along the line of steepest descent, where the phase $\operatorname{Im}[x\phi]$ stays constant so nothing oscillates. Near the saddle the height is a downturned parabola — a clean Gaussian. That is the whole method in one image.

Figure

The steepest-descent landscape: height is Re xφ(z), the original contour deforms onto the green path of constant phase through the saddle z₀ where φ′(z₀)=0, and near that pass the integrand becomes a Gaussian of width set by the curvature φ″(z₀)

Divergence is the second thing to see, and it is not a defect. Track the terms $a_n/x^n$ : they shrink at first, reach a smallest term near $n\approx x$ , then grow factorially. Truncating at that minimum extracts the best the series can offer, with a leftover error of the same order as the smallest term — for Stirling, astonishingly small, $e^{-2\pi x}$ . Beyond that floor the series actively lies. Two functions can therefore share a series yet differ by exponentially small amounts it can never resolve — the seed of the Stokes phenomenon.

Figure

Stirling's series is divergent: the term magnitudes for x = 2 fall to a smallest term ≈ 8×10⁻⁷ near k = 7, then grow without bound — truncating at that minimum pins the error near the floor e^{−2πx} ≈ 3.5×10⁻⁶

Finally, large argument and large order are distinct limits. For $J_\nu(z)$ with $\nu$ fixed and $z\to\infty$ the function oscillates — a propagating wave. Hold $z/\nu$ near $1$ and you probe a turning point, where the two saddles collide and Airy functions govern the crossover — a classical particle at a turning radius. Same integral, different saddle geometry, different physics.

The core equations

The subject rests on one definition and a handful of master formulas. Poincaré's definition makes a divergent series a legitimate approximation and fixes the tilde: for each fixed order $N$ the remainder must vanish faster than the last retained power as $x\to\infty$ — it says nothing about $N\to\infty$ at fixed $x$ , which is why the series may diverge.

f(x)\sim\sum_{n=0}^{\infty} a_n\,x^{-n}\ \ (x\to\infty) \iff f(x)-\sum_{n=0}^{N}a_n x^{-n}=o(x^{-N}) \qquad (1)

where: $f(x)$ — the function being approximated as $x\to\infty$ ; $a_n$ — expansion coefficients, fixed numbers independent of $x$ ; $o(x^{-N})$ — remainder smaller than the last kept term, at fixed $N$ . It is a promise about the error at fixed order, never an equality. Watson's lemma is the simplest engine, for integrals where the large parameter multiplies a linear exponent and the action sits at the endpoint. The whole large- $x$ expansion is read off from the power-law behavior of $g$ at $t=0$ ; everything else is exponentially suppressed.

\int_0^{\infty} e^{-x t}\,g(t)\,\mathrm{d}t \sim \sum_{n=0}^{\infty} \frac{a_n\,\Gamma(\lambda+n+1)}{x^{\lambda+n+1}},\qquad g(t)\sim t^{\lambda}\sum_{n}a_n t^{n} \qquad (2)

where: $g(t)$ — amplitude with a power-law expansion at $t=0$ , $\lambda>-1$ ; $x$ — large parameter in the exponent; $\Gamma$ — gamma function, from integrating each monomial term by term. Each power near the endpoint contributes one gamma factor; the resulting series is asymptotic, generally divergent. The general engine is the method of steepest descent. For an analytic exponent the dominant contribution comes from an interior saddle $\phi'(z_0)=0$ ; expanding to second order turns the integral into a Gaussian of width set by $\phi''(z_0)$ .

\int_C e^{x\,\phi(z)}h(z)\,\mathrm{d}z \sim \sqrt{\frac{2\pi}{x\,|\phi''(z_0)|}}\;e^{x\,\phi(z_0)}\,h(z_0)\,e^{i\theta},\quad \phi'(z_0)=0 \qquad (3)

where: $z_0$ — saddle point, root of $\phi'(z_0)=0$ ; $\phi''(z_0)$ — curvature at the saddle, sets the Gaussian width and amplitude; $\theta$ — orientation of the steepest-descent path, $\theta=-\tfrac12\arg\phi''(z_0)$ ; $h(z)$ — slowly varying amplitude at the saddle. A real saddle gives monotone exponential behavior; a complex-conjugate pair interferes to give oscillation — the distinction separating $I_\nu,K_\nu$ from $J_\nu,Y_\nu$ . Applying (3) to the gamma function's Euler integral to all orders produces Stirling's series, its corrections governed by the Bernoulli numbers — the archetype of a divergent asymptotic series.

\ln\Gamma(x)\sim \Big(x-\tfrac12\Big)\ln x - x + \tfrac12\ln(2\pi) + \sum_{k=1}^{\infty}\frac{B_{2k}}{2k(2k-1)\,x^{2k-1}} \qquad (4)

where: $B_{2k}$ — Bernoulli numbers, $B_2=\tfrac16,\ B_4=-\tfrac{1}{30},\ B_6=\tfrac{1}{42}$ ; $x$ — large argument, valid for $|\arg x|<\pi$ . Exponentiating the first three terms gives $\Gamma(x)\sim\sqrt{2\pi}\,x^{x-1/2}e^{-x}(1+\tfrac{1}{12x}+\tfrac{1}{288x^2}-\cdots)$ , the $\tfrac1{12x}$ correction coming from $B_2$ . The same method applied to the Sommerfeld integral for the Hankel function delivers the large-argument form of every Bessel-type function; the outgoing $H^{(1)}_\nu$ carries the oscillatory phase and $x^{-1/2}$ decay of a cylindrical wave.

H^{(1)}_\nu(x)\sim \sqrt{\tfrac{2}{\pi x}}\;e^{\,i\left(x-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)}\sum_{k=0}^{\infty}\frac{i^k a_k(\nu)}{x^k},\qquad a_k(\nu)=\frac{\prod_{j=1}^{k}\big(4\nu^2-(2j-1)^2\big)}{k!\,8^{k}} \qquad (5)

where: $H^{(1)}_\nu$ — Hankel function of the first kind, an outgoing cylindrical wave; $a_k(\nu)$ — Bessel coefficients, $a_0=1,\ a_1=\tfrac{4\nu^2-1}{8}$ ; $x$ — large argument, $x\gg\nu^2$ . Real and imaginary parts recover $J_\nu,Y_\nu$ ; the product terminates only for half-integer $\nu$ , why the spherical Bessel functions are elementary. Where two saddles coalesce — at a turning point — the Gaussian model fails and the uniform solution is the Airy function.

\operatorname{Ai}(x)\sim \frac{e^{-\frac{2}{3}x^{3/2}}}{2\sqrt{\pi}\,x^{1/4}}\ (x\to+\infty),\qquad \operatorname{Ai}(-x)\sim \frac{\sin\!\left(\frac{2}{3}x^{3/2}+\frac{\pi}{4}\right)}{\sqrt{\pi}\,x^{1/4}} \qquad (6)

where: $\operatorname{Ai}$ — Airy function, the universal turning-point solution; $x^{3/2}$ — the WKB phase/exponent integral $\int\sqrt{|x'|}\,\mathrm{d}x'$ . The exponential decay on one side and the oscillation with a $\pi/4$ phase on the other are the WKB connection formula — the same $\pi/4$ that becomes the Maslov index in Bohr–Sommerfeld quantization.

Derivation AStirling's formula by steepest descent

The goal is to obtain the leading Stirling approximation $\Gamma(x+1)\sim\sqrt{2\pi x}\,(x/e)^x$ directly from the Euler integral, by recognizing it as an exponential integral with a single real saddle and reducing the neighborhood of that saddle to a Gaussian. This is the cleanest demonstration of the steepest-descent philosophy: cast the object as $e^{f}$ , find where $f$ is stationary, and integrate the parabola there — then read off how the higher corrections build the Bernoulli series of equation (4).

Step 1. Write the factorial as an integral and pull the algebraic factor into the exponent: $\Gamma(x+1)=\int_0^\infty t^{x}e^{-t}\,\mathrm{d}t=\int_0^\infty e^{\,x\ln t-t}\,\mathrm{d}t$ , so the integrand is $e^{f(t)}$ with $f(t)=x\ln t-t$ . Why: with the large parameter multiplying $\ln t$ , the integrand is sharply peaked and its peak drifts with $x$ — the hallmark that saddle-point applies.

Step 2. Locate the maximum: $f'(t)=x/t-1=0\Rightarrow t_0=x$ , and $f''(x)=-1/x$ . Why: the peak sits at $t_0=x$ , so the integrand is concentrated within order $\sqrt{x}$ of $t=x$ , and the negative curvature sets the Gaussian width.

Step 3. Expand $f$ to second order about $t_0=x$ , extend the limits to the whole line (the tails are exponentially negligible), and do the Gaussian: $\Gamma(x+1)\approx e^{f(x)}\int_{-\infty}^{\infty}e^{-(t-x)^2/2x}\,\mathrm{d}t=e^{\,x\ln x-x}\sqrt{2\pi x}$ . Why: $f(x)=x\ln x-x$ is the leading exponent and the Gaussian supplies $\sqrt{2\pi x}$ — exactly formula (3) with $|\phi''|=1/x$ .

Step 4. Rearranging gives the boxed form; retaining the cubic and quartic terms of $f$ produces the $1+\tfrac1{12x}+\cdots$ corrections, i.e. the Bernoulli series of equation (4). Why: each higher derivative of $f$ at the saddle feeds one more term of the divergent tail.

\Gamma(x+1)=x!\sim \sqrt{2\pi x}\;\Big(\frac{x}{e}\Big)^{x}

Verified: $\Gamma(11)=10!=3\,628\,800$ ; the leading term $\sqrt{20\pi}(10/e)^{10}=3.5987\times10^6$ is 0.83% low, and $\times(1+\tfrac1{120})$ gives $3.62868\times10^6$ — four significant figures.

Go deeper · Derivation B · 1 minWKB form of Bessel's equation by the Liouville transformation

Derivation BWKB form of Bessel's equation by the Liouville transformation

The aim is to recover the large-argument amplitude $x^{-1/2}$ and the oscillatory phase of the Bessel functions straight from the differential equation, without any integral representation, by using the Liouville substitution to strip the first-derivative term and expose a constant-coefficient equation at leading order. This makes plain why every solution of Bessel's equation decays as $x^{-1/2}$ and oscillates with unit frequency far from the origin — the manipulation that underlies the whole WKB method.

Step 1. Start from Bessel's equation $x^2 y''+x y'+(x^2-\nu^2)y=0$ , with solutions $J_\nu,Y_\nu$ . Why: the $xy'$ term is what obscures the large- $x$ behaviour and must be removed first.

Step 2. Substitute $y=x^{-1/2}u$ ; the first-derivative term cancels, leaving $u''+\big(1-\tfrac{\nu^2-1/4}{x^2}\big)u=0$ . Why: this is Schrödinger normal form with effective "energy" $1$ and a potential dying like $1/x^2$ — the $x^{-1/2}$ prefactor is now visibly a WKB amplitude.

Step 3. For $x\to\infty$ the $1/x^2$ potential is negligible, leaving $u''+u\approx0$ , so $u\sim A\cos x+B\sin x$ . Why: at leading order it is a free oscillator of unit frequency, giving both the oscillation and the decay $y\sim x^{-1/2}(A\cos x+B\sin x)$ .

Step 4. Fixing the constants by matching to the integral representation yields the standard normalization and phase. Why: the ODE alone cannot set the absolute normalization or the $-\nu\pi/2-\pi/4$ phase; those come from the global connection of equation (5).

J_\nu(x)\sim\sqrt{\tfrac{2}{\pi x}}\cos\!\big(x-\tfrac{\nu\pi}{2}-\tfrac{\pi}{4}\big)\quad (x\to\infty)

Verified: at $x=20,\nu=0$ the leading form gives $\sqrt{2/(20\pi)}\cos(20-\pi/4)=0.16666$ against the true $J_0(20)=0.16702$ — three significant figures, about 0.2%.

Go deeper · Derivation C · 1 minOptimal truncation and the exponentially small error floor

Derivation COptimal truncation and the exponentially small error floor

The goal is to show quantitatively why an asymptotic series is best stopped at its smallest term, and to estimate the residual error there. The engine is the factorial growth of the Bernoulli numbers driving Stirling's series: since the numerators grow like a factorial while the denominators grow only geometrically in $x$ , the terms must turn around, and the location and size of that turnaround pin down both the optimal truncation order and the accuracy floor.

Step 1. Estimate the general term. With $|B_{2k}|\sim 2(2k)!/(2\pi)^{2k}$ , the $k$ -th term behaves as $2\,(2k-2)!/[(2\pi)^{2k}x^{2k-1}]$ . Why: the factorial $(2k-2)!$ is what makes the series diverge — no power of $x$ can outrun it as $k\to\infty$ .

Step 2. Minimize over $k$ : $\frac{\mathrm d}{\mathrm dk}\ln[(2k-2)!/(2\pi x)^{2k}]=0$ with $\ln(2k-2)!\approx 2k\ln(2k)$ gives the smallest term near $k\approx\pi x$ . Why: terms shrink while the geometric factor wins, then grow once the factorial takes over; the crossover is where to stop.

Step 3. Evaluate the smallest term at $k\approx\pi x$ ; the factorial collapses it to order $e^{-2\pi x}$ . Why: the optimally truncated Stirling series carries an irreducible error of order $e^{-2\pi x}$ — a hard floor no rearrangement of the same series can beat.

|R_{N^\ast}(x)|\;\sim\;\min_{N}\frac{|B_{2N}|}{2N(2N-1)x^{2N-1}}\;\approx\;e^{-2\pi x}

Verified: the term data for $x=2$ (second figure) fall to a smallest term $\approx 8\times10^{-7}$ near $k=7$ , of the same order as the floor $e^{-4\pi}\approx3.5\times10^{-6}$ , matching Bender–Orszag §3.6 and DLMF §5.11.

Limiting cases

The one saddle-point machinery specializes to each special function as the number, location, and reality of the saddles change — real saddles giving exponentials, complex pairs oscillation, coalescing saddles the Airy form.

Case	Result	Meaning
Gamma, large argument	$\Gamma(x)\sim\sqrt{2\pi}\,x^{x-\frac12}e^{-x}\big(1+\tfrac{1}{12x}+\tfrac{1}{288x^2}-\cdots\big)$	single real saddle at $t=x$ ; Stirling's divergent series
Bessel, large argument	$J_\nu(x)\sim\sqrt{\tfrac{2}{\pi x}}\cos\!\big(x-\tfrac{\nu\pi}{2}-\tfrac{\pi}{4}\big)$	complex-conjugate saddles interfere; oscillation, $x\gg\nu^2$
Modified Bessel	$I_\nu(x)\sim\dfrac{e^{x}}{\sqrt{2\pi x}},\quad K_\nu(x)\sim\sqrt{\tfrac{\pi}{2x}}\,e^{-x}$	real saddles; pure exponential growth/decay, no oscillation
Legendre, large degree	$P_n(\cos\theta)\sim\sqrt{\tfrac{2}{\pi n\sin\theta}}\cos\!\big[(n+\tfrac12)\theta-\tfrac{\pi}{4}\big]$	Laplace–Heine saddle of the generating integral; $0<\theta<\pi$ bulk
Hermite, oscillatory bulk	$e^{-x^2/2}H_n(x)\sim\big(\tfrac{2}{\pi}\big)^{1/4}(2n)^{n/2}n^{-1/4}e^{-n/2}\cos(\cdots)$	two saddles inside $\|x\|<\sqrt{2n}$ ; Plancherel–Rotach, Airy at the edge

Common trapTreating the tilde as an equals sign

The single most common error is to read $f(x)\sim\sum a_n x^{-n}$ as a convergent identity and "improve" the answer by summing more terms. An asymptotic series is generally divergent: the terms shrink only until the smallest one, near $n\approx x$ and of size $\sim e^{-x}$ (or $\sim e^{-2\pi x}$ for Stirling's optimally truncated series), after which they grow without bound and the partial sums walk away from the truth. The tilde constrains the error at fixed truncation order — never a licence to take $N\to\infty$ . A subtler trap: because the series is blind to exponentially small terms, it cannot by itself decide a Stokes jump — the coefficient of a subdominant exponential lives below the resolution of the dominant series, and fixing it needs resurgence, Borel summation, or a matching condition, not more terms.

Worked examples

1 · Stirling's approximation for 10! (fully worked)

Use Stirling's series to estimate $10!=\Gamma(11)$ (so $x=10$ ) and compare against the exact $3\,628\,800$ term by term.

Step 1. Leading term. At $x=10$ the prefactor is $\sqrt{20\pi}=7.9267$ and $(10/e)^{10}=4.5400\times10^{5}$ , giving $\sqrt{20\pi}(10/e)^{10}=\mathbf{3.5987\times10^{6}}$ . Against the exact $3.6288\times10^{6}$ this is 0.83% low — good for two lines of algebra.

Step 2. First correction. The $B_2$ term exponentiates to $(1+\tfrac{1}{12x})=1.008333$ ; multiplying gives $3.5987\times10^{6}\times1.008333=\mathbf{3.62868\times10^{6}}$ . The error collapses from 0.83% to $3\times10^{-5}$ — four significant figures from one term.

Step 3. Second correction. Including $B_4$ via $\tfrac{1}{288x^2}=3.472\times10^{-5}$ gives $\mathbf{3.628806\times10^{6}}$ , within $3$ parts in $10^{6}$ of the exact $3\,628\,800$ . More terms keep improving it — but only because $x=10$ puts the smallest term (near $k\approx\pi x\approx31$ ) far away. Push past it and the series would diverge.

2 · The large-argument Bessel function $J_0(20)$

Estimate $J_0(20)$ from the asymptotic form and see how one correction sharpens it; the true value is $J_0(20)=0.167025$ .

Step 1. Leading oscillatory term. With $\nu=0$ the phase is $\omega=x-\pi/4=19.2146$ rad and the amplitude $\sqrt{2/(\pi x)}=0.17841$ . Then $J_0(20)\approx0.17841\cos(19.2146)=0.17841\times0.9341=\mathbf{0.16666}$ , correct to three significant figures (0.2% high).

Step 2. One correction term. To next order $J_0(x)\sim\sqrt{2/\pi x}\,[\cos\omega\,(1-\tfrac{9}{128x^2})+\sin\omega\,\tfrac{1}{8x}]$ . With $\tfrac{1}{8x}=0.00625$ and $\sin\omega=0.3572$ , the bracket is $0.9341\times0.999824+0.3572\times0.00625=0.93617$ , so $J_0(20)\approx0.17841\times0.93617=\mathbf{0.167025}$ .

Step 3. Interpret. One correction term drives the estimate from three figures to six significant figures — why physicists reach for the asymptotic form over the convergent-but-hopelessly-slow power series at large argument. Leading: 0.16666; one correction: 0.167025; exact: 0.167025.

Physics of the portrait

Asymptotic analysis grew out of one stubborn puzzle — how to give meaning to divergent series that were nonetheless astonishingly accurate — and no one wrestled with it more consequentially than George Gabriel Stokes, who found that the constants in such an expansion could jump.

Origins · 1 minGeorge Gabriel Stokes and how this was discovered

George Gabriel Stokes

1819–1903 · asymptotic expansions, the Stokes phenomenon, and steepest descent

Stokes came to asymptotics through optics. Studying the intensity near a rainbow's edge, he took up Airy's integral for the diffraction pattern and needed its behavior for large argument on both sides of the caustic. The convergent power series was numerically useless, so Stokes worked with the divergent asymptotic expansion — and in 1857 noticed something disturbing. As the argument rotated through the complex plane, the constant multiplying the exponentially small (subdominant) solution changed discontinuously, jumping as the ray crossed certain lines, now called Stokes lines. Far from a paradox, this is how the local expansions of a single analytic function knit together globally: a series determines its function only up to exponentially small terms, and those hidden terms surface as the jump. The insight sat quietly for a century until, in the 1980s–90s, Écalle's resurgence and Berry's demonstration that the jump is a smooth error-function transition turned Stokes's puzzle into a precise, calculable structure. Everything here — Stirling's divergent tail, the WKB connection formulae, the way $J_\nu$ and $H^{(1)}_\nu$ trade dominance — is a chapter of the story Stokes opened at that rainbow.

History

1730

— de Moivre and Stirling fix the constant √(2π) in the factorial asymptotic

1857

— Stokes discovers the discontinuous jump of subdominant constants across Stokes lines

1886

— Poincaré gives the rigorous definition of an asymptotic (divergent) series

1909

— Debye develops the method of steepest descent for the large-order Bessel functions

1918

— Watson proves the endpoint lemma bearing his name and completes the Bessel treatise

Going further

Texts. Bender & Orszag, Advanced Mathematical Methods for Scientists and Engineers, ch. 6 (the definitive physics treatment) · Stone & Goldbart, Mathematics for Physics, asymptotic-methods ch. · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., chs. 12, 14 · Hassani, Mathematical Physics · Nakahara, Geometry, Topology and Physics (path-integral/steepest-descent) · Byron & Fuller, Mathematics of Classical and Quantum Physics; Riley, Hobson & Bence, Mathematical Methods, 3rd ed. · Olver, Asymptotics and Special Functions; NIST DLMF chs. 2, 5, 9, 10.

Provenance. A verifier should recompute: (1) $\Gamma(11)=3\,628\,800$ against $\sqrt{20\pi}(10/e)^{10}=3.5987\times10^6$ and the $(1+\tfrac1{120})$ correction $3.62868\times10^6$ ; (2) the leading Bessel term $\sqrt{2/(20\pi)}\cos(20-\pi/4)=0.16666$ and one-correction value $0.167025$ against $J_0(20)=0.167025$ ; (3) $B_2=\tfrac16,B_4=-\tfrac1{30}$ and Stirling coefficients $\tfrac1{12},\tfrac1{288}$ ; (4) $a_1(\nu)=(4\nu^2-1)/8$ ; (5) the divergence data for $x=2$ (smallest term $\approx8\times10^{-7}$ near $k=7$ , floor $e^{-4\pi}\approx3.5\times10^{-6}$ ). Forms follow Bender–Orszag chs. 3 & 6, Stone & Goldbart, Arfken–Weber–Harris chs. 12/14, and DLMF chs. 5, 9, 10; all numbers are standard tabulated values.

Prerequisites. Gamma Function · Bessel Functions of the First and Second Kind · Hankel Functions and Asymptotics · Orthogonal Polynomials · Generating Functions. Leads to. Modified Bessel Functions · Legendre Polynomials · Spherical Bessel Functions · Associated Legendre Functions.

Ask the tutor: Ask me to walk through any saddle-point evaluation step by step — Stirling, a Bessel or Hankel asymptotic, or the Airy turning-point connection — and I will show you where the Gaussian collapse happens.

What does the tilde in f(x) ∼ Σ aₙx⁻ⁿ actually mean?

Something precise but weaker than equality. For each fixed number of terms $N$ , the difference between $f$ and the partial sum vanishes faster than the last kept power as $x\to\infty$ : $f(x)-\sum_0^N a_n x^{-n}=o(x^{-N})$ . The statement is about $x\to\infty$ at fixed $N$ , and says nothing about $N\to\infty$ at fixed $x$ — which is why the series can be wildly accurate with a few terms yet diverge if you keep going. Read the tilde as "the error is smaller than the last term I kept", never as "=".

If the series diverges, how can it beat a convergent one?

Because accuracy and convergence are different things. A convergent series (say the Taylor series of $J_0$ ) reaches the answer eventually, but at large argument it may need hundreds of enormous, nearly-cancelling terms. An asymptotic series gives most of that accuracy in two or three terms, each tuned to the dominant large- $x$ structure. The catch is a floor: the error cannot go below the smallest term, of order $e^{-x}$ (or $e^{-2\pi x}$ for Stirling). For physics that floor sits far below measurement precision, so the divergent series is the practical tool and the convergent one the impractical curiosity.

When do I use Watson's lemma and when full steepest descent?

Use Watson's lemma when your integral already has the Laplace form $\int_0^\infty e^{-xt}g(t)\,\mathrm{d}t$ , with the large parameter multiplying a linear exponent and the action at the endpoint $t=0$ ; then the whole expansion drops out of the power-law behavior of $g$ near zero. Use full steepest descent when the exponent is a general analytic $x\phi(z)$ and the action is at an interior saddle $\phi'(z_0)=0$ ; there you deform onto the constant-phase path before expanding. In practice steepest descent often reduces to Watson's lemma once you straighten the contour near the saddle.

Why do Bessel functions oscillate but modified Bessel functions do not?

It comes down to whether the saddles are real or complex. For $J_\nu$ and $Y_\nu$ the saddles form a complex-conjugate pair; adding the two conjugate contributions gives a cosine — oscillation and the $x^{-1/2}$ decaying wave. For $I_\nu$ and $K_\nu$ the saddles sit on the real axis, so each gives a pure real exponential: $I_\nu$ grows like $e^{x}/\sqrt{2\pi x}$ , $K_\nu$ decays like $\sqrt{\pi/2x}\,e^{-x}$ , no oscillation. The factor of $i$ that turns Bessel's equation into the modified one moves the saddles from the imaginary to the real axis — the whole difference between a propagating and an evanescent wave.

What is the Stokes phenomenon, and does it make the function ill-defined?

No — the function stays perfectly analytic and single-valued; only its asymptotic description jumps. As you rotate the phase, one exponential is dominant and another subdominant (exponentially smaller). Across a Stokes line the coefficient of the subdominant exponential changes discontinuously — allowed because it controls a term below the resolution of the dominant series, in the $e^{-x}$ gap the leading asymptotics cannot see. Berry showed in 1989 that the "jump" is really a smooth error-function transition of width $\sqrt{x}$ in the phase, and resurgence goes further: the divergent tail of one series encodes exactly the subdominant exponential it seemed to ignore.

Can a divergent asymptotic series be summed to recover the exact function?

Often yes, with the right resummation. The standard tool is Borel summation: divide the $n$ -th coefficient by $n!$ to tame the factorial growth, sum the resulting convergent Borel series, then Laplace-transform back. When the Borel transform has no singularities on the positive axis this recovers the exact function; when it does, those singularities encode the non-perturbative, exponentially small terms — the subdominant exponentials of the Stokes phenomenon. Écalle's resurgence organizes this into "trans-series". This is how modern physics extracts instanton effects straight from a divergent expansion — the divergence is the signal, not the noise.

Asymptotic Expansions of Special Functions

Where this lives today

The physical picture

The core equations

Limiting cases

Worked examples

1 · Stirling's approximation for 10! (fully worked)

2 · The large-argument Bessel function $J_0(20)$

Physics of the portrait

Going further

Before this

Where this leads

Asymptotic Expansions of Special Functions

Where this lives today

The physical picture

The core equations

Limiting cases

Worked examples

1 · Stirling's approximation for 10! (fully worked)

2 · The large-argument Bessel function J0(20)J_0(20)J0​(20)

Physics of the portrait

Going further

Before this

Where this leads

2 · The large-argument Bessel function $J_0(20)$