Defines the order symbols and Poincare's notion of an asymptotic expansion, shows that such an expansion determines its coefficients but not its function, says which operations on it are legitimate, and gives the optimal-truncation rule that makes a divergent series into a numerical method.
Let f and g be defined on a neighbourhood of a point x0, which may be +∞, and suppose g does not vanish there. All three relations below are statements about the limit x→x0 and about nothing else; none of them says anything at a fixed x.
f(x)=O(g(x)) means that ∣f/g∣ is bounded near x0: there are C and a neighbourhood on which ∣f(x)∣≤C∣g(x)∣.
f(x)=o(g(x)) means f(x)/g(x)→0 as x→x0.
f(x)∼g(x) means f(x)/g(x)→1 as x→x0; equivalently f=g+o(g).
A sequence of functions φ0,φ1,φ2,…, none vanishing near x0, is an asymptotic sequence at x0 if
φn+1(x)=o(φn(x))(x→x0)for every n.
The two cases used in this part of the handbook are φn(x)=x−n as x→+∞ and φn(x)=(x−x0)n as x→x0. Each term of an asymptotic sequence is negligible beside the one before it, and that is the only ordering the symbol ∼ will ever rest on.
Definition(Poincare asymptotic expansion)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.6, Eqs. (12.95)-(12.99) and footnotes 4-5, p. 583; Ahlfors (1979), Ch. 5, §2.5, p. 205
Let f be defined for all large real x and let a0,a1,a2,… be constants. Write
sN(x)=n=0∑Nxnan,RN(x)=f(x)−sN(x).
We say that ∑n≥0anx−n is an asymptotic expansion of f as x→+∞, and write
f(x)∼n=0∑∞xnan(x→+∞),
when
x→∞limxNRN(x)=0for each fixed N.(P)
Equivalently: RN(x)=o(x−N) for each N, so the error made by stopping after the x−N term is smaller than the last term kept. The sign ∼ is not an equality. It is a statement about x→∞ at fixed N, and it carries no information at a fixed x.
Divergence is normal and is not part of the definition. Nothing in (P) asks the series ∑anx−n to converge, and for the expansions of interest it does not: the coefficients an typically grow like n!, so the series diverges for every x. Arfken calls such a series semiconvergent(Arfken, Weber & Harris (2013), p. 580).
A conflict of definition, and which we follow. Arfken's formal definition of an asymptotic series imposes (P) and a second requirement, that xNRN(x)→∞ as N→∞ at fixed x, which excludes convergent series of inverse powers (Arfken, Weber & Harris (2013), Eqs. (12.97)-(12.98), p. 583); his own footnote records that some authors regard that exclusion as artificial. Ahlfors, characterising the remainder in Stirling's formula, uses the fixed-order limit alone (Ahlfors (1979), p. 205). We follow Ahlfors: in this handbook (P) is the definition, a convergent inverse-power series is allowed to count as its own asymptotic expansion, and when divergence matters we say so. Readers taking a statement from Arfken should keep his extra clause in mind.
More generally, if {φn} is an asymptotic sequence at x0, then f∼∑anφn means
f(x)−n=0∑Nanφn(x)=o(φN(x))(x→x0)for each N.
Expansions in half-integer powers, and expansions carrying a prefactor such as e−x or exf(t0), are of this kind; see Watson's lemma and Laplace's method.
(The coefficients are unique.) If f∼∑n≥0anx−n and also f∼∑n≥0bnx−n as x→+∞, then an=bn for every n. A function has at most one asymptotic power-series expansion.
(The function is not.) If f∼∑n≥0anx−n and h(x)=o(x−N) for every N, then f+h has the same expansion. In particular f and f+e−x are indistinguishable at this level of description, so infinitely many functions share one expansion.
(Existence can fail.) Not every function has an asymptotic power series; Arfken's example is ex.
Riley, Hobson and Bence make the same two-sided point for expansions in the complex plane: at a fixed value of argz the expansion is unique, but one expansion can stand for more than one function, and one function can need different expansions in different sectors.
Idea of the proofOriginal
Both halves of the uniqueness statement come from the same observation: condition (P) pins down each coefficient as a limit, and a limit cannot see anything that decays faster than every power. Taking x→∞ in f∼a0+a1/x+⋯ gives a0=limf; subtracting a0, multiplying by x and taking the limit again gives a1; and so on. That recursion determines the an from f alone, which is part 1. It also shows what the recursion is blind to: a function h with xNh(x)→0 for every N contributes zero at every stage, which is part 2. Part 3 is the first step of the recursion failing, since ex has no finite limit.
ProofOriginal
The coefficients as limits. Suppose f∼∑anx−n. We show by induction on N that
aN=x→∞limxN(f(x)−n=0∑N−1anx−n),
the empty sum being 0. For N=0: condition (P) with N=0 says f(x)−a0→0, so a0=limf(x), which is the claim. For the step, assume the formula for 0,…,N−1. Condition (P) at order N says xN(f−sN)→0, that is
xN(f(x)−n=0∑N−1anx−n)−aN⟶0,
which is the claim at order N. Each aN is therefore a limit computed from f and the earlier coefficients, and the earlier coefficients are themselves so computed. If ∑bnx−n is also an asymptotic expansion of f, the same limits produce the bn; since a limit is unique, bn=an for all n by induction. This is part 1.
Adding an invisible function. Let h satisfy xNh(x)→0 for every N. With R~N=(f+h)−sN=RN+h we get xNR~N=xNRN+xNh→0+0=0 for each fixed N, so f+h satisfies (P) with the same coefficients. The function h(x)=e−x qualifies, because xNe−x→0 for every N. This is part 2.
A function with no expansion. If f∼∑anx−n then by step 1 the number a0=limx→∞f(x) exists and is finite. Since ex→∞, no such a0 exists, so ex has no asymptotic power-series expansion at +∞. This is part 3. □
Linear combinations.αf+βg∼∑n≥0(αan+βbn)x−n for all constants α,β.
Products.fg∼∑n≥0cnx−n with cn=∑k=0nakbn−k: the Cauchy product of two asymptotic expansions is the asymptotic expansion of the product.
Integration. If in addition a0=a1=0 and f is locally integrable on [x1,∞) for some x1, then ∫x∞f(t)dt converges for large x and
∫x∞f(t)dt∼n≥2∑n−1anx−(n−1).
Differentiation is not permitted. From f∼∑anx−n one may not conclude f′∼−∑nanx−n−1. Arfken states this: term-by-term differentiation holds only under special extra conditions. If it is known in advance that f′ has some asymptotic power-series expansion, then that expansion is the differentiated one.
Arfken states 2 and 3 without the proviso a0=a1=0 in 3; we have added it because without it ∫x∞f need not exist, and we prove the statement in the form given here.
ProofOriginal
Throughout, sN and tN denote the partial sums of the two expansions and RN=f−sN, SN=g−tN, so RN(x)=o(x−N) and SN(x)=o(x−N) for each fixed N.
(αf+βg)−(αsN+βtN)=αRN+βSN=o(x−N), and αsN+βtN is the N-th partial sum of ∑(αan+βbn)x−n.
Write uN(x)=∑n=0Ncnx−n with cn the Cauchy coefficients. Then
The product sNtN is a finite sum of terms akblx−(k+l) with k,l≤N; those with k+l≤N assemble to exactly uN, so sNtN−uN is a finite sum of powers x−m with m>N, hence O(x−N−1)=o(x−N). Next, sN and tN are bounded for large x (each tends to a0, resp. b0), so sNSN=o(x−N) and tNRN=o(x−N); and RNSN=o(x−N)o(x−N)=o(x−N). Adding the four pieces gives fg−uN=o(x−N).
Fix N≥2. Because a0=a1=0, the partial sum is sN(t)=∑n=2Nant−n, and f(t)=sN(t)+RN(t) with RN(t)=o(t−N). In particular f(t)=O(t−2) for large t, so ∫x∞f converges absolutely. Integrating the partial sum exactly,
∫x∞sN(t)dt=n=2∑Nn−1anx−(n−1),
which is the (N−1)-st partial sum of the claimed expansion. For the remainder, let ε>0 and choose X with ∣RN(t)∣≤εt−N for t≥X (possible since RN(t)=o(t−N) and N≥2). For x≥X,
∫x∞RN(t)dt≤ε∫x∞t−Ndt=N−1εx−(N−1),
so ∫x∞RN=o(x−(N−1)). Since N≥2 was arbitrary, every order of the claimed expansion is verified.
For the last sentence of part 4, suppose f is continuously differentiable for large x and f′∼∑n≥0dnx−n.
d0=0: by step 1 applied to f′, d0=limx→∞f′(x). If d0=0 then f(x)=f(x1)+∫x1xf′ grows like d0x, contradicting f(x)→a0 finite.
d1=0: with d0=0, step 1 at order 1 gives xf′(x)→d1. If d1=0 then f′(t)∼d1/t, so f(x)−f(x1)=∫x1xf′ grows like d1lnx, again contradicting convergence of f.
Hence the expansion of f′ starts at order 2, and part 3 applies to f′. Since f(X)→a0, ∫x∞f′(t)dt=a0−f(x), so
a0−f(x)∼n≥2∑n−1dnx−(n−1),that isf(x)∼a0−m≥1∑mdm+1x−m.
By the uniqueness of part 1, am=−dm+1/m for m≥1, that is dm+1=−mam: the expansion of f′ is the term-by-term derivative of that of f.
The counterexample showing that the hypothesis "f′ has an expansion" cannot be dropped is in the remark below. □
Remark(Where the hypotheses are used, and the standing counterexamples)Original
"For each fixed N" is the whole content of (P). Every statement above is proved one order at a time, with N frozen before x is sent to infinity. Nothing is ever summed over N, and nothing may be: the series usually diverges.
Why differentiation fails. Let f(x)=e−xsin(ex). Since ∣f(x)∣≤e−x, part 2 of the uniqueness statement gives f∼∑n≥00⋅x−n, the zero expansion. But
f′(x)=−e−xsin(ex)+cos(ex),
which oscillates between −1 and 1 for arbitrarily large x and so does not tend to a limit. Thus f′ has no asymptotic power-series expansion at all, let alone the differentiated one. The hypothesis that fails is smallness of f′: an asymptotically negligible function may have a large derivative. This is why part 4 is stated conditionally.
Why a0=a1=0 is needed for integration. With f(x)=1/x we have f∼x−1, but ∫x∞dt/t diverges. The expansion says nothing about convergence of the integral; that has to be imposed.
What the expansion cannot see. Part 2 is not a defect to be repaired but a permanent feature. Exponentially small terms, such as e−x at +∞ or e−1/ε as ε↓0, are invisible to every order of a power-series expansion, and two physical answers that differ by such a term have identical expansions. Arfken's footnote records that Poincare's definition "allows (or neglects)" exactly these terms, and that sharpening it matters for the theory in the complex plane (Arfken, Weber & Harris (2013), footnote 4, p. 583). In the complex plane the same function may need different expansions in different sectors (Riley, Hobson & Bence (2006), §25.7.3, point (iii), p. 902); the switching is the Stokes phenomenon, and The WKB approximation meets it.
An expansion is not a Laurent series. Ahlfors warns against the confusion explicitly: the function whose expansion he computes is not even defined in a neighbourhood of infinity, so it has no Laurent development, and the formal series he writes does not converge (Ahlfors (1979), p. 205). The resemblance of the two notations is an accident of writing powers of 1/x.
Given f∼∑n≥0anx−n and a fixed numerical value of x:
Form the terms, not the sum. Compute un=anx−n for n=0,1,2,… and tabulate ∣un∣. For an expansion with an growing like n! the ratio ∣un+1/un∣ grows linearly in n, so the terms fall, reach a minimum, and then rise without limit.
Find the smallest term. Let n∗(x) be the index minimising ∣un∣. Riley, Hobson and Bence make the point that n∗ depends on x: there is no "number of terms" attached to the series, only one attached to the series and the value of x.
Stop just before it. Take sn∗−1(x) as the numerical value. Adding further terms makes the answer worse, as Arfken's table shows for exE1(x) at x=5.
Quote the error as the smallest term. The achievable accuracy is of the order of ∣un∗∣; for the expansions of this chapter that is roughly e−cx for some c>0, so it improves rapidly with x and cannot be improved at all at fixed x.
Turn the rule into a theorem when you can. The step above is a rule of thumb unless the remainder has been bounded. When the series alternates and one can prove ∣RN∣≤∣uN+1∣ — the first omitted term — consecutive partial sums bracket the answer, and steps 3 and 4 become rigorous: the value lies between sn∗−1 and sn∗, and their half-difference is an honest error bar. This happens whenever the remainder is an integral of a function of one sign, as in the example below.
If that is not accurate enough, change the problem, not the truncation. Increase x (for a sum, by evaluating more terms explicitly first and expanding only the tail; see Euler-Maclaurin summation), or find an expansion in a different variable. Never add more terms.
Example(The Gaussian tail at three standard deviations)Original
Let
Q(x)=∫x∞e−t2/2dt(x>0),
the unnormalised upper tail of a standard Gaussian. We build its asymptotic expansion, with a rigorous remainder, and then run the optimal-truncation rule at x=3.
Step 1: integrate by parts, keeping the remainder exact. Write e−t2/2=t−1⋅te−t2/2 and note dtd(−e−t2/2)=te−t2/2. Then
The right-hand side is exactly the magnitude of the first omitted term. Two consequences: RN(x)=O(x−2N−1e−x2/2), so the series is an asymptotic expansion of ex2/2Q(x) in the sense of (P) (in the asymptotic sequence φn=x−2n−1); and, the remainder alternating in sign with N, consecutive partial sums bracket the answer, so step 5 of the method applies.
Step 3: run the truncation rule at x=3. The terms are un=(2n−1)!!/32n+1, with ∣un+1/un∣=(2n+1)/9. This is below 1 for n≤3 and equals 1 at n=4, so the terms fall until n=4, where u4 and u5 tie, and rise afterwards: n∗=4. Writing PN=e−9/2∑n<N(−1)nun for the partial sums and comparing with the value Q(3)=0.00338369257 obtained by quadrature:
N
first omitted term e−9/2uN
PN
PN−Q(3)
1
4.11×10−4
0.003702999
+3.19×10−4
2
1.37×10−4
0.003291555
−9.21×10−5
3
7.62×10−5
0.003428703
+4.50×10−5
4
5.93×10−5
0.003352509
−3.12×10−5
5
5.93×10−5
0.003411771
+2.81×10−5
6
7.24×10−5
0.003352509
−3.12×10−5
7
1.05×10−4
0.003424940
+4.12×10−5
8
1.74×10−4
0.003320318
−6.34×10−5
Every entry obeys the bracketing: the errors alternate in sign and are smaller in magnitude than the corresponding first omitted term. The best values are P4 and P5, and the rule of step 5 returns
Q(3)∈[0.0033525,0.0034118],midpoint 0.0033821,
against the true 0.0033837: a relative error of 4.6×10−4 from five terms of a series that diverges. From N=6 on the partial sums move away again, exactly as step 3 predicts.
Step 4: the same rule at x=6. Now ∣un+1/un∣=(2n+1)/36, so the terms fall until n∗=18 and the attainable accuracy is far higher. Truncating at N=18 gives Q(6)=2.473008486×10−9 with a relative error of 1.1×10−8 — about eight correct figures. The lesson of the two tables together is the one the definition insists on: the quality of an asymptotic expansion is a statement about x, not about the number of terms.
Application(Tail probabilities, and the expansions of this part of the handbook)Original
Exceedance probabilities. If X is standard normal, Pr(X>x)=Q(x)/2π. Engineering and climate work routinely needs this at x between 3 and 6 — a design wave height, a 5σ detection threshold, a once-in-a-century return level — where the convergent Taylor series ∫0xe−t2/2dt=∑n(−1)nx2n+1/(2nn!(2n+1)) is numerically hopeless — at x=6 its largest term is 1.06×106 and the answer it must produce after cancellation is 2.5×10−9, fifteen digits down — while the expansion above is immediate. The five-term estimate of the example gives Pr(X>3)≈0.0033821/2π=1.3493×10−3 against the true 1.3499×10−3, and the eighteen-term estimate at x=6 gives Pr(X>6)=9.86588×10−10, correct to the figures shown. The same series is the standard large-argument representation of the complementary error function; see Error function and incomplete Gamma function and Watson's lemma, where it is derived again in one line.
Where else these series come from. Almost every closed-form answer a physicist uses in a limiting regime is an asymptotic expansion in this sense and not a convergent one:
the large-argument forms of the Bessel, Airy and Hankel functions, obtained from their integral representations (Asymptotic forms of Bessel functions);
Stirling's series for Γ, whose terms eventually grow (Stirling's formula);
perturbation series in a small coupling, where divergence is generic rather than exceptional (Regular perturbation expansions).
In all of these the practical questions are the ones this unit answers: how many terms to keep, what error to quote, and what the series cannot tell you however many terms are computed.
Turns a sum over integers into an integral plus a short, explicitly bounded series of boundary derivative terms built from the Bernoulli numbers, so that a slowly convergent series can be evaluated to eight figures from a handful of terms.
both valid for 0<∣t∣<2π, the radius of convergence being set by the zeros of et−1 nearest the origin at t=±2πi. Setting x=0 in the second gives Bn(0)=Bn. The first few are
Three properties are all the proof below needs. Stein and Shakarchi show that they characterise the family.
(B1) Derivative.Bn′(x)=nBn−1(x) for n≥1.
(B2) Endpoints.Bn(1)=(−1)nBn(0); hence Bn(1)=Bn(0)=Bn for every even n, and Bn(1)=Bn(0)=Bn=0 for every odd n≥3. The exception is n=1: B1(0)=−21 and B1(1)=+21.
(B3) Mean zero.∫01Bn(x)dx=0 for n≥1. This follows from (B1) and (B2): the integral equals (Bn+1(1)−Bn+1(0))/(n+1), and n+1≥2.
For (B1), differentiate the second generating function with respect to x: the left side becomes t times itself, which shifts the index by one. For (B2), replace t by −t in the first generating function and use −t/(e−t−1)=tet/(et−1), which is the second one at x=1.
A warning about conventions. Ahlfors writes the Laurent development of (ez−1)−1 with coefficients Bk that are all positive and carry half our index: his Bk is (−1)k−1 times our B2k(Ahlfors (1979), Ch. 5, §1.3, Ex. 4, p. 186). Arfken notes that a third convention, Bn=2(2n)!ζ(2n)/(2π)2n with halved index and all signs positive, is also in use (Arfken, Weber & Harris (2013), Eq. (12.39), p. 564). Everything in this handbook uses the convention displayed above, in which B1=−1/2 and the signs of B2n alternate. Formulas lifted from a book must be checked against its own definition before being combined with ours. The arithmetic properties of these numbers, and the closed form of ζ(2n) they produce, are collected in Riemann zeta function and Bernoulli numbers.
the 1-periodic extension of Bk restricted to [0,1). By (B2), Bk(1)=Bk(0) for every k≥2, so Bk is continuous for k≥2; B1 is the sawtooth x−⌊x⌋−21, which jumps at the integers. On [0,1) we have Bk=Bk, so Bk(m)=Bk at every integer m.
Theorem(Euler-Maclaurin summation formula)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.3, Eqs. (12.56)-(12.57), pp. 568-569
Let a<b be integers, let q≥1, and let f∈C2q([a,b]). Then
Only even-index Bernoulli numbers occur: the odd ones from B3 on vanish, and the single surviving odd term, from B1, is the 21(f(a)+f(b)) written separately.
Half-open form. If f∈C2q([a,∞)), if ∫a∞f and ∫a∞∣f(2q)∣ converge, and if f(x)→0 and f(2p−1)(x)→0 as x→∞ for p=1,…,q, then letting b→∞ gives
Arfken writes the formula the other way round, as the integral over [0,n] expressed through the trapezoidal sum, and reads the p-sum as the correction to the trapezoidal rule. The two arrangements differ only by moving terms across the equals sign.
Idea of the proofOriginal
Integrating by parts turns ∫01f into boundary values minus ∫01f′ — and the freedom in choosing the antiderivative of 1 is exactly where the Bernoulli polynomials come in. Take the antiderivative B1(x)=x−21, chosen to be odd about the midpoint; its boundary values ∓21 assemble the trapezoidal average 21(f(0)+f(1)). Repeat: B1 has antiderivative B2/2, whose boundary values are equal by (B2), so the boundary term collapses to B2[f′(1)−f′(0)]/2. Repeat again: B2 has antiderivative B3/3, which vanishes at both ends, so that step produces no boundary term at all. That alternation — one boundary term, then none, then one, then none — is why only even-index Bernoulli numbers appear, and it lets the two steps be bundled into a single induction from q to q+1. Finally, the one-interval identity is applied to each of [a,a+1],…,[b−1,b] and the results added: the half-weights at interior nodes combine into full weights, the derivative terms telescope to the two endpoints, and the leftover integrals glue into one integral against the periodised Bernoulli function.
Let g∈C2q+2([0,1]) and transform the last integral in two integrations by parts. First, B2q=B2q+1′/(2q+1) by (B1), and B2q+1(1)=B2q+1(0)=0 by (B2) because 2q+1≥3 is odd; so the boundary term vanishes and
Substituting into (∗q) gives (∗q+1). With (∗1) from Step 1, (∗q) holds for every q≥1.
Step 3 (add the intervals). Fix q and f∈C2q([a,b]). For each integer ν with a≤ν≤b−1 apply (∗q) to g(x)=f(x+ν), whose derivatives are g(k)(x)=f(k)(x+ν), and sum over ν:
∑ν∫01f(x+ν)dx=∫abf, by the change of variable x↦x−ν on each piece.
∑ν21(f(ν)+f(ν+1))=∑j=abf(j)−21(f(a)+f(b)): each interior node is counted with weight 21 twice, the two end nodes once.
∑ν[f(2p−1)(ν+1)−f(2p−1)(ν)]=f(2p−1)(b)−f(2p−1)(a), a telescoping sum.
∑ν∫01B2q(x)f(2q)(x+ν)dx=∫abB2q(x)f(2q)(x)dx, because on [ν,ν+1) the periodised function B2q is B2q(x−ν).
Step 5 (the half-open form). Let b→∞ through the integers. The three limits assumed in the statement dispose of f(b) and of f(2p−1)(b), the integral ∫abf converges by assumption, and the remainder converges because ∫a∞∣f(2q)∣ does, with the same bound. □
Note(Why ∣B2q∣≤∣B2q∣)Original
Let k≥1 and let cj=∫01Bk(x)e−2πijxdx be the Fourier coefficients of the 1-periodic function Bk (see Fourier series). By (B3), c0=0. For k=1, integration by parts gives
cj=∫01(x−21)e−2πijxdx=−2πij1(j=0).
For k≥2, integrating by parts and using Bk′=kBk−1 together with Bk(1)=Bk(0) kills the boundary term and leaves cj(Bk)=2πijkcj(Bk−1). By induction,
cj(Bk)=−(2πij)kk!(j=0,k≥1).
For k≥2 these coefficients are O(∣j∣−2), so the Fourier series converges absolutely and uniformly; its sum is continuous and has the same coefficients as the continuous function Bk, so the two agree:
Bk(x)=−k!j=0∑(2πij)ke2πijx(k≥2).
Taking absolute values with k=2q and using ∣(2πij)2q∣=(2π∣j∣)2q,
the last equality being Arfken's closed form B2q=(−1)q+12(2q)!ζ(2q)/(2π)2q(Arfken, Weber & Harris (2013), Eq. (12.38), p. 564). The bound is attained at x=0, where B2q=B2q, so it is sharp.
Remark(Where the hypotheses are used, and why the formula is asymptotic)Original
Integer endpoints. Step 3 adds whole unit intervals. With non-integer a or b the telescoping fails and extra fractional-interval terms appear; the formula as stated is simply false, as f≡1 on [0,21] already shows.
f∈C2q, and nothing more. The proof integrates by parts 2q times and never differentiates beyond that order, so 2q continuous derivatives is exactly what is needed for the formula at that q. A function with few derivatives still gets a valid formula at small q.
The bound needs f(2q) integrable, not small. No smallness of the derivatives is assumed anywhere; the quality of the result is decided afterwards, by how large ∫∣f(2q)∣ turns out to be.
Increasing q usually makes things worse eventually. Nothing in the theorem says Rq→0 as q→∞ at fixed a,b,f, and Arfken records that in applications the remainder integral, though small, generally does not tend to zero (Arfken, Weber & Harris (2013), p. 569). The reason is visible in the bound: ∣B2q∣/(2q)!=2ζ(2q)/(2π)2q shrinks like (2π)−2q, but ∫∣f(2q)∣ normally grows faster, like (2q)! divided by a fixed power. The q-expansion is therefore an asymptotic series in the strict sense, and the optimal-truncation rule of that unit applies: stop at the smallest correction term.
When it is exact. If f is a polynomial of degree d, then f(2q)≡0 for 2q>d and the formula is an identity with no remainder. Taking f(x)=x2, a=0, b=N, q=1 gives ∑j=0Nj2=N3/3+N2/2+N/6, which is N(N+1)(2N+1)/6: the classical closed forms for power sums are the degenerate case.
What goes wrong with a rapidly oscillating summand. Take f(x)=cos(ωx) on [a,b]. Then ∣f(2p−1)∣≤ω2p−1 and ∣B2p∣/(2p)!=2ζ(2p)/(2π)2p, so a correction term is of size at most (2/ω)(ω/2π)2p and the remainder bound is of size (b−a)ω(ω/2π)2q at worst. For ω<2π both shrink geometrically and the formula is excellent; for ω>2π both grow with p and q from the first term on, and no truncation is useful. The threshold ω=2π is exactly the radius of convergence of the generating function that defines the Bn, and it is the statement that the summand must be resolved by the integer grid: one sample per period is not enough. Euler-Maclaurin is a tool for summands that vary slowly from one integer to the next.
Method(Summing a slowly convergent series)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.3, Example 12.3.1 and Table 12.4, pp. 569-570
To evaluate S=∑n=1∞f(n) for a smooth, eventually monotone f decaying too slowly for direct summation:
Split the series. Choose a cut N and write S=∑n=1N−1f(n)+∑n=N∞f(n). The first part is evaluated exactly, term by term.
Check the limits. Confirm that ∫N∞f converges and that f and the odd derivatives f′,f′′′,…,f(2q−1) tend to 0 at infinity, so the half-open form applies.
Apply the formula to the tail only. Add ∫N∞f+21f(N) and then the corrections −B2pf(2p−1)(N)/(2p)! for p=1,2,…
Bound the error before deciding where to stop. Compute (2q)!∣B2q∣∫N∞∣f(2q)∣. If f(2q) has constant sign on [N,∞) — the usual case for a monotone decaying f — this integral is just ∣f(2q−1)(N)∣, so the bound equals the magnitude of the last correction term retained. The error is then never worse than the last term kept.
Choose N and q together. Raising N costs one evaluation of f per step and shrinks every correction; raising q costs one derivative and helps only while the terms are still falling. Arfken's table makes the first point vividly: moving the cut from N=1 to N=4 took his estimate of ζ(3) from 1.1667 to 1.202057 against a true value of 1.202057, that is from one correct figure to seven.
Stop when the terms stop falling. The corrections form an asymptotic series in 1/N; truncate at the smallest one, as in Asymptotic series.
Example(ζ(3/2) to eight figures from four terms and three corrections)Original
The series ∑n≥1n−3/2 converges, but hopelessly slowly: its tail beyond N is about 2N−1/2, so eight-figure accuracy by direct summation would need some 4×1016 terms. Euler-Maclaurin does it with seven numbers.
Set-up. Take f(x)=x−s with s=23, and cut at N=5. All the derivatives are available in closed form:
f(k)(x)=(−1)k(s)kx−s−k,(s)k:=s(s+1)⋯(s+k−1),
so f and every derivative tend to 0 at infinity and ∫5∞x−3/2dx=2⋅5−1/2 converges: step 2 of the method is satisfied for every q. The half-open formula gives, for the tail,
so ∣Rq∣≤(2q)!∣B2q∣(s)2q−1N−s−2q+1, which is exactly the magnitude of the p=q correction term. This is the situation described in step 4: the error is bounded by the last term kept.
The arithmetic. With s=3/2, N=5, and B2=1/6, B4=−1/30, B6=1/42:
contribution
value
∑n=14n−3/2
1.67100348032
∫5∞x−3/2dx=2/5
0.89442719100
21⋅5−3/2
0.04472135955
p=1: 2!B2s5−5/2=12⋅55/2s
+0.00223606798
p=2: −720⋅59/2(s)3
−0.00001304373
p=3: +30240⋅513/2(s)5
+0.00000030746
Running totals: after the first three rows, 2.61015203087; after p=1, 2.61238809885; after p=2, 2.61237505512; after p=3,
ζ(3/2)≈2.61237536258,∣R3∣≤3.1×10−7.
The true value is ζ(3/2)=2.61237534869, so the actual error is +1.4×10−8, comfortably inside the guaranteed bound and about an order of magnitude better than it. The corrections are still falling at p=3 — the ratio of successive ones is about 1/40 here — so a fourth would help; that is the optimal-truncation test of step 6, applied in the direction of "keep going".
Reading the table. The whole calculation used four exact terms of the series, one elementary integral, and three derivative corrections at the single point x=5. Against 4×1016 terms of direct summation, that is the point of the formula.
Application(Bose-Einstein condensation, lattice sums and Romberg extrapolation)Original
The condensation temperature of an ideal Bose gas. For a non-interacting gas of bosons of mass m in a box, the number of particles in excited states at temperature T and zero chemical potential is obtained by summing the Bose occupation numbers over momentum states. Replacing the sum by an integral — the leading Euler-Maclaurin term, legitimate because the level spacing is tiny compared with kBT — gives the density
nexc=λT3gζ(3/2),λT=2πmkBTh,
with g the spin degeneracy. The gas cannot hold more than this in excited states, so condensation begins when nλT3=gζ(3/2), that is at
Tc=mkB2πℏ2(gζ(3/2)n)2/3.
The number computed in the example is the one that appears here. For spin-polarised 87Rb (m=1.4432×10−25kg, g=1) at a typical trapped density n=1×1014cm−3=1×1020m−3,
Tc=3.98×10−7K≈398nK,
in the few-hundred-nanokelvin range where these experiments are done. The corrections in the example matter too: for a gas in a finite box the discreteness of the levels is exactly the difference between the sum and the integral, and the 21f(N) and derivative terms are the leading finite-size shifts of Tc.
Lattice sums. The same split — exact summation over a core of near neighbours, Euler-Maclaurin for the tail — is the standard way to evaluate slowly convergent lattice sums in solid-state physics, such as Madelung constants and dispersion-energy sums ∑∣R∣−6 over a Bravais lattice, where direct summation converges far too slowly to be useful.
Quadrature. Read in Arfken's direction, the formula says that the trapezoidal rule's error on a smooth integrand is a series in even powers of the step h with coefficients built from boundary derivatives (Arfken, Weber & Harris (2013), Eq. (12.57), p. 569). Two consequences are standard in numerical analysis: the error expansion in h2 is what Richardson extrapolation eliminates order by order, giving the Romberg scheme (Numerical quadrature); and for a periodic integrand integrated over a whole period, every boundary derivative term cancels, so the trapezoidal rule converges faster than any power of h — the reason spectral methods use equally spaced quadrature on periodic domains (Spectral methods).
Stirling's series. Applying the formula to f(x)=lnx turns lnn!=∑j=1nlnj into nlnn−n plus boundary terms and an asymptotic series in 1/n with Bernoulli coefficients; see Stirling's formula, where the constant 21ln(2π) that Euler-Maclaurin leaves undetermined is pinned down by other means.
Shows that the large-parameter behaviour of a Laplace-type integral is fixed entirely by the behaviour of its integrand at the endpoint, so that expanding the integrand near the origin and replacing each power by a Gamma function produces, term for term, a valid asymptotic expansion.
with x a large positive parameter and f measurable and locally integrable on (0,∞). Arfken and Riley, Hobson and Bence call F the Laplace transform of f; see The Laplace transform.
The function f is of exponential order if there are constants M≥0, c∈R and T>0 with
∣f(t)∣≤Mect(t≥T).
Arfken states this condition in the equivalent form ∣e−ctf(t)∣≤M and notes that it guarantees convergence for x>c; his standing counterexample is f(t)=et2, which satisfies no such bound and has no transform at all. Riley, Hobson and Bence put the same fact as the existence of an abscissa s0 with convergence for x>s0 and divergence below it.
which is the substitution u=xt in Euler's integral for Γ (The Gamma function). Arfken writes it as L{tn}=Γ(n+1)/sn+1 and records that the integral diverges at the origin once the exponent reaches −1; that is the restriction μ>0.
Theorem(Watson's lemma)Original
Let λ>0, let f be measurable and locally integrable on (0,∞), and suppose:
(W1) Exponential order. There are M≥0, c∈R and T>0 with ∣f(t)∣≤Mect for t≥T.
(W2) An asymptotic expansion at the origin. There are constants a0,a1,a2,… and a δ∈(0,T] such that for every N≥1 there is a constant CN with
f(t)−n=0∑N−1antn+λ−1≤CNtN+λ−1(0<t≤δ).
Then F(x)=∫0∞e−xtf(t)dt converges absolutely for every x>c, and
F(x)∼n=0∑∞anxn+λΓ(n+λ)(x→+∞),
in the sense of Asymptotic series: for every N≥1 there are KN and xN with
F(x)−n=0∑N−1anxn+λΓ(n+λ)≤xN+λKN(x≥xN).
In words: expand the integrand near the endpoint and integrate the expansion term by term against e−xt over the whole half-line. No convergence of ∑nantn+λ−1 is assumed, and none is used.
Idea of the proofOriginal
Because e−xt decays on the scale t∼1/x, the integral is deaf to everything except an ever-shrinking neighbourhood of t=0. Two estimates make this precise and together they are the whole proof. First, whatever happens beyond a fixed δ>0 contributes at most O(e−δx), which is smaller than every power of 1/x and so cannot disturb any order of the expansion. Second, on (0,δ] the integrand may be replaced by the first N terms of its expansion at the cost of a remainder bounded by CNtN+λ−1, and that remainder integrates against e−xt to exactly CNΓ(N+λ)x−(N+λ) — one order beyond the last term kept. The individual powers then integrate to Gamma functions, up to another exponentially small correction from extending the range from δ to ∞. Nothing is ever summed over n; the argument is run once for each fixed N, which is precisely what the Poincare definition asks for.
ProofOriginal
Throughout fix N≥1, write μn=n+λ, and assume x≥x0:=max(c+1,2).
Step 0 (convergence). Taking N=1 in (W2) gives ∣f(t)∣≤∣a0∣tλ−1+C1tλ on (0,δ], and λ>0 makes both powers integrable at the origin; f is locally integrable on [δ,T]; and (W1) gives ∫T∞e−xt∣f∣≤M∫T∞e−(x−c)tdt<∞ for x>c. So F(x) converges absolutely for x>c.
Step 1 (everything past δ is beyond all orders). Split at T and use δ≤T:
Since xKe−xδ→0 for every K, this term is O(x−K) for every K, and in particular O(x−(N+λ)).
Step 2 (the monomials, with the range extended). For μ>0 the substitution u=xt gives ∫0∞e−xttμ−1dt=Γ(μ)x−μ. The piece we must discard is small: for x≥2,
By Steps 1-3 the three error contributions are respectively O(x−K) for every K, at most CNΓ(N+λ)x−(N+λ), and O(x−K) for every K. Hence
F(x)−n=0∑N−1anxμnΓ(μn)≤xN+λKN
for all x beyond some xN, with KN built from CNΓ(N+λ) and the (exponentially small) constants of Steps 1 and 2. As N was arbitrary, this is the assertion. □
Remark(Where the hypotheses are used, and what each one rules out)Original
λ>0 is needed in Steps 2 and 3, where Γ(μ) is finite only for μ>0. If f(t)∼a0/t at the origin, F(x) does not exist for any x: Arfken makes the same point for integer powers, noting that ∫0∞e−sttndt diverges at the origin once n≤−1(Arfken, Weber & Harris (2013), p. 1008). The lemma never repairs a divergent integral; it describes one that converges.
(W1) is used only in Step 1, to make the far field exponentially small. Dropping it costs everything: for f(t)=et2 there is no x for which F(x) exists (Arfken, Weber & Harris (2013), p. 1008). Note how weak the hypothesis is, though: f may grow like ect with c as large as one likes, since the conclusion is a statement about x→∞.
(W2) is used only in Step 3. It is an asymptotic hypothesis, not a convergence hypothesis: the series ∑antn+λ−1 may diverge for every t>0, and Step 3 never adds up more than N of its terms. Replacing (W2) by "f is analytic near 0" would be a strictly stronger and quite unnecessary assumption.
Only the germ of f at 0 is seen. Alter f on [δ,∞), keeping (W1): by Step 1 the value of F changes by O(e−δx), so the whole expansion is unchanged. This is the non-uniqueness built into every asymptotic expansion (Asymptotic series): one expansion, many functions. A corollary worth stating is that no amount of work with the expansion can recover f away from the origin, or detect a term such as e−δx in F.
The direction matters. The lemma is Abelian: small-t information about f yields large-x information about F. The converse inference — large-x behaviour of F forcing small-t behaviour of f — is a Tauberian statement and is false without further hypotheses, again because F cannot distinguish f from f plus something supported away from the origin.
General exponents. Nothing in the proof uses that the exponents advance by 1. If (W2) holds in the form f(t)−∑n<Nantαn−1≤CNtαN−1 with 0<α0<α1<⋯→∞, the identical argument gives F(x)∼∑nanΓ(αn)x−αn. Expansions in powers of t1/2 or t1/3, which arise whenever a change of variable has a square-root or cube-root branch, are covered by this.
Method(Using the lemma)Original
Given an integral depending on a large parameter x:
Bring it to Laplace form. Substitute so that the exponent is −xt with t running over [0,∞) (or [0,δ], which by Step 1 of the proof is the same thing up to exponentially small terms), and so that the decay starts at t=0. If the exponent is xϕ(t) with ϕ decreasing and ϕ′(0)=0, the substitution u=ϕ(0)−ϕ(t) does it; if ϕ′ vanishes at the maximum, the right tool is Laplace's method instead.
Check (W1). Verify that the surviving f grows at most exponentially. Usually it is bounded or algebraically growing, and this is immediate.
Expand f about t=0. Pull out the leading power: write f(t)=tλ−1h(t) with λ>0 and h smooth (or analytic) at 0, and expand h in a Taylor or binomial series. Then an=h(n)(0)/n!, and Taylor's theorem with remainder (Taylor's theorem) supplies the constants CN of (W2) at once.
Substitute Gamma functions. Replace each antn+λ−1 by anΓ(n+λ)x−(n+λ). Nothing else changes.
Restore any prefactor removed in step 1, such as exϕ(0) or a Jacobian constant.
Truncate optimally. The result is an asymptotic series and normally divergent; stop at the smallest term and quote that term as the error (Asymptotic series).
Example(An integral with a half-integer exponent)Original
Take
F(x)=∫0∞t(1+t)e−xtdt(x>0).
Step 1: hypotheses. Here f(t)=t−1/2(1+t)−1/2 is positive, locally integrable on (0,∞), and bounded by 1 for t≥1, so (W1) holds with M=1, c=0, T=1. The integral converges for every x>0.
Step 2: expand at the origin. With λ=21, write f(t)=tλ−1h(t) with h(t)=(1+t)−1/2. The binomial series gives
h(t)=n≥0∑(n−1/2)tn,(n−1/2)=2nn!(−1)n(2n−1)!!,
convergent for ∣t∣<1; truncating it is Taylor's theorem with remainder, which supplies (W2) on (0,21] with explicit CN. So an=(n−1/2).
Step 3: Gamma functions. Using Γ(n+21)=(2n−1)!!π/2n,
using sinh2v=(cosh2v−1)/2, the substitution w=2v, and the integral representation of the modified Bessel function (Modified Bessel functions). At x=8 this gives
F(8)=e4K0(4)=0.609297669257,
which agrees with direct numerical quadrature to all figures shown.
Step 5: the series at x=8. The terms are un=π/8[(2n−1)!!]2/(4nn!8n), with ratio ∣un+1/un∣=(2n+1)2/(32(n+1)); this is below 1 until n=8, so the smallest term is u8 and the optimal truncation keeps n=0,…,7.
N
∑n<N(−1)nun
error
1
0.626657069
+1.74×10−2
2
0.607074035
−2.22×10−3
3
0.609827899
+5.30×10−4
4
0.609110747
−1.87×10−4
5
0.609385282
+8.76×10−5
6
0.609246299
−5.14×10−5
7
0.609333887
+3.62×10−5
8
0.609267805
−2.99×10−5
9
0.609325885
+2.82×10−5
10
0.609267603
−3.01×10−5
The errors alternate in sign and shrink until N=9, then start to grow: the signature of a divergent asymptotic series. Averaging the two best partial sums,
21(0.609267805+0.609325885)=0.609296845,
against the true 0.609297669: a relative error of 1.4×10−6 from nine terms of a series that converges for no value of x.
Application(Large-argument expansions of special functions, and early-time response)Original
The complementary error function. Watson's lemma is the standard machine for turning an integral representation into a large-argument expansion. Take erfc(x)=π2∫x∞e−s2ds for x>0 and substitute s2=x2(1+u), so that sds=21x2du and s=x1+u:
This is a Laplace-type integral with large parameter X=x2, with λ=1 and an=(n−1/2). Watson's lemma replaces un by n!X−n−1, and restoring the prefactor,
The coefficients are (−1)n(2n−1)!!/2n, which is exactly the series Arfken records for 1−erf(x)(Arfken, Weber & Harris (2013), Exercise 12.6.3, p. 584) — a check on the whole calculation. Numerically at x=3, where erfc(3)=2.20904970×10−5: the terms of the bracket fall until n=9 (the ratio of consecutive ones is (2n+1)/(2x2)=(2n+1)/18), and truncating after six of them gives 2.208627×10−5, a relative error of 1.9×10−4. See Error function and incomplete Gamma function.
The same route, with the same three lines, produces the large-argument expansions of the incomplete Gamma function, of the exponential integral, and — after the substitution of the example above — of K0 and the other modified Bessel functions (Asymptotic forms of Bessel functions).
Early-time response from a transform. In linear-response problems one often knows a material function only through its short-time behaviour. If the creep compliance of a viscoelastic solid behaves as J(t)=J0+J1t1/2+J2t+O(t3/2) as t↓0 — the half-power being the signature of diffusive relaxation — then Watson's lemma applied to Jˉ(s)=∫0∞e−stJ(t)dt gives, with no further information about the material,
sJˉ(s)=J0+2πs1/2J1+sJ2+O(s−3/2)(s→∞),
using Γ(3/2)=π/2. Since the high-frequency modulus is 1/(sJˉ(s)) at s=iω, the measured high-frequency behaviour determines J0,J1,J2 and conversely. The same correspondence — small t in the solution, large s in the transform — is what makes Laplace-transform solutions of diffusion problems (The Laplace transform, Causal Green's functions for initial-value problems) easy to read at early times and hard to read at late ones; the late-time end is the Tauberian direction that the remark above warns about.
Gives the leading large-parameter behaviour of an integral of the form ∫gexf whose exponent has a single non-degenerate interior maximum, by showing that the integral is controlled by a shrinking neighbourhood of that maximum on which the exponent is a Gaussian.
with f:(a,b)→R real-valued, g:(a,b)→C, and x→+∞. The exponent is real, so exf is positive and the large-x competition is purely one of size: the integrand is largest where f is largest, and raising x sharpens that preference without moving it.
A point t0∈(a,b) is an interior non-degenerate maximum of f if f is twice continuously differentiable near t0 with
f′(t0)=0,f′′(t0)<0.
We write A=∣f′′(t0)∣=−f′′(t0)>0 for the curvature of the exponent at the top.
Evans states the method for the equivalent form e−k/ε with ε↓0, which is the same thing with x=1/ε and f=−k. The name under which the method travels is not uniform. Hassani develops it in the complex plane, where method of steepest descent means the deformation of a contour so that it crosses a saddle point of f along the direction in which Ref falls fastest. Grimmett and Stirzaker attach the same name, "Laplace's method of steepest descent", to the real-axis statement itself — their Exercise 5.9.6 is about ∫0∞efn(x)dx for a real differentiable fn, with no contour anywhere in sight. So the two names overlap, and the real case treated here is the one in which no contour has to be moved; The method of steepest descents takes up the complex case.
Method(Laplace's method, leading order)AdaptedAdapted from Hassani (2013), §12.5, Eq. (12.34), p. 385, and Eq. (12.36), p. 387; Grimmett & Stirzaker (2020), Exercise 5.9.6, p. 216
Let t0∈(a,b) and suppose:
(L1)f is of class C2 on a neighbourhood of t0, with f′(t0)=0 and A:=−f′′(t0)>0;
(L2)g is continuous at t0, and g(t0)=0;
(L3)t0 is the strict global maximum of f with a gap: for every δ>0,
sup{f(t):t∈(a,b),∣t−t0∣≥δ}<f(t0);
This is Hassani's steepest-descent formula (12.34) in the case of a real maximum on the real axis, with ∣f′′(z0)∣ our A and his phase θ1 equal to 0. That phase records the direction in which the contour leaves the saddle: Hassani writes 21f′′(z0)=r2eiθ2 and fixes θ1 by 2θ1+θ2∈{π,3π} together with the convention 0≤θ1<π. Here f′′(t0)=−A<0, so θ2=π and the two candidate directions are θ1=0 and θ1=π; the convention excludes the second, which is why θ1=0 — at a real maximum the descent direction is along the real axis, taken in the direction of increasing t. The formula is also the precise form of the replacement Grimmett and Stirzaker state "under mild conditions", conditions (L1)-(L4) being one sufficient set.
The correction series. When f and g are smooth near t0, the leading term is the first of an asymptotic expansion in integer powers of 1/x,
which is Hassani's series in half-integer powers of the large parameter, specialised to the real axis. The pointer there should be to his corrected form (12.36), p. 387: Hassani first writes the expansion as (12.33), p. 385, then observes on p. 386 that "that infinite sum is in reality inconsistent", because only the first power of t was kept in the expansion of z(t) inside g(z)dz, and restores consistency by expanding z(t) as well. The powers α−k−1/2 with coefficients Γ(k+21) are the same in both, so the structure claimed here is unaffected; only the coefficients change. The remark below says where the expansion comes from.
Idea of the proofOriginal
Write the integrand as exf(t0) times ex(f(t)−f(t0)). The second factor is 1 at t0 and strictly less than 1 elsewhere, and raising x pushes it towards zero everywhere except at t0. Two things follow. Away from t0 the integrand is smaller than its peak value by a factor e−xη for some fixed η>0, so that region contributes exponentially little. Near t0 Taylor's theorem says f(t)−f(t0)≈−21A(t−t0)2, so the integrand is a Gaussian of width 1/xA; as x grows the Gaussian narrows, the window over which g is sampled shrinks onto t0, and g may be replaced by the constant g(t0). What remains is g(t0)exf(t0)∫e−xAu2/2du=g(t0)exf(t0)2π/(xA). The proof is this picture with ε's: a tolerance ε is fixed, a δ is chosen within which g varies by less than ε and f′′ by less than ε, the integral is squeezed between two Gaussians of curvature A±ε, and ε is then sent to zero afterx.
ProofOriginal
Set J(x)=x1/2e−xf(t0)I(x); the claim is J(x)→g(t0)2π/A.
Step 1 (choose δ for a given tolerance). Let 0<ε<A/2. By (L1) f′′ is continuous at t0, so there is δ1>0 with ∣f′′(τ)+A∣≤ε for ∣τ−t0∣≤δ1. For such t, Taylor's theorem with the Lagrange remainder (Taylor's theorem) and f′(t0)=0 give f(t)−f(t0)=21f′′(τ)(t−t0)2 for some τ between t0 and t, whence
By (L2) there is δ2>0 with ∣g(t)−g(t0)∣≤ε for ∣t−t0∣≤δ2. Put δ=min(δ1,δ2), small enough that [t0−δ,t0+δ]⊂(a,b).
Step 2 (the outside is exponentially small). By (L3) there is η>0 with f(t)≤f(t0)−η for t in the outer set O={t∈(a,b):∣t−t0∣≥δ}. For x≥x0, split the exponent as xf=x0f+(x−x0)f and bound the second factor on O:
The right-hand side is a fixed number depending only on ε (and on A, g(t0)), and it tends to 0 as ε↓0. Since the left-hand side does not depend on ε at all, it is zero. That is J(x)→g(t0)2π/A, which is the assertion. □
Remark(Where the hypotheses are used, and what each failure does to the answer)Original
Each hypothesis controls one feature of the answer, and dropping it changes the answer in a specific, computable way rather than merely breaking the proof.
(L1), f′′(t0)=0, fixes the power of x. If f′′(t0)=0 and the first non-vanishing derivative is of order 2m, the local shape is f(t0)−γ(t−t0)2m and the width of the peak is x−1/(2m), not x−1/2. For instance
∫−∞∞e−xt4dt=21Γ(41)x−1/4,
obtained from u=xt4. The general degenerate leading term is g(t0)exf(t0)m1Γ(2m1)(γx)−1/(2m), with γ=∣f(2m)(t0)∣/(2m)!; the Gaussian is replaced by a different one-dimensional integral, and the Gamma function appears for the same reason as in Watson's lemma.
An interior maximum, as opposed to an endpoint one, fixes the constant. If the maximum sits at t0=a with f′(a)=0, only half the Gaussian lies in the range and the leading term is halved. If instead f′(a)<0, so that f simply decreases from the endpoint, the peak is one-sided and exponential rather than Gaussian, and integration by parts (or Watson's lemma with λ=1) gives
∫abgexfdt∼x∣f′(a)∣g(a)exf(a),
with x−1 in place of x−1/2. Endpoint and interior contributions must both be examined before one is discarded.
(L2), g(t0)=0, is what makes the leading term leading. If g(t0)=0 the formula gives 0, correctly but uselessly: the proof then yields only I(x)=o(exf(t0)x−1/2), and one must go to the next term of the correction series. If g vanishes to order k at t0, which Gaussian moment survives decides the order. For evenk the moment ∫uke−xAu2/2du is a positive multiple of x−(k+1)/2 and that is the true order of e−xf(t0)I(x); for oddk that moment vanishes by symmetry, the next term of the expansion takes over, and the order drops to x−(k+2)/2. It can be smaller still — zero — when the symmetry is exact: for f(t)=−t2/2 on R and g odd, every moment vanishes and I(x)=0 identically.
(L3) is global, and no local hypothesis can replace it. The condition compares f(t0) with f everywhere else, with a gap. Without it a distant region can dominate: f(t)=−t2 on (−∞,∞) has a fine maximum at 0, but take instead a continuous f that equals −t2 on (−∞,10] and then climbs, with f(t)→+∞ as t→+∞. It still has a local maximum at 0 satisfying (L1), and (L2) still holds there with g≡1, but it has no global maximum, so (L3) fails — and so does (L4): I(x)=+∞ for every x>0. If two separated points share the maximum value, each contributes its own Gaussian and the two contributions are added — with their relative phase, if g is complex.
(L4) is the hypothesis that the integral exists at all, for at least one x; once it holds at x0 it holds for every x≥x0, since f≤f(t0).
Where the correction series comes from. Substitute s=f(t0)−f(t) separately on each side of t0. Near t0 this is a two-to-one change of variable with t−t0=±2s/A(1+O(s)), and it converts the integral into
I(x)=exf(t0)∫0σe−xsΦ(s)ds+(exponentially small),Φ(s)=∣f′(t+(s))∣g(t+(s))+∣f′(t−(s))∣g(t−(s)),
where t±(s) are the two solutions. Each branch on its own contributes every power s−1/2,s0,s1/2,s1,…. But Φ is the sum over the two branches, and t+−t0 and t−−t0 differ precisely in the sign of s, so when f and g are smooth the terms carrying an odd power of s — those coming from the odd Taylor coefficients of g and of the map s↦t — cancel between the branches and only the powers sn−1/2, n≥0, survive, beginning with 2/Ag(t0)s−1/2. That is exactly the shape (W2) of Watson's lemma requires with λ=21, and the lemma then returns a series in x−n−1/2 — integer powers of 1/x relative to the leading x−1/2. This is why Watson's lemma is listed as the prerequisite of this unit: Laplace's method is Watson's lemma plus a change of variable.
Note(Evans's ratio form: no smoothness needed)AdaptedAdapted from Evans (2010), §4.5.2, Lemma (Asymptotics) and its proof, p. 206
Many applications want not the integral itself but a weighted average against the same exponential, and for those much less is needed. Evans states and proves the following. Let k,l:R→R be continuous, with l growing at most linearly and k growing at least quadratically, and suppose k attains its minimum over R at a unique point y0. Then
There is no hypothesis of differentiability, let alone of non-degeneracy, on k — the price being that nothing is claimed about the size of either integral, only about their ratio. Evans's proof is worth remembering as the probabilistic reading of this whole unit: the normalised integrand
με(y)=∫e(k(y0)−k(z))/εdze(k(y0)−k(y))/ε
is a probability density on R which tends to 0 exponentially fast at every y=y0, so the family {με} concentrates all its mass at y0 and the average of any continuous l against it tends to l(y0). Laplace's method is concentration of measure; the statement proved above adds the extra information, bought with (L1), of how fast the mass concentrates and how large the integrals are.
Put the integral in the form ∫gexf with a single large parameter x multiplying a fixed real function f, and with any x-dependence of the prefactor removed. An integrand such as txh(t) becomes h(t)exlnt: all of the x-dependence must be moved into the exponent before the exponent is maximised.
Locate the maximum. Solve f′(t)=0 and compare the values of f at the roots and at the endpoints. Keep the global maximum; if several points tie, treat each and add.
Check the three things that change the answer: that the maximum is interior (otherwise halve, or use the endpoint formula of the remark), that f′′(t0)<0 (otherwise the power of x changes), and that g(t0)=0.
Shift and expand. Put t=t0+u and expand the exponent, f(t0+u)=f(t0)−21Au2+61f′′′(t0)u3+⋯. This is the step Riley, Hobson and Bence carry out explicitly when they set x=n+y in the Gamma integral.
Extend the limits to ±∞ and integrate. The added region is exponentially small (Step 2 of the proof), and Hassani makes the same remark about extending the range. The leading result is
I(x)≈g(t0)exf(t0)xA2π,A=−f′′(t0).
If a correction is wanted, keep the cubic and quartic terms of the exponent and the linear and quadratic terms of g. The whole exponent is xf, so the factor left over beside the Gaussian carries the large parameter too:
exp{x(61f′′′u3+241f′′′′u4)}=1+6xf′′′u3+24xf′′′′u4+72x2(f′′′)2u6+⋯.
Integrate term by term against the Gaussian using ∫−∞∞u2me−Axu2/2du=(2m−1)!!(Ax)−m2π/(Ax). Odd powers integrate to zero, killing the u3 term; the surviving u4 and u6 terms contribute x⋅x−2 and x2⋅x−3 relative to the leading Gaussian, so both land at relative order 1/x and the first correction is of that order. (Dropping the x from the exponential here is the standard slip: it would make the u4 term look like a relative x−2 and lose half of the genuine first correction.)
Sanity-check the scaling. The answer must be exf(t0) times x−1/2 times something independent of x; if a stray power of x survives, some x-dependence was left outside the exponent at step 1.
This is the normalising constant of a symmetric Beta distribution and the continuous analogue of the central binomial coefficient; the closed form gives exact values to check against.
Step 1: all the x-dependence into the exponent.[t(1−t)]x=exf(t) with
f(t)=lnt+ln(1−t),g(t)≡1on (0,1).
Step 2: the maximum.f′(t)=1/t−1/(1−t)=0 gives t0=21, and
Step 3: the hypotheses.f is smooth on (0,1), strictly increasing on (0,21) and strictly decreasing on (21,1), and f(t)→−∞ at both ends, so t0=21 is the strict global maximum with a gap: (L3) holds. g≡1 gives (L2), and (L4) holds with x0=1 since ∫01t(1−t)dt=61<∞. (L1) is Step 2.
Step 4: the leading term.
I(x)∼1⋅e−2xln28x2π=4x1⋅21xπ.
Step 5: the first correction. With u=t−21 one has t(1−t)=41(1−4u2), so the exponent is exactly
f(t)−f(t0)=ln(1−4u2)=−4u2−8u4−364u6−⋯(∣u∣<21).
Keeping the quartic term as a factor, exln(1−4u2)=e−4xu2(1−8xu4+O(xu6,x2u8)), and with ∫−∞∞e−4xu2du=21π/x and ∫−∞∞u4e−4xu2du=64x23⋅21π/x,
I(x)=4x1⋅21xπ(1−8x3+O(x−2)).
Step 6: the numbers. Exact values from Γ(x+1)2/Γ(2x+2):
x
exact B(x+1,x+1)
leading term
relative error
with 1−8x3
relative error
10
2.57740196×10−7
2.67266808×10−7
+3.70%
2.57244303×10−7
−0.192%
25
1.55112326×10−16
1.57425526×10−16
+1.49%
1.55064143×10−16
−0.031%
100
5.49442959×10−62
5.51500370×10−62
+0.374%
5.49432243×10−62
−0.0020%
At x=10, for instance, the leading term alone is 2.67266808×10−7 against the exact 2.57740196×10−7. The relative error of the leading term falls as 0.375/x — exactly the coefficient 3/8 predicted in Step 5 — and including that one correction improves each row by roughly a further factor x, as an expansion in 1/x should.
Application(Vanishing viscosity, Bayesian asymptotics and the saddle point of statistical mechanics)Original
The vanishing-viscosity limit of Burgers' equation. The Cole-Hopf formula solves the viscous conservation law utε+uεuxε=εuxxε explicitly as a ratio of two integrals,
with h an antiderivative of the initial datum. This is precisely the ratio Evans's lemma evaluates, with l(y)=(x−y)/t growing linearly and K growing quadratically in y. Whenever y↦K(x,y,t) has a unique minimiser y(x,t) — which Evans shows holds for all but countably many x at each fixed t>0 — the limit is (x−y(x,t))/t. That expression is the Lax-Oleinik formula, and it picks out exactly the entropy solution — the unique admissible one — of the inviscid law (Evans (2010), §4.5.2, Example 2 and Eq. (26), pp. 205-206). Laplace's method is therefore the mechanism by which an artificial viscosity selects the physically correct shock (Conservation laws and shocks). Note which version of the method is needed: the exponent K is only continuous, not smooth, at the points where the minimiser jumps, so the ratio form of the note is doing the work and the non-degeneracy hypothesis (L1) would be false.
Bayesian asymptotics and the normal approximation. With n independent observations, a posterior density is proportional to π(θ)enℓn(θ), where ℓn is the average log-likelihood (Maximum likelihood). Applying the theorem with x=n, f=ℓn and g=π about the maximum-likelihood estimate θ^ gives the model evidence
∫π(θ)enℓn(θ)dθ∼π(θ^)enℓn(θ^)n∣ℓn′′(θ^)∣2π,
and, applied to numerator and denominator at once, the statement that the posterior itself is asymptotically Gaussian with variance 1/(n∣ℓn′′(θ^)∣) — the observed-information formula. Taking logarithms of the evidence and keeping the terms that grow with n produces the −21logn per parameter of the Bayesian information criterion. See Bayesian inference.
The saddle point in statistical mechanics. A partition function written as an integral over an extensive variable, Z(N)=∫eNs(m)dm with s an entropy per particle, is evaluated by this method with x=N: the free energy per particle is maxms(m) and the corrections are O(logN/N), which is the statement that the thermodynamic limit exists and that fluctuations about the most probable state are Gaussian of width N−1/2. Two separated maxima with equal s is a first-order phase transition, and the remark above says what then happens to the formula: the two contributions add. The large-deviation rate function of Large deviations and Cramer's theorem is the same object seen from probability theory, and the complex-plane version needed when the exponent is not real is The method of steepest descents.
Gives the large-parameter behaviour of an oscillatory integral ∫geiλf: each interior point where the phase is stationary and its second derivative is non-zero contributes a term of size λ−1/2 carrying a fixed ±π/4 twist, each endpoint contributes a term of size λ−1, and everywhere else the oscillation cancels to all orders.
Definition(Oscillatory integrals, stationary points, degeneracy)AdaptedAdapted from Evans (2010), §4.5.3, Eq. (40) and (52)-(53), pp. 209, 214
Fix a real-valued phasef and an amplitudeg, both defined on an interval [a,b]⊆R, and a large real parameter λ>0. The object of study is the oscillatory integral
I(λ)=∫abg(t)eiλf(t)dt.
A point c with
f′(c)=0
is a stationary point (or critical point) of the phase; it is non-degenerate when f′′(c)=0, and degenerate otherwise.
In n dimensions the same integral is written, as Evans (Evans (2010), §4.5.3, Eq. (40), p. 209) does, with a small parameter ε in place of a large one:
Iε=∫Rna(y)eiϕ(y)/εdy,
with a∈Cc∞(Rn), ϕ smooth and real. The two notations are related by ε=1/λ, and both appear below. The stationary set is now {Dϕ=0}, and a stationary point yk is non-degenerate when the Hessian D2ϕ(yk) is non-singular. For a non-singular real symmetric matrix A, its signature
sgnA=(number of positive eigenvalues)−(number of negative eigenvalues)
records how many of the n directions curve upwards and how many downwards; in one dimension sgnf′′(c) is +1 at a minimum of the phase and −1 at a maximum.
Two features distinguish this from Laplace's method. The integrand here has modulus ∣g∣ everywhere, so nothing is exponentially small and nothing is exponentially large: the whole interval looks equally important to the eye. What makes most of it irrelevant is cancellation. And the answer decays only like a power of λ, not exponentially. Evans (Evans (2010), §4.5.3b, p. 208) describes the method as the variant of Laplace's method got by replacing the −1 in the exponent by i, and the formulas below bear that out, with one extra ingredient: a fixed phase shift of ±π/4 that has no counterpart in the real case.
Method(Stationary phase)AdaptedAdapted from Evans (2010), §4.5.3, Lemma 1 p. 209, Lemma 2 (41) p. 210, Eq. (54) p. 215; Riley, Hobson & Bence (2006), §25.8.3, Eq. (25.73), p. 916
1. The model integral (Fresnel). For σ∈{+1,−1} and λ>0,
∫−∞∞eiσλu2/2du=λ2πeiσπ/4,
the integral being understood as the limit of ∫−RR as R→∞.
2. No stationary point. Let g∈Cc∞((a,b)) and let f be smooth with f′=0 on the support of g. Then for every m∈N,
∫abg(t)eiλf(t)dt=O(λ−m)(λ→∞).
3. One non-degenerate stationary point. Let f be smooth and real on [a,b], let g∈Cc∞((a,b)), and suppose f′ vanishes at exactly one point c of the support of g, with f′′(c)=0. Put σ=sgnf′′(c). Then
5. The general case. Let f and g be smooth on [a,b] with f real. If f has finitely many stationary points c1,…,cN, all in the open interval (a,b) and all non-degenerate, and if f′(a)=0 and f′(b)=0, then the leading behaviour is the sum of one term of part 3 for each ck, of size λ−1/2, plus the two endpoint terms of part 4, of size λ−1:
6. In n dimensions.Evans (Evans (2010), §4.5.3, Eq. (54), p. 215) states, for a∈Cc∞(Rn), ϕ smooth, Dϕ vanishing on the support of a only at y1,…,yN and every D2ϕ(yk) non-singular,
as ε→0. He presents this as an informal assembly of his Lemmas 1-3 by a partition of unity rather than as a theorem with a proof. Part 3 is the case n=1, ε=1/λ, in which sgnD2ϕ is ±1 and the signature factor becomes e±iπ/4.
Idea of the proofOriginal
Away from a stationary point the integrand turns through a full cycle in a t-interval of length about 2π/(λ∣f′∣), so neighbouring cycles cancel in pairs and only the mismatch between them survives; each integration by parts makes that mismatch smaller by a factor λ−1, and with no boundary to stop the process the contribution is smaller than every power of 1/λ. At a stationary point the phase is flat to first order, so the integrand stays nearly constant over an interval of length about λ−1/2 — the width at which 21λ∣f′′∣(t−c)2 first reaches order one — and that whole interval adds up in phase. This is why the answer is of size λ−1/2: one stationary interval of width λ−1/2, counted without cancellation.
Three steps turn this picture into a formula. Cut the integral into a piece near c and a piece away from it, and dispose of the second by the integration by parts just described. In the first piece change the variable so that the phase becomes exactly quadratic — possible because f′′(c)=0, which is the one-dimensional Morse lemma. What is left is the Fresnel integral, the oscillatory twin of the Gaussian, and it is evaluated by rotating the contour through π/4 onto a ray where it is a Gaussian. That rotation is the source of the eiσπ/4: it is not a property of f or g at all, only of the sign of the curvature of the phase.
At an endpoint the cancellation argument still works but the integration by parts leaves a boundary term, and that term, of size λ−1, is all that an endpoint ever contributes.
Proof(Part 1: the Fresnel integral)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.8.3, Eqs. (25.67)-(25.69), pp. 912-913
Take σ=+1; the case σ=−1 follows by conjugating. Rescaling u=2/λs reduces the claim to ∫−∞∞eis2ds=πeiπ/4.
Consider e−z2, which is entire, and the closed contour made of the segment [0,R] of the real axis, the circular arc z=Reiθ with 0≤θ≤π/4, and the segment of the ray argz=π/4 back to the origin. By Cauchy's theorem the integral around it vanishes.
On the arc, ∣e−z2∣=e−R2cos2θ, and cos2θ≥1−4θ/π for 0≤θ≤π/4 (both sides agree at the two ends and the left side is concave there), so
Letting R→∞ in the vanishing contour integral therefore gives ∫0∞e−z2dz along the real axis equal to the same integral along the ray. Parametrising the ray by z=seiπ/4, so that z2=is2 and dz=eiπ/4ds, and using ∫0∞e−z2dz=π/2 from Gaussian integrals,
2π=eiπ/4∫0∞e−is2ds.
Conjugating and doubling (the integrand is even) gives ∫−∞∞eis2ds=πeiπ/4, which is Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §25.8.3, Eq. (25.68), p. 913) identity; they read off from it the values C(∞)=S(∞)=21 of the Fresnel integrals. Restoring the scaling, ∫eiλu2/2du=2π/λeiπ/4. □
Proof(Parts 2 and 4: integration by parts)AdaptedAdapted from Evans (2010), §4.5.3, Lemma 1, p. 209
Both rest on the single identity
eiλf(t)=iλf′(t)1dtdeiλf(t),
legitimate exactly where f′=0 — this is the step that uses the absence of a stationary point, and the only one.
Part 4. With f′=0 on [a,b], insert the identity and integrate by parts once:
The new amplitude g1=(g/f′)′ is continuous, because f′ does not vanish on the compact interval and is therefore bounded away from zero there; so ∣∫abg1eiλfdt∣≤∫ab∣g1∣dt, a constant, and the remainder is already O(λ−1). To reach the O(λ−2) claimed, integrate by parts once more, which needs g1∈C1 and hence f,g∈C3:
Under the weaker hypothesis f,g∈C2 the same first step gives the boundary term plus o(λ−1), the remainder tending to zero by the Riemann-Lebesgue lemma applied to the continuous amplitude g1 rather than by a second integration by parts.
Part 2. Now g has compact support inside (a,b), so the boundary term vanishes identically and the same step gives
∫geiλfdt=−iλ1∫g1eiλfdt,g1=(f′g)′,
with g1 again smooth and compactly supported in the set where f′=0. Repeating m times produces (−iλ)−m∫gmeiλfdt with gm smooth and compactly supported, and the last integral is bounded independently of λ. This is Evans's (Evans (2010), §4.5.3, Lemma 1, p. 209) Lemma 1, which he states for a linear phase p⋅y with p=0 and proves by differentiating m times along the coordinate in which p does not vanish; the version here is the same computation with a variable f′ in the denominator. □
Proof(Part 3: the stationary contribution)AdaptedAdapted from Evans (2010), §4.5.3, Lemma 3(ii) pp. 212-214, Lemma 2 p. 210, Eq. (54) p. 215
Write σ=sgnf′′(c) and μ=∣f′′(c)∣>0.
Step 1: localise. Choose δ>0 so small that f′ has no zero other than c in [c−δ,c+δ] and f′′ keeps its sign there, and a cut-off χ∈Cc∞((c−δ,c+δ)) with χ≡1 on [c−δ/2,c+δ/2]. Split g=χg+(1−χ)g. The amplitude (1−χ)g is smooth, compactly supported in (a,b), and supported where f′=0, so by part 2 its integral is O(λ−m) for every m. Only ∫χgeiλf remains.
Step 2: make the phase exactly quadratic. Taylor's theorem with integral remainder gives, for t near c and using f′(c)=0,
with h smooth and h(c)=f′′(c). Shrinking δ if necessary, σh>0 on [c−δ,c+δ]. Define
u(t)=(t−c)σh(t).
Since σh is smooth and strictly positive, its square root is smooth, so u is smooth, u(c)=0 and u′(c)=μ=0. By the inverse function theorem u is a diffeomorphism of a neighbourhood of c onto a neighbourhood of 0, and by construction
f(t)=f(c)+2σu(t)2.
This is the one-dimensional case of the Morse lemma, Evans's (Evans (2010), §4.5.3, Lemma 3(ii), pp. 212-214) Lemma 3(ii), whose hypothesis detD2ϕ(0)=0 is here exactly f′′(c)=0; in one variable the induction over coordinates that he needs collapses to the explicit formula above. Changing the variable,
The first integral is part 1. For the second, the phase σu2/2 has derivative σu=0 on the support of 1−ζ; integrating by parts as in part 2, with no boundary term because (1−ζ)/u vanishes at infinity, and noting that each step improves the decay of the amplitude by two powers of u, gives O(λ−m) for every m.
Step 4: the remainder. The function G−G(0)ζ is smooth, compactly supported and vanishes at u=0, so G(u)−G(0)ζ(u)=uH(u) with H smooth and compactly supported. Using ueiσλu2/2=iλσdudeiσλu2/2 and integrating by parts (no boundary term, H having compact support),
∫uH(u)eiσλu2/2du=−iλσ∫H′(u)eiσλu2/2du.
The integral on the right is of the same shape as the one we started with, with amplitude H′ and phase σu2/2; it is O(1) at once, which already gives an absolute error O(λ−1) — a relative error O(λ−1/2) at this stage, the leading term being of size λ−1/2. Feeding it back through steps 3-4 once more shows it is in fact O(λ−1/2), so the remainder is O(λ−3/2) — the relative error O(λ−1) asserted in part 3 and in Evans's (Evans (2010), §4.5.3, Eq. (54), p. 215) formula (54).
Step 5: collect. Combining steps 1-4,
∫abgeiλfdt=eiλf(c)μg(c)λ2πeiσπ/4+O(λ−3/2),
which is part 3. Part 5 follows by applying step 1 with one cut-off per stationary point, part 4 to what is left over (whose support now reaches the endpoints, so the boundary term survives), and adding. □
Remark(Where the hypotheses are used)Original
f real is what makes this a different method. The whole argument uses ∣eiλf∣=1: there is no large or small factor to exploit, only cancellation. If f has an imaginary part the integrand has peaks and troughs, and the right tool is Laplace's method or, in the complex plane, The method of steepest descents. The two are not merely analogous: the steepest-descents formula contains the stationary-phase formula as the case where the saddle sits on a level line, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.3, p. 917) verify that their two formulas, (25.66) and (25.73), coincide.
f′=0 is used exactly once, in the identity eiλf=(iλf′)−1(eiλf)′, and everything in parts 2 and 4 follows from it. Where f′ vanishes the factor 1/f′ blows up, which is the analytic signature of the fact that a stationary point contributes more than every endpoint.
f′′(c)=0 fixes both the power and the phase. It is the hypothesis behind the change of variable in step 2 — without it σh is not smooth and u is not a diffeomorphism. If f′′(c)=0 but f′′′(c)=0, the phase is cubic and the stationary interval has width λ−1/3 instead of λ−1/2, so the contribution is of size λ−1/3 and is given by an Airy function rather than by a Fresnel integral; see Airy functions. The transition between the two, as two stationary points merge, is where a uniform approximation is needed and the formula of part 3 fails: it predicts a contribution that blows up as f′′(c)→0.
Smoothness of g controls the size of what is left over, not the leading term. For part 2 we need to differentiate the amplitude m times, so g∈Cm buys O(λ−m) and no more. A jump in g acts like an endpoint: g=χ[0,1] in part 4 produces λ−1 terms at 0 and 1 even though the phase is perfectly smooth there.
Compact support in (a,b) is what removes the endpoint terms in parts 2 and 3, and nothing else. Part 5 is the honest statement for an integral over a finite interval: the endpoints always contribute at order λ−1, one order below the stationary points, so they matter whenever two-term accuracy is wanted, and they are the entire answer when there is no stationary point at all.
Non-degeneracy is also what makes the count finite. If f′ vanished on a whole subinterval, that subinterval would contribute at order 1: the integral ∫01eiλ⋅0dt=1 does not decay. The hypothesis of finitely many non-degenerate stationary points in part 5 rules this out, and Evans (Evans (2010), §4.5.3, (52)-(53), p. 214) imposes it in the same two pieces, (52) and (53).
The error estimate is not uniform in the data. The constant hidden in O(λ−3/2) depends on bounds for f′′′,f′′′′ and for the derivatives of g near c, and on how far f′′(c) is from zero. Two stationary points a distance d apart may be treated separately only when λd2∣f′′∣≫1, which is the condition Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.3, pp. 918-919) express by saying that the first Cornu spiral must have wound up before the second begins.
Method(Applying the method)Original
Put the integral in the standard form∫abg(t)eiλf(t)dt with f real and λ the large parameter. A real integrand ∫gcos(λf) is the real part of this, so work with the complex form and take real parts at the end. If the large parameter sits inside the integrand in some other way, rescale the variable first until it appears as an overall factor in the exponent, as Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, p. 911) do with s=x1/2t for the Airy integral.
Check non-degeneracy by computing f′′ at each root. If some f′′(ck)=0, or if two roots are closer than λ−1/2, stop: the formula below does not apply and a uniform (Airy-type) approximation is needed.
Write down one term per stationary point,
λ∣f′′(ck)∣2πg(ck)exp[iλf(ck)+i4πsgnf′′(ck)],
remembering that the sign of the π/4 follows the sign of the curvature of the phase, not the sign of anything else.
Add the endpoint termsiλ1[f′geiλf]ab if the interval is finite and g does not vanish there, and one such term at every jump of g or of f′. They are one power of λ smaller than the stationary terms and are the whole answer when there is no stationary point.
Sanity-check the result. The leading term must have modulus ∣g(c)∣2π/(λ∣f′′(c)∣), which has the dimensions of g times the dimensions of t; if the original integral is real, the stationary terms must combine into something real, which for two stationary points related by a symmetry means their phases are negatives of each other. Finally compare against a numerical quadrature at one moderate λ: the stationary terms on their own carry a relative error falling like λ−1/2, and once the endpoint terms are added the relative error should fall like λ−1, as part 5 asserts.
Example(An integral with one interior stationary point and two endpoints)Original
Take
I(λ)=∫02exp[iλ(31t3−t)]dt,
so g≡1 and f(t)=31t3−t.
Stationary points.f′(t)=t2−1 vanishes at t=±1, and only t=c=1 lies in (0,2). There f(1)=31−1=−32 and f′′(1)=2tt=1=2, which is non-zero, so the point is non-degenerate and σ=+1.
The stationary term. By step 4 of the method,
S(λ)=2λ2πei(−2λ/3+π/4)=λπei(π/4−2λ/3).
The endpoint terms.f′(0)=−1, f(0)=0, f′(2)=3, f(2)=38−2=32, so
Numerical quadrature of the original integral, split at t=1 so that the oscillation is resolved on each side, gives
I(100)=−0.178593013−0.023923248i,
so the two-term prediction is correct to 1.7×10−4 in absolute terms and to 0.097% relative to ∣I(100)∣=0.18019. That is consistent with the predicted remainder O(λ−3/2)=O(10−3).
What each piece is worth. The stationary term alone misses by 7.9×10−3, forty-five times the error of the two-term answer: at λ=100 the endpoints are not a refinement one can skip. Their ratio to the stationary term is of order λ−1/2=0.1, exactly as the orders in part 5 predict, and at λ=104 it would be 0.01.
ApplicationOriginal
Group velocity of a dispersive wave train. A one-dimensional linear wave problem with dispersion relation ω(k) has solutions
u(x,t)=∫−∞∞A(k)ei(kx−ω(k)t)dk.
Fix a ray x/t=V and let t→∞. Then u=∫A(k)eitφ(k)dk with φ(k)=kV−ω(k), so the large parameter is the elapsed time and the stationary points are the wavenumbers with
ω′(k)=V.
Part 3 then says u∼2π/(t∣ω′′(kV)∣)A(kV)ei(kVx−ω(kV)t)∓iπ/4: an observer travelling at speed V sees, after a long time, only the wavenumber whose group velocity is V, with an amplitude decaying like t−1/2 because the energy originally present spreads over a region growing like t. The t−1/2, the selection rule ω′(k)=V and the π/4 are all read straight off the formula, and the degenerate case ω′′(kV)=0 — where the rule selects a band of wavenumbers rather than one, the decay slows to t−1/3 and the Airy function appears — is exactly the failure described in the remark above.
High-frequency limit of the wave equation, and geometric optics.Evans (Evans (2010), §4.5.3, pp. 207-216) uses stationary phase to justify what the geometric-optics ansatz uε=aεeipε/ε only suggests. Solving utt=Δu with rapidly oscillating data gε(x)=a(x)eip(x)/ε by the Fourier representation and applying his formula (54) to the resulting double integral, he finds that the stationary set is precisely the ray geometry: the stationary conditions reproduce the characteristics of the Hamilton-Jacobi equation pt±∣Dp∣=0 for the phase. The amplitude of the answer carries ∣detDy,z2ϕ∣−1/2, which on the stationary set equals ∏i=1n−1∣1−tκi∣−1/2, since his (56), Evans (Evans (2010), §4.5.3, Eq. (56), p. 216), gives detDy,z2ϕ+=(−1)n∏i=1n−1(1−tκi), where κ1,…,κn−1 are the principal curvatures of the surface on which the initial phase is constant — so the approximation blows up exactly at the focal times t=1/κi, where neighbouring rays cross and the stationary point becomes degenerate. Stationary phase thus produces both geometrical optics and its caustics, which the formal ansatz cannot see.
Diffraction and the Cornu spiral. In Fresnel diffraction the field at a point behind an edge is ∫eiπu2/2du over the unobstructed aperture, and the amplitude-phase diagram that Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.3, pp. 914-915) draw — the Cornu spiral — is the graph of the partial integrals. Their picture is the method made visible: the tightly wound ends of the spiral are the regions of rapid phase variation, which cancel, and the straight central stretch is the stationary region, which does not. The oscillations of intensity in the illuminated region and the smooth decay in the geometrical shadow are the two ends of that spiral.
Asymptotics of special functions. The integral representations of the Bessel and Airy functions are oscillatory integrals with a parameter, and stationary phase gives their large-argument behaviour directly: for Jn(x)=π1∫0πcos(nθ−xsinθ)dθ with x→∞ at fixed n, the phase xsinθ−nθ has (for x>n) exactly one stationary point in (0,π), at cosθ=n/x; writing the cosine as 21(ei(nθ−xsinθ)+e−i(nθ−xsinθ)) turns that one point into a conjugate pair of contributions, which combine into the familiar 2/(πx)cos(x−nπ/2−π/4) of Asymptotic forms of Bessel functions — the π/4 being the signature factor of this unit. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.3, worked example, pp. 919-920) carry out the companion calculation for Ai(x) with x<0, where the two real stationary points at ±−x give the oscillatory branch π−1/2(−x)−1/4sin(32(−x)3/2+4π), in contrast with the single complex saddle that governs x>0 and that belongs to The method of steepest descents.
Why path integrals produce classical mechanics. The stationary points of the action are the classical trajectories, and the weight eiS/ℏ is an oscillatory integrand with λ=1/ℏ. The formal application of part 6 to the path integral — one term per classical path, with amplitude ∣detδ2S∣−1/2 and a phase π/4 times the signature of the second variation — is the semiclassical (WKB) approximation, with the signature appearing as the Maslov index. The second variation and its degeneracies are the subject of The second variation and Legendre's condition and Conjugate points and Jacobi's condition, and the one-dimensional case is The WKB approximation.
Evaluates ∫Cg(z)eλf(z)dz for large λ by deforming the contour onto a path through a saddle point of f on which the imaginary part of f is frozen, so that the oscillation disappears and the integral becomes a real Gaussian; gives the leading-order formula g(z0)eλf(z0)eiθ2π/(λ∣f′′(z0)∣) and the rule that fixes the direction θ.
Definition(Saddle points, level lines, steepest paths, valleys)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.8.1, Eqs. (25.58)-(25.63), pp. 905-908; Arfken, Weber & Harris (2013), §12.7, Eqs. (12.100)-(12.103), pp. 586-587; Hassani (2013), §12.5, Eqs. (12.23)-(12.28), pp. 382-383
Let f and g be analytic in a region containing a contour C, let λ>0 be a large real parameter, and consider
I(λ)=∫Cg(z)eλf(z)dz.
Hassani (Hassani (2013), §12.5, p. 382) notes that a complex λ may always be made real and positive by absorbing its phase into f; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, p. 585) and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8, p. 905) instead absorb λ into f altogether and write the integrand as gef. We keep λ visible, since it is the quantity the expansion is in.
Write
f(z)=φ(x,y)+iψ(x,y),z=x+iy,
so that ∣eλf∣=eλφ and argeλf=λψ. Since f is analytic, φ and ψ are harmonic and conjugate, so:
The curves φ=const, on which the modulus of the integrand is constant, are the level lines. The curves ψ=const, on which its phase is constant, are the lines of steepest descent or ascent. Because ∇φ⊥∇ψ, the two families cross at right angles, and the direction in which ∣eλf∣ changes fastest is the direction in which the phase does not change at all. That coincidence is the whole point of the method.
φ has no interior maximum or minimum where f is analytic: at one, φxx and φyy would have the same sign, contradicting φxx+φyy=0. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, p. 586, Exercise 12.7.1) call the corresponding statement for ∣F∣ Jensen's theorem.
A point z0 with
f′(z0)=0
is a saddle point of f. There ∇φ=∇ψ=0: both the modulus and the phase of the integrand are stationary. Writing
f′′(z0)=Aeiα,A=∣f′′(z0)∣>0,z−z0=ρeiθ,
the Taylor expansion f(z)=f(z0)+21f′′(z0)(z−z0)2+O(ρ3) gives
the level lines through z0 are the two perpendicular lines made up of the four directions with cos(2θ+α)=0;
the lines of constant phase are the two perpendicular lines made up of the four directions with sin(2θ+α)=0, bisecting the level lines at π/4;
of the latter, the pair with cos(2θ+α)=−1, that is
θ=2π−αorθ=2π−α+π,
are the directions in which ∣eλf∣ falls off fastest. These two are antiparallel, so together they form a single smooth curve through z0: the line of steepest descents. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.1, p. 908) point out that the plural is deliberate — the modulus decreases on both sides of the saddle.
A valley is a region in which φ→−∞, so that the integrand tends to zero; the level lines through the saddles divide the plane into valleys and the hills between them. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, pp. 908-909) call these regions sectors; earlier, in §25.8.1 (Riley, Hobson & Bence (2006), §25.8.1, p. 907), the same authors observe that a steepest path cannot close on itself, since φ decreases monotonically along it until it meets a singularity of f.
Let f and g be analytic in a region containing the contour C, let z0 be a saddle point of f with f′′(z0)=Aeiα and A=0, and suppose that
C may be deformed, without crossing a singularity of f or g and without changing its endpoints or the valleys in which they lie, into a contour C0 that crosses z0 along the line of steepest descents;
g is continuous and non-zero at z0 and varies slowly compared with eλf near z0;
the contribution to ∫C0 from outside any fixed neighbourhood of z0 is exponentially smaller than the contribution from inside it.
Then, as λ→+∞,
∫Cg(z)eλf(z)dz∼g(z0)eλf(z0)eiθλA2π,
where θ is the direction, in the sense of traversal of C0, in which the deformed contour passes through z0:
θ=2π−αor2π−α+π,
the choice being fixed by the orientation of C. If several saddle points must be crossed, the contributions are added.
The three books state this identically up to bookkeeping. Hassani (Hassani (2013), §12.5, Eq. (12.34), p. 385) writes I(α)≈eαf(z0)2π/αeiθ1g(z0)∣f′′(z0)∣−1/2 with the large parameter α explicit, as here. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Eqs. (12.106), (12.108), pp. 588-589) absorb λ into their exponent w and write g(z0)F(z0)eiθ2π/∣w′′(z0)∣ with θ=−21argw′′+(2π or 23π), which is the same pair of directions. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, Eq. (25.66) p. 909, sign rule p. 910) write ±2π/Ag(z0)ef0ei(π−α)/2 and resolve the sign, on the following page, by the rule that it is + when the direction of traversal satisfies −2π<θ≤2π; since −ei(π−α)/2=ei[(π−α)/2+π], their ± is exactly the choice between the two values of θ above.
The two admissible directions differ by π and hence by an overall minus sign in the answer, so a convention is needed and the books choose differently.
Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, pp. 909-910) keep the factor ei(π−α)/2 and attach ± in front, taking + when the traversal direction θ lies in (−2π,2π] — in their words, when the path crosses the saddle in the direction of increasing real part.
Hassani (Hassani (2013), §12.5, Eq. (12.31), p. 384) fixes 0≤θ1<π and orients the real parameter t accordingly, so that no ± appears in his formula; the price is that the orientation of the contour must be checked against his choice, and reversed with a minus sign if it disagrees.
Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Eq. (12.106), p. 588) list both values of θ and say only that the choice is settled by the sense of traversal at z0, and that nothing but the sign depends on it.
This unit follows the third: θis the traversal direction, and no separate sign appears. Whichever convention is used, the orientation of the original contour is what decides, and getting it wrong changes the sign of the answer — in the worked example below it would flip the sign of a cosine.
Idea of the proofOriginal
Two obstacles stand between ∫geλfdz and Laplace's method. The exponent is complex, so the integrand oscillates and the contributions cancel rather than concentrate; and even its modulus eλφ has no peak, because a harmonic function has no interior maximum. Cauchy's theorem removes both at once. The contour may be moved, and we move it to a curve on which Imf is constant. On such a curve the oscillating factor is a single constant phase that comes outside the integral, and what remains is a real integral with a genuine maximum of eλφ.
Where can that maximum be? Only at a point where the curve of constant phase has a critical point of φ along it, and the no-maximum theorem forces that point to be a saddle: f′(z0)=0. Near it f is quadratic, so on the steepest path the integrand is a Gaussian in the arclength, of width λ−1/2, and Laplace's method finishes the job. The two factors in the answer are then visible in advance: eλf(z0) is the height of the pass, and 2π/(λA) is the width of the Gaussian, with A=∣f′′(z0)∣ measuring how sharply the ground falls away. The one quantity the local picture cannot supply is the direction eiθ in which the path crosses: that depends on where the contour came from and where it is going.
Step 1: move the contour.geλf is analytic in the region between C and C0, so by Cauchy's theorem the two integrals agree, provided the deformation sweeps no singularity and the endpoints are held fixed. When the endpoints are at infinity, it is enough that they stay inside the same valleys throughout, since the integrand is exponentially small there and the joining arcs contribute nothing in the limit. This is hypothesis 1, and it is a global condition: the local analysis below says nothing about whether it holds.
Step 2: freeze the phase. On C0, Imf=ψ(z0) is constant by construction, and φ≤φ(z0) with equality only at z0, so
λ(f(z)−f(z0))=−λs(z),s(z)≥0real.
Put s=τ2 with τ real, negative before the saddle and positive after it. The integral becomes
∫C0geλfdz=eλf(z0)∫−∞∞g(z(τ))dτdze−λτ2dτ,
an integral of exactly the shape treated in Laplace's method, with a maximum of the exponent at the interior point τ=0.
Step 3: the local change of variable. Near z0, f(z)−f(z0)=21Aeiα(z−z0)2+O((z−z0)3). Dropping the cubic remainder and writing z−z0=ρeiθ with θ the steepest-descent direction, so that ei(2θ+α)=−1,
This is Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §25.8.2, Eq. (25.65), p. 909) substitution, written there as dz=±2/Aei(π−α)/2ds; the sign they carry is our choice of which of the two antiparallel directions θ names, and it is fixed by the orientation of C because τ must increase along the contour.
Step 4: evaluate. By hypothesis 2, g may be replaced by g(z0) over the width ρ∼(λA)−1/2 of the Gaussian. Counting powers crudely this costs O(λ−1/2) relative to the leading term, but the bound is not sharp. Expanding g(z(τ))dz/dτ in powers of τ, the term of first order is odd in τ and integrates to zero against e−λτ2, so the first surviving correction is the quadratic one and the true relative error is O(λ−1). This is why Hassani's (Hassani (2013), §12.5, Eq. (12.33), p. 385) expansion runs in the Gaussian moments λ−k−1/2, with nothing in between. By hypothesis 3 the limits may be extended to ±∞, the added tails being exponentially small. With ∫−∞∞e−λτ2dτ=π/λ from Gaussian integrals,
Step 5: several saddles. If the deformation of step 1 can only be carried out through a chain of saddles joined by valleys, step 2 applies on each stretch separately and the results add; the error in neglecting the overlap is exponentially small provided the saddles are well separated. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, p. 589) close their treatment by warning that hypothesis 3 — that the saddle neighbourhood really does dominate — has to be verified case by case, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, p. 909) say the same of the multi-saddle case. □
Remark(Where the hypotheses are used)Original
Analyticity is used only in step 1, and it is used globally. Everything else is Laplace's method on a curve. But step 1 is where the method can silently fail: if the deformation from C to C0 crosses a pole of g or a branch cut of f, the answer acquires a residue or a cut integral that the saddle-point formula knows nothing about. The deformation must also respect the valleys containing the ends of C; a path that would end on a hill is not a legitimate deformation however much better its saddle looks.
f′′(z0)=0 sets the width and the direction. If f′′(z0)=0 as well, two saddles have coalesced; the local exponent is cubic, the width of the pass is λ−1/3 rather than λ−1/2, and the answer involves an Airy function (Airy functions) instead of a Gaussian. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.1, p. 907) note that the Taylor expansion must then be carried to the first non-vanishing term. The same failure appears as an infinity in the formula above, since A→0.
"g slowly varying" is a real restriction, not a formality.Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, p. 909) require that g have no zeros or singularities near the saddle and that rapidly varying factors be moved into the exponent, usually as a logarithm. A factor zν with ν comparable to λmoves the saddle, and treating it as a constant prefactor gives the wrong answer; written as eνlnz it correctly shifts f. A zero of g at z0 makes the leading term vanish and the true leading behaviour is O(λ−3/2), from the next term of the expansion.
The contribution of the rest of the contour is assumed, not proved. Hypothesis 3 is what makes the extension of the integration limits to ±∞ harmless, and it is the hypothesis most often left unchecked. It holds automatically when φ decreases to −∞ along C0 on both sides, since the discarded tails are then smaller than eλφ(z0)−λδ for some δ>0.
With several saddles, the largest Ref(z0) wins — but only among those the contour actually crosses. A saddle with a bigger Ref that is not on any admissible deformation contributes nothing at all; the worked example below has one such, and if it contributed the integral would be exponentially large instead of exponentially small. When the deformation does cross two saddles of equal Ref, their contributions interfere and the answer oscillates in λ. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, p. 912) discuss the other reason to prefer one saddle over another: a nearby second saddle spoils the quadratic approximation even when it is subdominant.
The formula is the first term of a series, and the series is asymptotic.Hassani (Hassani (2013), §12.5, Eq. (12.33), p. 385) continues the expansion by developing g in a Taylor series about z0 and integrating term by term against the Gaussian moments ∫e−λτ2τ2kdτ=λ−k−1/2Γ(k+21). In that calculation the change of variable of step 3 is used in its quadratic form, so his higher coefficients account for the variation of g; the contributions of the cubic and higher terms of f enter the same orders and must be expanded as well, as is done explicitly for the Gamma integral in Stirling's formula. The resulting series is asymptotic, not convergent (Asymptotic series).
The real-axis case is Laplace's method and the level-line case is stationary phase. If C already passes through z0 along the steepest path — as it does for the Gamma integral — no deformation is needed and the method reduces to Laplace's method. If instead the contour is taken along a level line, where the modulus rather than the phase is constant, the same final formula is obtained by a different route: that is The method of stationary phase, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.3, p. 917) verify that their two formulas agree.
Method(Carrying out a saddle-point calculation)Original
Get the integral into the form ∫Cg(z)eλf(z)dz, with λ large and real and with every rapidly varying factor inside the exponent. Powers such as zν with large ν, and factorials, belong in f as logarithms; only genuinely slowly varying factors stay in g. If the large parameter is buried, rescale the variable first to expose it, as Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, p. 911) do with s=x1/2t for the Airy integral.
Find the saddles by solving f′(z)=0, and compute f(z0) and f′′(z0)=Aeiα at each.
Map the landscape. Determine the valleys, that is the sectors at infinity in which Ref→−∞, and locate the endpoints of C in them. Then decide which chain of saddles joins the starting valley to the finishing valley. This is the step that decides the answer, and it is not local: a saddle with a larger Ref that lies on no admissible chain does not contribute.
Fix the direction. At each saddle on the chain, the steepest-descent line runs along θ=(π−α)/2 and θ=(π−α)/2+π; take θ to be whichever of the two points along the direction of traversal, which is determined by which valley the path is coming from.
Write down one term per saddle, g(z0)eλf(z0)eiθ2π/(λA), and add them.
Check. Compare Ref at the saddles used with any crude bound on the integral: a prediction that grows faster than the bound means the wrong chain was chosen. If the original integral is real, the sum must come out real, which for a conjugate pair of saddles means the two terms are conjugate — a sensitive test of the directions chosen in step 4. Finally evaluate the integral numerically at one moderate λ; for analytic f and g at a non-degenerate saddle the correction to each saddle's term is O(λ−1) relative to it, there being no half-integer order between.
Example(A quartic phase with two contributing saddles)Original
Take g≡1 and
I(λ)=∫−∞∞exp[λ(it−41t4)]dt,f(z)=iz−41z4.
On the real axis the integrand has modulus e−λt4/4, so the integral converges absolutely and satisfies the crude bound
∣I(λ)∣≤∫−∞∞e−λt4/4dt=21Γ(41)(4/λ)1/4,
which we will use in step 6. Replacing t by −t and conjugating leaves the integrand unchanged, so I(λ) is real.
Saddles.f′(z)=i−z3=0 gives z3=i, so
z1=eiπ/6,z2=ei5π/6,z3=ei3π/2=−i.
At a saddle z4=z⋅z3=iz, so f(zk)=izk−41izk=43izk, which saves all the arithmetic:
giving A=3 at both, with α=4π/3 at z1 and α=2π/3 at z2.
The landscape. For large ∣z∣ the quartic dominates and Ref≈−41∣z∣4cos4θ, so the valleys are the four sectors with cos4θ>0, centred on θ=0,π/2,π,3π/2; call them V0,V1,V2,V3. The contour runs from −∞, in V2, to +∞, in V0.
The steepest path through z3=−i is the imaginary axis: on z=iτ the exponent f=−τ−41τ4 is real with a maximum at τ=−1, so that path descends from z3 into V3 one way and into V1 the other. It therefore joins V3 to V1 and cannot be part of any path from V2 to V0. The chain that is available is
V2⟶z2⟶V1⟶z1⟶V0,
so z1 and z2 contribute and z3 does not. The bound above confirms this: Ref(z3)=3/4>0, so a contribution from z3 would be of size e3λ/4, impossible for a quantity bounded by O(λ−1/4).
Directions. At each saddle θ=(π−α)/2 or that plus π:
at z2 (reached first, coming from the lower left and leaving to the upper right): θ2=(π−2π/3)/2=π/6;
at z1 (reached second, coming from the upper left and leaving to the lower right): θ1=(π−4π/3)/2=−π/6.
Both point into the right half plane, consistent with a contour traversed from left to right, which is Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §25.8.2, p. 910) rule of thumb.
The two terms are complex conjugates, so the sum is real, as it must be:
I(λ)∼23λ2πe−3λ/8cos(833λ−6π)
The number. At λ=20: e−3λ/8=e−7.5, 22π/60=0.6472086, 833⋅20=12.99038, so the argument of the cosine is 12.46678 and the cosine is 0.9950452. The prediction is
I(20)≈3.561873×10−4.
Numerical quadrature of the original real integral gives
I(20)=3.541078×10−4,
a relative error of 0.59%. At λ=40 the prediction is 1.326920×10−7 against the exact 1.325479×10−7, a relative error of 0.11%. The O(λ−1) accuracy the method claims is accuracy relative to each saddle's own contribution, and not relative to their sum: the error of the sum, measured against the exact value, is +5.9×10−3 at λ=20, +1.1×10−3 at 40, −9.3×10−4 at 60 and −4.3×10−3 at 240, so it does not fall like λ−1 at all. The reason is that the first correction shifts the argument of the cosine as well as the amplitude, and a shift of size λ−1 in the argument costs little near a maximum of the cosine and arbitrarily much, in relative terms, near a zero of it, where the sum itself is small.
What the example shows. The answer is exponentially small, of size e−3λ/8, even though the integrand has modulus 1 at t=0: the smallness is pure cancellation, and no amount of care on the real axis would reveal it. It is also oscillatory in λ, and the oscillation comes from the imaginary parts of f at the two saddles — a single saddle can never produce it. Finally, had the directions been taken as θ1=5π/6 and θ2=7π/6, the result would have been the same cosine with the opposite sign, and at λ=20 the prediction would have been wrong not by 0.59% but by 200%.
ApplicationOriginal
Stirling's formula and the Gamma function. The oldest and still the most used application is Γ(λ+1)=∫0∞e−ttλdt. Writing tλ=eλlnt — step 1 of the method, moving the rapidly varying factor into the exponent — and rescaling t=λu puts the integral in standard form with a saddle at u=1 on the real axis, where the contour already runs along the steepest path, so θ=0 and the formula gives 2πλλλe−λ at once. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Example 12.7.1, p. 589) and Hassani (Hassani (2013), §12.5, Example 12.5.1, p. 385) both present exactly this calculation, and Arfken remarks that the method is probably the easiest route to the leading term. The corrections, and the extension to complex argument, are the subject of Stirling's formula.
Asymptotics of Bessel, Hankel and Airy functions. Every special function with a contour-integral representation yields its large-argument behaviour this way, and the structure of the answer is read off the saddle geometry. Hassani (Hassani (2013), §12.5, Example 12.5.2, p. 386) applies the formula to Hν(1)(λ)=iπ1∫Ce(λ/2)(z−1/z)z−ν−1dz, where f′=0 gives z0=±i and the contour selects z0=+i, with f(i)=i producing the travelling-wave factor eiλ and θ1=3π/4 producing the phase −π/4 that appears in every Bessel asymptotic formula (Asymptotic forms of Bessel functions, Integral representations of Bessel functions). Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.8.2, pp. 911-912) do the same for Ai(x) with x>0: the saddles sit at t=±i after rescaling, the one at +i is chosen, and the result (2πx1/4)−1e−2x3/2/3 is the decaying branch that identifies the function. The oscillatory branch, for x<0, comes from two real stationary points and belongs to The method of stationary phase; the fact that the same function is governed by different saddles in different sectors of the complex plane is the Stokes phenomenon of The WKB approximation.
Statistical mechanics: the partition function and the equivalence of ensembles.Hassani (Hassani (2013), §12.5, p. 382) introduces the method precisely as Debye's technique for the contour integral that represents the partition function. The canonical pattern is that the microcanonical density of states is recovered from the canonical partition function by an inverse Laplace transform,
with N the particle number playing the role of λ. The saddle-point condition is the thermodynamic relation that fixes the temperature given the energy per particle, the exponent at the saddle is N times the entropy, and the Gaussian width gives the O(N−1/2) relative fluctuations. That the two ensembles agree to leading order, and differ by a factor 2πNσ2 in the prefactor, is the saddle-point formula. The same calculation in probabilistic language is the Legendre transform at the heart of Large deviations and Cramer's theorem, and the Gaussian correction is the local central limit refinement of The central limit theorem.
Inverting Laplace transforms. The Bromwich integral f(t)=2πi1∫c−i∞c+i∞fˉ(s)estds is in standard form with λ=t, and the long-time behaviour of f follows from the saddle of st+lnfˉ(s) — or, when fˉ has a branch point, from the steepest path wrapped around the cut. Diffusion fronts, dispersive pulse arrivals and the late-time tails of Green functions are all computed this way; see The Laplace transform and Integrals around branch cuts.
Numerical mathematics: quadrature along the steepest path. The formula's own construction suggests an algorithm. Deforming an oscillatory integral onto the steepest-descent contour and then applying Gauss-Hermite quadrature turns a problem whose cost grows with λ into one whose cost falls with λ, because the integrand becomes a smooth Gaussian bump of width λ−1/2 with no oscillation at all. This is the basis of numerical steepest descent for highly oscillatory integrals, and it is the practical counterpart of Numerical quadrature failing on ∫eiλf for large λ. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Example 12.7.2, p. 590) make the point in miniature: for H(t)=∫e−t(z2−1/4)cos(tz)(1+z2)−1dz at t=100 the integrand on the real axis oscillates so fast that ordinary quadrature is painful, while the saddle at z0=i/2 gives 0.23633 against the exact 0.23284 with no numerical work at all.
Derives Γ(x+1)∼2πxxxe−x from the Euler integral by Laplace's method, carries the expansion of the exponent far enough to produce the first correction 1+1/(12x), records the full Stirling series in Bernoulli numbers together with its range of validity and its divergence, and tests the result on a binomial coefficient.
Definition(Stirling's formula, Stirling's series, and the notation)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.4, Eqs. (13.60) and Ex. 13.4.1, pp. 624-625; Ahlfors (1979), §2.5, Eqs. (37)-(38), p. 204
so that Γ(n+1)=n! for non-negative integers n. The symbol a(x)∼b(x) means a/b→1; it says nothing about a−b, which here grows without bound.
Stirling's formula is the leading-order statement Γ(x+1)∼2πxxxe−x. Stirling's series is the asymptotic expansion of which it is the first term. The series is normally written for the logarithm, where the coefficients are simplest,
with the Bernoulli numbers B2=61, B4=−301, B6=421 of Riemann zeta function and Bernoulli numbers. Ahlfors (Ahlfors (1979), §2.5, Eqs. (40)-(41), p. 205) writes the same coefficients as Cν=(−1)ν−1Bν/[(2ν−1)2ν] in the older numbering, in which B1=61 and B2=301 are positive and carry half the index; both conventions give 12z1 for the first term, and this unit uses the modern one throughout.
The word series is a courtesy: as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.4, p. 623) point out, the Bernoulli numbers grow so fast that the sum diverges for every z. It is an asymptotic series in the sense of Asymptotic series, useful when truncated and useless when summed.
Example(Stirling's formula and its first corrections)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.7 Example 12.7.1 p. 589, §13.4 Eqs. (13.60) and Ex. 13.4.1, pp. 624-625; Ahlfors (1979), §2.5, Eqs. (36)-(38), pp. 203-204; Grimmett & Stirzaker (2020), §5.9, Eqs. (6)-(7), p. 214
1. Leading order. As x→+∞ through real values,
Γ(x+1)∼2πxxxe−x,equivalentlyn!∼2πnnne−n.
Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §5.9, Eqs. (6)-(7), p. 214) state the second form, attributing it in a footnote to de Moivre, and the equivalent Γ(t)∼tt−1/2e−t2π for the Gamma function itself.
2. First correction. More precisely,
Γ(x+1)=2πxxxe−x[1+12x1+O(x21)],
and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.4, Exercise 13.4.1, p. 625) continue the bracket as 1+12z1+288z21−51840z3139+⋯.
valid in the right half plane, with J(z)→0 as z→∞ in any half plane Rez≥c>0, and shows that expanding the integrand in powers of 1/z generates the series of part 3. Since lnΓ(z+1)=lnz+lnΓ(z), parts 3 and 4 are the same statement. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Exercise 12.7.3, p. 591) likewise record that the leading term holds for complex argument provided the real part is large and positive.
Idea of the proofOriginal
Two things have to happen before Laplace's method can be used, and they are the whole derivation. The integrand txe−t has its maximum at t=x, which runs away as x grows, so the first move is to rescale, t=xu, which pins the maximum at u=1 for every x. The second move is to recognise that ux is not a slowly varying prefactor but a rapidly varying one, and to put it where it belongs, in the exponent, as exlnu. After both moves the integral reads xx+1e−x∫0∞e−xφ(u)du with φ(u)=u−1−lnu, a fixed function, independent of x, with a single non-degenerate minimum at u=1 where φ=0 and φ′′=1. Laplace's method then returns 2π/x, and the three factors multiply to Stirling's formula.
The correction comes from the same picture pushed one step. The Gaussian that dominates has width v∼x−1/2 about the minimum, so in the expansion φ=21v2−31v3+41v4−⋯ the cubic term contributes xv3∼x−1/2 and the quartic xv4∼x−1. Expanding the exponential of these small quantities and integrating against the Gaussian, the cubic term dies by symmetry at first order but survives squared, at the same order x−1 as the quartic. Two contributions, −43x−1 from the quartic and +65x−1 from the squared cubic, leave +121x−1. The famous 1/12 is that near-cancellation.
One thing the local analysis cannot produce is 2π itself — it comes from the Gaussian integral, that is from the global normalisation, which is why every proof of Stirling's formula has to pass through some identity that knows about π.
This is the step that The method of steepest descents insists on: a factor varying as fast as exlnu may not be treated as a slowly varying prefactor, because it moves the maximum. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Example 12.7.1, p. 589) and Hassani (Hassani (2013), §12.5, Example 12.5.1, p. 385) make the same two moves, Hassani keeping x inside the exponent as lnu−u/x with the saddle at u=x, which is the same calculation before rescaling.
Step 3: check the hypotheses of Laplace's method. On (0,∞),
φ′(u)=1−u1,φ′′(u)=u21,
so φ has exactly one stationary point, at u=1; there φ(1)=0 and φ′′(1)=1>0, so it is a non-degenerate interior minimum; and φ>0 elsewhere, with φ(u)→+∞ at both ends (−lnu at 0, u at ∞). These are exactly the hypotheses of Laplace's method, which therefore gives
∫0∞e−xφ(u)du∼xφ′′(1)2π=x2π.
Step 4: leading order. Combining steps 2 and 3,
Γ(x+1)∼xx+1e−xx2π=2πxxxe−x,
which is part 1. This is the computation Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §5.9, Exercise 5.9.6, p. 216) set as an exercise, with the maximum of fn(u)=nlnu−u at u=n before rescaling.
Step 5: the first correction. Put u=1+v, so that φ(1+v)=v−ln(1+v) and, for ∣v∣<1,
φ(1+v)=2v2−3v3+4v4−⋯.
Hence
∫0∞e−xφdu=∫−1∞e−xv2/2exp[x(3v3−4v4+⋯)]dv.
By the localisation that underlies Laplace's method, only ∣v∣≲x−1/2+δ matters: outside that range the Gaussian factor is smaller than any power of 1/x, so the lower limit may be moved from −1 to −∞ at an exponentially small cost, and inside it the bracket is small and its exponential may be expanded. Counting powers with v∼x−1/2: the term xv3 is of order x−1/2, the terms xv4 and (xv3)2 are of order x−1, and everything else is smaller. Keeping these,
exp[x(3v3−4v4+⋯)]=1+3xv3−4xv4+18x2v6+⋯.
Step 6: integrate term by term. The Gaussian moments are ∫−∞∞e−xv2/2v2kdv=2π/x(2k−1)!!x−k and the odd moments vanish, so ⟨v3⟩=0, ⟨v4⟩=3x−2 and ⟨v6⟩=15x−3, each relative to 2π/x. Therefore
and −43+65=−129+1210=121, so the bracket is 1+12x1+O(x−2). Multiplying by xx+1e−x gives part 2. The vanishing of the odd moments also explains why the next correction is O(x−2) and not O(x−3/2): every term of half-integer order carries an odd power of v.
Step 7: the series and the constant. Repeating steps 5 and 6 to higher order reproduces 288z21 and the later coefficients, but the labour grows quickly, and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.4, p. 622) recommend obtaining the full series instead from the Euler-Maclaurin formula of Euler-Maclaurin summation, applied to ∫0∞(z+t)−2dt and integrated twice; that route produces (13.57) and (13.58) with two unknown constants of integration, which are then fixed at C1=0 and C2=21ln2π by requiring consistency with the Legendre duplication formula. Ahlfors (Ahlfors (1979), §2.5, pp. 203-204) fixes the same two constants differently, the first from Γ(z+1)=zΓ(z) and the second from the reflection formula Γ(z)Γ(1−z)=π/sinπz at z=21+iy. □
Remark(Where the hypotheses are used, and what the derivation does not prove)Original
Rescaling is not cosmetic. Laplace's method requires the maximum to sit at a fixed interior point while the large parameter grows; in the unrescaled integral the maximum is at t=x and the width of the peak is x, so there is no small parameter at all. Step 1 is what creates one.
φ′′(1)=1=0 is what makes the answer 2π/x rather than something else. A degenerate minimum would give a different power of x and a different constant, as in the Airy-type cases of The method of stationary phase. Here non-degeneracy is automatic, since φ′′(u)=u−2 never vanishes.
The uniqueness of the minimum and the growth of φ at both ends are what let the tails be discarded.φ(u)→∞ as u→0+ only logarithmically, so the contribution from near t=0 is O(x−x) rather than merely exponentially small; it is still negligible, but the estimate is worth making rather than assuming.
Step 5 is formal where it rearranges an asymptotic expansion, and the justification is borrowed. Expanding ex(v3/3−⋯) and integrating term by term over all of R is not a convergent manipulation; it is legitimate as an asymptotic statement because of the localisation argument quoted at the start of the step, which is the same one that proves Laplace's method and is proved in Laplace's method and Watson's lemma. Carrying it to all orders produces a divergent series, so no amount of care turns the derivation into an exact identity; for that one needs the remainder formula of Ahlfors (Ahlfors (1979), §2.5, Eq. (36), p. 203).
The constant 2π is global, not local. Here it enters through the Gaussian integral. Every other derivation must import it from somewhere equally non-local: Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.4, p. 624) take it from the Legendre duplication formula, Ahlfors (Ahlfors (1979), §2.5, p. 204) from the reflection formula, and Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §5.9, p. 215) and Folland (Folland (1999), §10.3, Exercise 21, p. 327) from the normalisation of the standard normal distribution via the central limit theorem. An elementary argument comparing lnn! with ∫ln gives n!∼Cnn+1/2e−n easily and stops there, unable to evaluate C.
Real x→+∞ is what the derivation above covers; the complex statement needs the contour. Nothing in steps 1-6 is complex, so they establish the formula for large positive real argument only. Its extension to Rez≥c>0 is Ahlfors's (Ahlfors (1979), §2.5, pp. 203-204) theorem, with the bound on J(z) obtained by splitting the defining integral at η=∣z∣/2; the saddle-point derivation extends the same way, as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.7, Exercise 12.7.3, p. 591) indicate. The formula fails near the negative real axis for the obvious reason that Γ has poles there.
The series diverges, so "more terms" is not always "more accuracy". The term B2n/[2n(2n−1)z2n−1] eventually grows for fixed z, because ∣B2n∣ grows faster than any power. The best accuracy attainable at a given z is reached at the smallest term, and is roughly its size; see the optimal-truncation discussion in Asymptotic series. At z=10 the smallest term is far below machine precision, which is why Stirling's series is the standard way to compute lnΓ for large arguments; at z=1 it is useless after two terms.
Method(Using Stirling's formula in practice)Original
Work with logarithms.n! overflows for n beyond a few hundred while lnn! does not, and almost every use of the formula is inside a ratio or an exponent anyway. Use
lnΓ(x+1)≈21ln(2πx)+xlnx−x+12x1−360x31.
Shift the argument up before applying it. The error falls like a power of 1/x, so for small x apply the recurrence Γ(x+1)=xΓ(x) a few times to move the argument to x+m with m around 10, evaluate the series there, and divide back out. This is how library implementations of lnΓ are organised.
Cancel before approximating. In a ratio such as (n2n) or n!/[(n−k)!k!], assemble the logarithms first and let the xlnx terms cancel algebraically; approximating each factorial separately and then dividing loses most of the accuracy to subtractive cancellation between large numbers.
Decide the number of terms from the target accuracy. The relative error of the m-term truncation is about the size of the first omitted term. One term gives ≈1/(12x), two give ≈1/(288x2): for three significant figures, one term suffices from x≈100, two from x≈6.
Keep track of which form is being used.Γ(x+1), Γ(x) and x! differ by factors of x, and the exponents x+21 and x−21 are easy to interchange by accident. The safest check is the leading behaviour of the logarithm, lnΓ(x+1)−lnΓ(x)=lnx.
Verify once numerically at a moderate argument, where the exact value is still computable, before trusting the formula at the argument of interest.
Example(Accuracy at small arguments, and a central binomial coefficient)Original
Part A: how good is it already at x=5 and x=10? Writing Sm(x) for the m-term truncation of the product form of part 2, direct computation gives
x
Γ(x+1)
S1/Γ
S2/Γ
S3/Γ
5
120
0.983493
0.999885
1.000021
10
3628800
0.991704
0.999968
1.000003
One term is within 1.7% at x=5 and 0.83% at x=10; two terms bring this to 1.2×10−4 and 3.2×10−5; three terms to 2.1×10−5 and 2.7×10−6. The one- and two-term truncations fall short of Γ and the three-term truncation overshoots it, so the error changes sign between the second and the third term — a reminder that the truncations of an asymptotic series do not bracket the answer in any regular way. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.4, Table 13.1, p. 625) tabulate the first two of these columns for s=1,…,10; our values agree with theirs everywhere except at s=10, where they print 0.99998 for the two-term ratio and we compute 0.9999682, which rounds to 0.99997.
Part B: the number of ways to split 100 objects into two equal halves. The exact value is
a number large enough that computing the three factorials separately is not an option in floating point. Following step 3 of the method, approximate the ratio rather than its parts. With n=50 and the two-term form Γ(m+1)=2πmmme−m(1+12m1+288m21+⋯),
The algebraic part collapses to 4n/πn. For the correction, the denominator squares to 1+6n1+72n21 to this order, whose reciprocal is 1−6n1+72n21; multiplying by the numerator,
since 241−61=−81 and 11521+721−1441=11521+16−8=11529=1281. Hence
(n2n)=πn4n[1−8n1+128n21+O(n−3)].
At n=50: 450=1.2676506×1030 and 50π=12.533141, so the leading term is 1.011438842×1029, high by 0.250%. The bracket is 1−4001=0.9975, giving
(50100)≈1.008910245×1029
against the exact 1.008913445×1029: a relative error of 3.2×10−6, or about eight correct significant figures from a formula with two terms. Including the third term, +128⋅25001, gives 1.008913406×1029, with relative error 3.9×10−8.
What the numbers show. The successive relative errors 2.5×10−3, 3.2×10−6, 3.9×10−8 fall by factors of about 790 and 81: each new term is smaller than the last by a factor of order n, which is what an asymptotic series does while it is still far from its optimal truncation. At n=1 the same three approximations give 2.257, 1.975 and 1.992 for (12)=2, so even there the first correction is worth having, but the series has stopped being a reliable source of further digits.
ApplicationOriginal
Statistical mechanics: entropy and the counting of microstates. Boltzmann's S=kBlnW requires the logarithm of a multinomial coefficient, and Stirling's formula in the crudest form lnN!≈NlnN−N is what turns it into a usable function of the occupation numbers. For W=N!/∏ini! with ni=Npi,
kBS=lnW≈−Ni∑pilnpi,
the Gibbs entropy: the NlnN terms cancel exactly, which is why the crude form suffices and why the answer is extensive. The same cancellation makes the 21ln(2πN) corrections subextensive, of relative order lnN/N, and therefore invisible in the thermodynamic limit — but they are precisely what is kept in the fluctuation theory. The N! of indistinguishability in the ideal-gas partition function enters the Sackur-Tetrode entropy the same way, and the resolution of the Gibbs paradox depends on treating it with Stirling rather than dropping it.
Random walks and the recurrence of diffusion.Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §3.10, Eq. (20), p. 86) apply the formula to the symmetric simple random walk: the probability of being back at the origin after 2k steps is u2k=(k2k)2−2k, and the central binomial asymptotics computed in the example above give
u2k∼πk1.
Since ∑kk−1/2 diverges, the walk returns to the origin infinitely often: the recurrence of the one-dimensional walk is a corollary of Stirling's formula. In two dimensions the same computation gives u2k∼1/(πk), still divergent, and in three u2k∼ck−3/2, convergent — which is Pólya's theorem, and the reason a diffusing particle in space escapes. See Random walks.
The de Moivre-Laplace theorem and the central limit theorem. Historically the formula was invented for this: to find the limiting shape of (kn)2−n. Writing k=n/2+yn/2 and applying the two-term formula to all three factorials gives the Gaussian 2/(πn)e−y2/2 directly, constant included. The logic also runs backwards, and both Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §5.9, pp. 214-215 and Problem 5.12.26, p. 403) and Folland (Folland (1999), §10.3, Exercise 21, p. 327) derive Stirling's formula from the central limit theorem — Grimmett and Stirzaker by Fourier inversion of the characteristic function of a standardised Gamma variable at the single point x=0, and again from a Poisson sum, and Folland from the expectation of the positive part of a standardised Poisson sum. The constant 2π is in all cases the normalisation of the normal density. See The central limit theorem.
Information theory and large deviations. With k=pn, the same two-term computation gives
so the number of binary strings of length n with a prescribed fraction of ones is 2nH(p) up to a factor O(n−1/2). This is the counting estimate behind the source-coding theorem, behind the typical-set argument, and behind the rate function of Large deviations and Cramer's theorem, where the exponent nH is the leading term and the Stirling prefactor supplies the polynomial correction that distinguishes a local from a global limit theorem.
Normalising approximate identities and special-function ratios. Ratios of factorials appear wherever a family of functions has to be normalised. Stakgold and Holst (Stakgold & Holst (2011), §2.2, Exercise 2.2.7(a), p. 122) use k!∼kke−k2πk to show that Ik=22k+1(k!)2/(2k+1)! satisfies kIk→π, which is exactly the statement that the Landau kernels (1−x2)k/Ik form a delta family — the first step of their proof of the Weierstrass approximation theorem. The same ratio, with Γ(k+21)/Γ(k+1)∼k−1/2, controls the normalisation of Legendre and Gegenbauer polynomials and the large-order behaviour of Bessel coefficients in Asymptotic forms of Bessel functions.