Laplace's Method for Integrals

Reach for almost any integral in physics that carries a large parameter — an action in units of $\hbar$ , a temperature written as an inverse $\beta$ , a particle number $N$ — and it tends to arrive in the shape $I(\lambda)=\int g(x)\,e^{\lambda f(x)}\,\mathrm{d}x$ with $f$ real and $\lambda$ large. Such integrals rarely have closed forms, yet are almost always dominated by one point: as $\lambda$ climbs, the weight $e^{\lambda f(x)}$ towers wherever $f$ is largest and collapses to nothing everywhere else, so the whole integral is decided in a shrinking neighbourhood of the maximum. Laplace's method is the disciplined way to cash this in — locate the dominant maximum $x_0$ , expand the exponent to second order, recognise the resulting Gaussian bump, and integrate it exactly. The payoff is a formula of striking economy — area equals peak value times peak height times an effective width $\sim\lambda^{-1/2}$ — plus a systematic ladder of corrections in powers of $1/\lambda$ built from the higher derivatives of $f$ and $g$ . That ladder does not converge; it is asymptotic in Poincaré's sense, superb when truncated near its optimal order and divergent if pushed past it. The method is also the gateway to a whole family: rotate the exponent to imaginary values and it becomes the method of stationary phase; lift the contour into the complex plane and deform it across a saddle and it becomes the method of steepest descent, of which the real-axis Laplace case is the parent. Underneath statistical mechanics, the central limit theorem, semiclassical quantum mechanics, and the asymptotics of every special function, the same Gaussian-peak idea is quietly doing the work.

Key takeaways

Laplace's method gives the leading large- $\lambda$ asymptotics of $I(\lambda)=\int_a^b g(x)\,e^{\lambda f(x)}\,\mathrm{d}x$ when the real exponent $f$ has a sharp interior maximum at $x_0$ .
Near the maximum the exponent is quadratic, $f(x)\approx f(x_0)-\tfrac12|f''(x_0)|(x-x_0)^2$ , so the integrand is a Gaussian of width $\sim\lambda^{-1/2}$ .
The leading result is $I(\lambda)\sim g(x_0)\,e^{\lambda f(x_0)}\sqrt{2\pi/(\lambda|f''(x_0)|)}$ — peak value, peak height, and effective width multiplied together.
Endpoint maxima give $\lambda^{-1}$ scaling (Watson's lemma), a vanishing $f''$ gives $\lambda^{-1/(2p)}$ ; the correction series in $1/\lambda$ is asymptotic and generally divergent.
Worked out for the gamma function it delivers Stirling's formula, $\Gamma(z+1)\sim\sqrt{2\pi z}\,z^z e^{-z}(1+\tfrac{1}{12z}+\cdots)$ , accurate to 0.003% at $z=10$ .
Key idea: for large $\lambda$ the weight $e^{\lambda f(x)}$ is a Gaussian peak of width $\sim 1/\sqrt{\lambda|f''(x_0)|}$ at the maximum $x_0$ ; replace the integrand by that Gaussian and integrate exactly.

The core insight

For large $\lambda$ the weight $e^{\lambda f(x)}$ is a Gaussian peak of width $\sim 1/\sqrt{\lambda|f''(x_0)|}$ centred at the maximum $x_0$ ; replace the integrand by that Gaussian and integrate exactly.

Where this lives today

Laplace's method is not a museum piece; it is the everyday reflex for extracting numbers from stubborn integrals, and it sits at the base of several live research frontiers. In quantum field theory, path integrals with a large action are evaluated by saddle-point expansion about classical instanton solutions, and the Gaussian fluctuation determinant around the saddle — a descendant of the $\sqrt{2\pi/|f''|}$ factor — sets tunnelling rates and vacuum decay probabilities. In statistics and machine learning the "Laplace approximation" is a workhorse for Bayesian inference: a posterior is approximated by a Gaussian at its mode, turning the intractable evidence integral into a Hessian determinant — the engine behind Gaussian-process and neural-network uncertainty estimates. Large-deviation theory, quantifying the exponentially small probabilities of rare events from finance to statistical physics, is at heart a saddle-point analysis of a generating function whose rate function is a Legendre transform. And the asymptotics of the special functions — Bessel, Airy, parabolic-cylinder — are still delivered to arbitrary order by equation (5) and its steepest-descent generalisation. An open direction is resurgence and Borel–Écalle summation, which reconnects the divergent tail of the series to the exponentially small non-perturbative contributions of subdominant saddles.

The physical picture

Watch $e^{\lambda f(x)}$ grow with $\lambda$ . The maximum rises exponentially faster than any neighbouring value, so a broad hump at $\lambda=1$ becomes a needle at $\lambda=100$ . Its height is $e^{\lambda f(x_0)}$ , and how fast it falls away is set by the curvature $f''(x_0)$ : a sharply curved peak is narrow, a flat one wide. The area under a spike is height times width times a shape factor, and for a Gaussian that factor is $\sqrt{2\pi}$ — the entire leading formula in disguise. The first figure makes the collapse visible: three height-normalised weights whose half-widths shrink exactly as $1/\sqrt{\lambda}$ .

Figure

As λ grows the weight e^(λf) collapses to a needle at the maximum x₀. The three curves are height-normalised so only their widths differ, and the half-width shrinks as 1/√λ — the λ = 16 peak is 4× narrower than the λ = 1 peak, so only an O(1/√λ) neighbourhood of x₀ survives

Why replace the true bump by a Gaussian at all? Because the region that matters is only $\sim\lambda^{-1/2}$ wide, and over that shrinking window the higher Taylor terms of $f$ are self-defeatingly small: a cubic contribution scales as $\lambda\cdot(\lambda^{-1/2})^3=\lambda^{-1/2}\to0$ , a quartic as $\lambda^{-1}$ , and only the quadratic term stays $O(1)$ . That is why the Gaussian is not an approximation of convenience but the honest leading behaviour, and why the corrections organise into powers of $1/\lambda$ . The second figure shows the machine on the gamma-function integral, where the payoff is a real number: $10!$ recovered to a fraction of a percent from a single Gaussian.

Figure

Laplace's method applied to the gamma integral for z = 10. The reduced integrand e^(z(ln s − s + 1)) is a sharp peak at s = 1 (solid blue); its quadratic Gaussian approximation e^(−z(s−1)²/2) is almost indistinguishable (dashed amber). Integrating the Gaussian gives Γ(11) ≈ 3.599×10⁶, and the +1/(12z) correction gives 3.6287×10⁶, versus the exact 10! = 3 628 800 — an error of 0.003%

The core equations

Everything begins with a single archetype: a real exponent multiplied by a large parameter, weighted by a slowly varying amplitude. The behaviour as $\lambda\to+\infty$ is decided by where $f$ is largest, not by the fine detail of $g$ .

I(\lambda)=\int_a^b g(x)\,e^{\lambda f(x)}\,\mathrm{d}x,\qquad \lambda\to+\infty \qquad (1)

where: $\lambda$ — large positive parameter setting the sharpness of the peak; $f(x)$ — real exponent, the "phase", whose maximum dominates; $g(x)$ — slowly varying amplitude or prefactor; $x_0$ — location of the global maximum of $f$ on $[a,b]$ . For an interior maximum with $f'(x_0)=0$ and $f''(x_0)<0$ , the quadratic expansion turns the integrand into a Gaussian, and the answer is the product of three transparent factors: peak height, amplitude, and width.

I(\lambda)\sim g(x_0)\,e^{\lambda f(x_0)}\sqrt{\dfrac{2\pi}{\lambda\,|f''(x_0)|}} \qquad (2)

where: $e^{\lambda f(x_0)}$ — peak height; $g(x_0)$ — amplitude at the peak; $\sqrt{2\pi/(\lambda|f''(x_0)|)}$ — effective width times $\sqrt{2\pi}$ ; $|f''(x_0)|=-f''(x_0)$ since the curvature is negative. A sharper peak (larger $|f''|$ ) is narrower and contributes less. The engine behind both the leading term and all corrections is the Taylor expansion of the exponent: the linear term is absent because $x_0$ is stationary, the quadratic term builds the Gaussian, and the cubic and quartic terms feed the $1/\lambda$ series.

f(x)=f(x_0)+\tfrac12 f''(x_0)(x-x_0)^2+\tfrac16 f'''(x_0)(x-x_0)^3+\tfrac1{24}f^{(4)}(x_0)(x-x_0)^4+\cdots \qquad (3)

where: $f''(x_0)<0$ — curvature that fixes the width; $f'''(x_0),f^{(4)}(x_0)$ — higher derivatives that generate the corrections, each weighted by an extra half-power of the small width $\lambda^{-1/2}$ . To turn those Taylor terms into numbers we need the moments of a Gaussian against $t^{2n}$ . Only even powers survive — odd powers integrate to zero by symmetry, which is why the first correction is $O(1/\lambda)$ and not $O(\lambda^{-1/2})$ .

\int_{-\infty}^{\infty} t^{2n}\,e^{-\frac{a}{2}t^{2}}\,\mathrm{d}t=\dfrac{(2n-1)!!}{a^{n}}\sqrt{\dfrac{2\pi}{a}} \qquad (4)

where: $a=\lambda|f''(x_0)|$ — the inverse variance; $(2n-1)!!=1\cdot3\cdot5\cdots(2n-1)$ — the double factorial; odd moments vanish. The normalised averages are $\langle t^2\rangle=1/a$ , $\langle t^4\rangle=3/a^2$ , $\langle t^6\rangle=15/a^3$ . Assembling those moments term by term promotes the leading formula to a full asymptotic expansion — a power series in $1/\lambda$ multiplying the Gaussian.

I(\lambda)\sim g(x_0)\,e^{\lambda f(x_0)}\sqrt{\dfrac{2\pi}{\lambda|f''(x_0)|}}\left[1+\dfrac{c_1}{\lambda}+\dfrac{c_2}{\lambda^2}+\cdots\right] \qquad (5)

where: $c_1,c_2,\dots$ — pure numbers built from the derivatives of $f$ and $g$ at $x_0$ . The series is asymptotic, not convergent: it is meant to be truncated, and its accuracy at fixed truncation improves as $\lambda$ grows. The first coefficient is worth memorising, since it fixes the leading fractional error of the bare Laplace estimate. Working carefully with the signs — recall $f''(x_0)<0$ — the moments of equation (4) combine into a compact expression.

c_1=\dfrac{f^{(4)}}{8\,(f'')^{2}}-\dfrac{5\,(f''')^{2}}{24\,(f'')^{3}}-\dfrac{g''}{2\,g\,f''}+\dfrac{g'\,f'''}{2\,g\,(f'')^{2}}\Bigg|_{x_0} \qquad (6)

where: all quantities are evaluated at $x_0$ , with $f''<0$ ; for constant amplitude $g\equiv1$ only the first two terms survive. Because $(f'')^3<0$ the middle term is positive, delivering the gamma integral's $+1/(12z)$ — a sign that is essential and easy to get backwards. The single most important instance is the gamma function, whose integral representation is tailor-made for the method and yields Stirling's formula in three lines.

\Gamma(z+1)\sim \sqrt{2\pi z}\;z^{z}\,e^{-z}\left(1+\dfrac{1}{12z}+\cdots\right)\quad (z\to\infty) \qquad (7)

where: $\Gamma(z+1)=\int_0^\infty t^{z}e^{-t}\,\mathrm{d}t$ equals $n!$ at integer $z=n$ . The bracket is equation (5) for $f(s)=\ln s-s$ , and the $1/(12z)$ is exactly $c_1/\lambda$ with $\lambda=z$ .

Derivation ALeading-order Laplace formula from scratch

The goal is to derive the $\lambda^{-1/2}$ result rigorously for an interior nondegenerate maximum while tracking honestly what is thrown away, so that "replace the integrand by a Gaussian" becomes a controlled limit rather than a hopeful gesture. Assume $f\in C^2$ has a unique global maximum at an interior point $x_0$ with $f''(x_0)<0$ , and $g$ is continuous with $g(x_0)\neq0$ .

Step 1. Factor out the peak value, $I(\lambda)=e^{\lambda f(x_0)}\int_a^b g(x)\,e^{\lambda[f(x)-f(x_0)]}\,\mathrm{d}x$ . Outside any fixed neighbourhood $|x-x_0|>\delta$ we have $f(x)-f(x_0)\le-\eta<0$ , so that region contributes at most $O(e^{-\lambda\eta})$ . Why: this quarantines the answer into an arbitrarily small window around $x_0$ , and the discarded tails are exponentially smaller than any power of $1/\lambda$ .

Step 2. Rescale to the natural width with $x=x_0+t/\sqrt{\lambda}$ , giving $\lambda[f(x)-f(x_0)]=\tfrac12 f''(x_0)t^2+O(\lambda^{-1/2}t^3)$ and $\mathrm{d}x=\mathrm{d}t/\sqrt{\lambda}$ . Why: this is the only scaling that keeps the quadratic term finite while sending every higher term to zero — what makes $\lambda^{-1/2}$ "the" width.

Step 3. Take the limit: the cubic remainder vanishes pointwise and $g(x_0+t/\sqrt\lambda)\to g(x_0)$ , so dominated convergence gives $I(\lambda)\sim \lambda^{-1/2}e^{\lambda f(x_0)}g(x_0)\int_{-\infty}^{\infty} e^{\frac12 f''(x_0)t^2}\,\mathrm{d}t$ . Why: with the tails already negligible, extending the limits to $\pm\infty$ is free, and freezing $g$ at its peak value is legitimate because $g$ is slowly varying over the shrinking window.

Step 4. Evaluate the Gaussian, $\int_{-\infty}^{\infty} e^{-\frac12|f''(x_0)|t^2}\,\mathrm{d}t=\sqrt{2\pi/|f''(x_0)|}$ , and reassemble to reach $I(\lambda)\sim g(x_0)\,e^{\lambda f(x_0)}\sqrt{2\pi/(\lambda|f''(x_0)|)}$ . Why: the $\lambda^{-1/2}$ from the rescaling times the $\sqrt{2\pi/|f''|}$ from the integral is the whole formula — three physical factors, one line. *For $f(s)=\ln s-s$ with $\lambda=z$ , $x_0=1$ , $|f''|=1$ , the formula reproduces $\Gamma(z+1)\sim\sqrt{2\pi z}\,z^z e^{-z}$ — Stirling's formula.*

Go deeper · Derivation B · 1 minEndpoint maximum and the crossover to Watson's lemma

Derivation BEndpoint maximum and the crossover to Watson's lemma

The goal is to show what changes when the maximum sits at a boundary with nonzero slope rather than at an interior stationary point: the Gaussian gives way to a one-sided exponential, and the scaling drops from $\lambda^{-1/2}$ to $\lambda^{-1}$ — precisely the regime Watson's lemma governs. Let the maximum lie at the left endpoint $x=a$ with $f'(a)<0$ so that $f$ decreases across $[a,b]$ , and $g(a)\neq0$ .

Step 1. Because $a$ is not stationary, the linear term now leads: $f(x)\approx f(a)+f'(a)(x-a)$ with $f'(a)<0$ . Why: at an ordinary endpoint the derivative does not vanish, so the exponent falls off linearly, and that slope — not a curvature — sets the decay.

Step 2. Substitute $u=\lambda|f'(a)|(x-a)$ , so the natural width is $\sim1/\lambda$ rather than $1/\sqrt\lambda$ . Why: a linear exponent decays over a window $\sqrt\lambda$ times narrower than a quadratic one, which is the whole origin of the $\lambda^{-1}$ law.

Step 3. Integrate the one-sided exponential, $I(\lambda)\sim g(a)e^{\lambda f(a)}\int_0^{\infty}\lambda^{-1}|f'(a)|^{-1}e^{-u}\,\mathrm{d}u=g(a)e^{\lambda f(a)}/(\lambda|f'(a)|)$ . Why: only half a peak lies inside the domain, and $\int_0^\infty e^{-u}\,\mathrm{d}u=1$ , leaving the leading term Watson's lemma reproduces. *This result is identical to Watson's lemma applied to $\int_0^\infty g(t)\,e^{-\lambda t}\,\mathrm{d}t\sim g(0)/\lambda$ after the change of variable $t=f(a)-f(x)$ .*

Go deeper · Derivation C · 1 minFirst correction in 1/λ for a simple model

Derivation CFirst correction in 1/λ for a simple model

The goal is to compute the subleading term explicitly, exposing how the higher Taylor coefficients of $f$ feed the asymptotic series and pinning down the sign that is so easy to mishandle. Take $g\equiv1$ and shift the interior maximum to $x_0=0$ with $f(0)=0$ and $f'(0)=0$ ; write $a=|f''(0)|>0$ , $b=f'''(0)$ , $c=f^{(4)}(0)$ , so that $f(x)=-\tfrac{a}{2}x^2+\tfrac{b}{6}x^3+\tfrac{c}{24}x^4+\cdots$ .

Step 1. Scale $x=t/\sqrt\lambda$ so that $\lambda f(x)=-\tfrac{a}{2}t^2+\tfrac{b}{6\sqrt\lambda}t^3+\tfrac{c}{24\lambda}t^4+\cdots$ . Why: the substitution sorts the Taylor terms by powers of $\lambda^{-1/2}$ — the cubic carries $\lambda^{-1/2}$ , the quartic $\lambda^{-1}$ — so every contribution's order is read off at a glance.

Step 2. Expand the non-Gaussian part of the exponential, $\exp(\tfrac{b}{6\sqrt\lambda}t^3+\tfrac{c}{24\lambda}t^4)=1+\tfrac{c}{24\lambda}t^4+\tfrac{b^2}{72\lambda}t^6+\cdots$ . Why: the bare cubic $t^3$ is odd and integrates to zero, so the survivors are the quartic $t^4$ and the cubic-squared $t^6$ , both at $O(1/\lambda)$ — this is why $f'''$ enters only through $(f''')^2$ .

Step 3. Take Gaussian averages with $\langle t^4\rangle=3/a^2$ and $\langle t^6\rangle=15/a^3$ , giving the $1/\lambda$ bracket $\tfrac{c}{24}\tfrac{3}{a^2}+\tfrac{b^2}{72}\tfrac{15}{a^3}=\tfrac{c}{8a^2}+\tfrac{5b^2}{24a^3}$ . Why: both surviving terms carry a plus sign — the exact point where careless algebra flips a sign.

Step 4. Assemble the boxed result, $I(\lambda)\sim\sqrt{2\pi/(\lambda a)}\left[1+\tfrac{1}{\lambda}\left(\tfrac{c}{8a^2}+\tfrac{5b^2}{24a^3}\right)+O(\lambda^{-2})\right]$ . Why: this is equation (5) with $c_1=c/(8a^2)+5b^2/(24a^3)$ , the $g\equiv1$ case of equation (6). *For $f=\ln s-s$ about $s=1$ one has $a=1$ , $b=f'''(1)=2$ , $c=f^{(4)}(1)=-6$ , so $c_1=-6/8+5\cdot4/24=-0.75+0.833=+1/12$ — exactly Stirling's $\Gamma(z+1)\sim\sqrt{2\pi z}\,z^z e^{-z}(1+1/12z+\cdots)$ .*

Limiting cases

The power of $\lambda$ in the prefactor is dictated entirely by the local geometry of the maximum — whether it is interior or at an edge, whether the curvature is finite or vanishing, whether the amplitude survives at the peak.

Case	Result	Meaning
Interior nondegenerate maximum	$I\sim g(x_0)\,e^{\lambda f(x_0)}\sqrt{\dfrac{2\pi}{\lambda\,\vert f''(x_0)\vert }}$	$f'(x_0)=0$ , $f''(x_0)<0$ , $g(x_0)\neq0$ . The generic case; the prefactor scales as $\lambda^{-1/2}$ .
Endpoint maximum (simple slope)	$I\sim g(a)\,e^{\lambda f(a)}\,\dfrac{1}{\lambda\,\vert f'(a)\vert }$	Maximum at $x=a$ with $f'(a)<0$ on $[a,b]$ . Half a peak and a linear exponent give $\lambda^{-1}$ scaling — the Watson's-lemma regime.
Degenerate maximum, $f''=0$	$I\sim g(x_0)\,e^{\lambda f(x_0)}\,\dfrac{1}{p}\,\Gamma\!\left(\tfrac{1}{2p}\right)\left(\dfrac{(2p)!}{\lambda\,\vert f^{(2p)}(x_0)\vert }\right)^{1/2p}$	Lowest nonvanishing derivative is order $2p\ge4$ and negative. The peak is flatter, so it is wider: scaling becomes $\lambda^{-1/(2p)}$ .
Vanishing prefactor, $g(x_0)=0$	$I\sim \dfrac{g''(x_0)}{2}\,e^{\lambda f(x_0)}\dfrac{1}{\lambda\vert f''(x_0)\vert }\left(\dfrac{2\pi}{\lambda\vert f''(x_0)\vert }\right)^{1/2}$	If $g(x_0)=g'(x_0)=0$ the first surviving Gaussian moment is the second, giving an extra $1/\lambda$ suppression.

Common trapThe maximum you expand about must be interior, nondegenerate, and global — check all three

Three assumptions hide inside the tidy $\lambda^{-1/2}$ formula, each with its own failure mode. First, the maximum must be interior: at an endpoint with nonzero slope you are in the Watson's-lemma regime and the scaling is $\lambda^{-1}$ , so equation (2) overcounts by a factor $\sqrt\lambda$ . Second, it must be nondegenerate: if $f''(x_0)=0$ the quadratic approximation is empty and you must keep the lowest nonvanishing even derivative, softening the power to $\lambda^{-1/(2p)}$ . Third, and most treacherous, it must be the global maximum: a taller peak anywhere else dominates completely, so a purely local search can return an answer that is exponentially wrong (equally high maxima both contribute and add).

Worked examples

1 · Stirling's formula, with the factorial 10! recovered to numbers (fully worked)

Apply Laplace's method to $\Gamma(z+1)=\int_0^\infty t^z e^{-t}\,\mathrm{d}t$ , obtain Stirling's formula, and test it on $z=10$ against the exact $10!=3\,628\,800$ .

Step 1. Put the integrand into Laplace form. Scale $t=zs$ so that $t^z e^{-t}=z^z e^{z\ln s}e^{-zs}$ and $\mathrm{d}t=z\,\mathrm{d}s$ , giving $\Gamma(z+1)=z^{z+1}\int_0^\infty e^{z(\ln s-s)}\,\mathrm{d}s$ . The large parameter is $\lambda=z$ and the exponent is $f(s)=\ln s-s$ .

Step 2. Locate and characterise the maximum. $f'(s)=1/s-1=0$ gives the interior maximum $s_0=1$ , with $f(1)=-1$ and $f''(s)=-1/s^2$ , so $f''(1)=-1$ and $|f''(1)|=1$ . The maximum is interior and nondegenerate — equation (2) applies.

Step 3. Assemble the leading term. With $g\equiv1$ , $I\sim e^{z f(1)}\sqrt{2\pi/(z|f''|)}=e^{-z}\sqrt{2\pi/z}$ , so $\Gamma(z+1)\sim z^{z+1}e^{-z}\sqrt{2\pi/z}=\mathbf{\sqrt{2\pi z}\,z^z e^{-z}}$ . This is Stirling's formula.

Step 4. Add the first correction. From equation (6) with $g\equiv1$ , $f'''(1)=2$ , $f^{(4)}(1)=-6$ : $c_1=\dfrac{-6}{8}-\dfrac{5\cdot4}{24\cdot(-1)}=-0.75+0.833=\mathbf{+\tfrac{1}{12}}$ , so $\Gamma(z+1)\sim\sqrt{2\pi z}\,z^z e^{-z}(1+\tfrac{1}{12z}+\cdots)$ .

Step 5. Put in numbers at $z=10$ . The leading term gives $\sqrt{20\pi}\,(10/e)^{10}=\mathbf{3.599\times10^6}$ , which is $0.83\%$ below the exact $10!$ ; multiplying by $1+1/120$ gives $\mathbf{3.6287\times10^6}$ , now within $\mathbf{0.003\%}$ of $3\,628\,800$ . One Gaussian plus one correction reproduces a seven-digit factorial to three decimal places.

2 · The modified Bessel function I₀(λ) at large argument (fully worked)

Use Laplace's method on the integral representation $I_0(\lambda)=\dfrac{1}{2\pi}\int_{-\pi}^{\pi} e^{\lambda\cos\theta}\,\mathrm{d}\theta$ to find its large- $\lambda$ asymptotics, and test at $\lambda=5$ against the exact $I_0(5)=27.2399$ .

Step 1. Identify the exponent and its maximum. Here $f(\theta)=\cos\theta$ , maximal at the interior point $\theta_0=0$ with $f(0)=1$ , $f'(0)=0$ , and $f''(0)=-1$ , so $|f''|=1$ . The prefactor is $g\equiv1/(2\pi)$ .

Step 2. Apply the leading formula. Equation (2) gives $\int_{-\pi}^{\pi}e^{\lambda\cos\theta}\,\mathrm{d}\theta\sim e^{\lambda}\sqrt{2\pi/\lambda}$ , so $I_0(\lambda)\sim \dfrac{1}{2\pi}e^{\lambda}\sqrt{2\pi/\lambda}=\mathbf{\dfrac{e^{\lambda}}{\sqrt{2\pi\lambda}}}$ — the standard asymptotic of $I_0$ .

Step 3. Compute the correction. Expanding $\cos\theta=1-\theta^2/2+\theta^4/24-\cdots$ gives $a=1$ , $b=f'''(0)=0$ , $c=f^{(4)}(0)=+1$ , so $c_1=c/(8a^2)+5b^2/(24a^3)=\mathbf{\tfrac18}$ , and $I_0(\lambda)\sim \dfrac{e^{\lambda}}{\sqrt{2\pi\lambda}}\left(1+\dfrac{1}{8\lambda}+\cdots\right)$ .

Step 4. Put in numbers at $\lambda=5$ . The leading term is $e^5/\sqrt{10\pi}=\mathbf{26.48}$ , $2.8\%$ below the exact $27.2399$ ; times $1+1/40$ gives $\mathbf{27.14}$ ( $0.36\%$ low), and the next term $9/(128\lambda^2)$ brings it to $27.22$ ( $0.09\%$ ). The successive partial sums close in on the true value — the signature of a good asymptotic series at moderate $\lambda$ .

Physics of the portrait

Laplace's method is the mathematical distillate of a physical instinct: in a system with an overwhelmingly large parameter, the whole answer is carried by the most probable configuration plus Gaussian jitter around it — and it traces to a single figure who was as much a physicist of the heavens as a mathematician of chance.

Origins · 1 minPierre-Simon Laplace and how this was discovered

Pierre-Simon Laplace

1749–1827 · celestial mechanics, probability theory, and the birth of asymptotic analysis

Pierre-Simon Laplace, the "Newton of France", spent his life showing that the visible order of the world was the overwhelmingly probable outcome of underlying law. In a 1774 memoir on the probability of causes he faced integrals $\int e^{\lambda f(x)}\,\mathrm{d}x$ with $\lambda$ large, arising when many independent observations are combined, and introduced the device that now bears his name: expand the exponent about its maximum, keep the quadratic term, integrate the resulting Gaussian. In the monumental Théorie analytique des probabilités of 1812 it became a systematic tool for the asymptotics of the theory of errors — the same analysis that underlies the central limit theorem he did so much to establish. The method was not born in a vacuum: Stirling and de Moivre had already found, around 1730, that $n!\approx\sqrt{2\pi n}\,(n/e)^n$ , the prototype result his technique reproduces in a few strokes. Its later life was richer still. Peter Debye lifted the idea into the complex plane as the method of steepest descent in 1909; G. N. Watson gave the rigorous endpoint lemma in 1918; and in the twentieth century the same Gaussian-peak philosophy became the saddle-point evaluation of partition functions and the stationary-phase backbone of semiclassical quantum mechanics. What Laplace built to tame the arithmetic of chance turned out to be the universal grammar of "large parameter, sharp peak".

History

1730

— Stirling and de Moivre establish n! ≈ √(2πn)(n/e)ⁿ, the prototype asymptotic

1774

— Laplace approximates large-parameter integrals by expanding about the maximum

1812

— Laplace's Théorie analytique des probabilités systematises the method

1909

— Debye recasts the idea in the complex plane as the method of steepest descent

1918

— Watson's lemma gives the rigorous endpoint and Laplace-transform version

Going further

Texts. Bender & Orszag, Advanced Mathematical Methods for Scientists and Engineers, §6.3–6.4 (the definitive treatment of Laplace's method and its corrections) · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 12 (asymptotic series) · Stone & Goldbart, Mathematics for Physics, ch. on asymptotic expansions · Hassani, Mathematical Physics, ch. on asymptotic methods · Byron & Fuller, Mathematics of Classical and Quantum Physics, on integral asymptotics · Olver, Asymptotics and Special Functions (rigorous error bounds). Key numerical checks: (1) the Gaussian moment identity $\int t^{2n}e^{-at^2/2}\,\mathrm{d}t=(2n-1)!!\,a^{-n}\sqrt{2\pi/a}$ ; (2) Stirling with $c_1=+1/12$ — the sign is positive because the middle term of equation (6) has $(f'')^3<0$ — reproducing $10!$ to $0.003\%$ at $z=10$ ; (3) the endpoint scaling $\sim1/\lambda$ consistent with Watson's lemma; (4) the Bessel asymptotic $I_0(\lambda)\sim e^\lambda/\sqrt{2\pi\lambda}\,(1+1/8\lambda)$ , giving $27.14$ versus the exact $27.2399$ at $\lambda=5$ .

Builds on: Order Symbols and Gauge Functions · Asymptotic Series and Poincaré Expansions · Regular Perturbation Theory. Leads to: Watson's Lemma · Method of Steepest Descent and Saddle Points · Method of Stationary Phase · WKB Approximation.

Ask the tutor: Bring me any integral with a large parameter — a partition function, a Bessel or gamma asymptotic, a tunnelling amplitude — and we will find its dominant peak and read off the leading behaviour together.

When exactly is Laplace's method the right tool, as opposed to steepest descent or stationary phase?

Use Laplace's method when the exponent $\lambda f(x)$ is real and has a genuine maximum on a real interval, with $\lambda\to+\infty$ . If the exponent is purely imaginary, $e^{i\lambda\phi(x)}$ , the integrand oscillates rather than peaks and you need the method of stationary phase, whose dominant points are where $\phi'=0$ and which carries a phase $e^{\pm i\pi/4}$ . If the exponent is complex, or the real maximum is inaccessible on the given contour, you deform the contour through a saddle point of $f$ — the method of steepest descent, of which Laplace's real-axis case is the special instance. In short: real peak is Laplace; oscillation is stationary phase; complex contour is steepest descent.

Why is the correction series divergent, and how can a divergent series be useful?

The coefficients $c_n$ are built from ever-higher Gaussian moments, which grow factorially, so $c_n/\lambda^n$ eventually increases without bound no matter how large $\lambda$ is. The series is asymptotic in Poincaré's sense: for a fixed truncation order $N$ , the error is of order the first omitted term and shrinks like $\lambda^{-(N+1)}$ as $\lambda\to\infty$ . In practice you truncate near the optimal order, roughly where consecutive terms stop shrinking ( $n\sim\lambda$ ), and there the accuracy is exponentially good; summing past that point only makes things worse. To go further you need resummation — Borel summation reconstructs a function whose asymptotic expansion is the given series, with resurgence tying the ambiguities to exponentially small contributions from other saddles.

How do I know I have found the global maximum and not just a local one?

This is the assumption most often violated in practice, and getting it wrong produces answers that are not merely inaccurate but exponentially wrong. The safeguard is to evaluate $f$ at every candidate — all interior stationary points and both endpoints — and compare the values of $f$ itself, not of $f''$ or $g$ , because the leading exponential $e^{\lambda f}$ swamps everything else. The largest value of $f$ wins outright; a competitor lower by a finite gap $\Delta$ is suppressed by $e^{-\lambda\Delta}$ , and only exactly tied maxima are summed together. When $f$ depends on parameters, watch for a value where a subdominant maximum overtakes the leading one — a Stokes phenomenon, a discontinuous change in which saddle dominates.

Where does the √(2π) actually come from — is it fundamental or a coincidence?

It is the area under a unit Gaussian, $\int_{-\infty}^\infty e^{-t^2/2}\,\mathrm{d}t=\sqrt{2\pi}$ , and it is as fundamental as the Gaussian's role as the universal shape of a smooth peak. Any function with a smooth interior maximum looks locally like an inverted parabola in its exponent, and exponentiating an inverted parabola gives a Gaussian; the constant $\sqrt{2\pi}$ is the once-and-for-all price of that shape. This is why the same $\sqrt{2\pi}$ turns up in Stirling's formula, in the normal distribution, and in the fluctuation determinant of a path integral — all the same integral in different clothes. Change the curvature and you rescale the width by $1/\sqrt{|f''|}$ , but the $\sqrt{2\pi}$ is invariant.

What happens if the second derivative vanishes at the maximum?

Then the peak is flatter than a parabola and the standard formula breaks down, because its width estimate $1/\sqrt{\lambda|f''|}$ diverges. You must keep the lowest nonvanishing derivative, say $f^{(2p)}(x_0)<0$ with $2p\ge4$ , so that locally $f(x)-f(x_0)\approx f^{(2p)}(x_0)(x-x_0)^{2p}/(2p)!$ . The integral becomes a generalised Gaussian $\int e^{-\beta u^{2p}}\,\mathrm{d}u$ that evaluates to a gamma function, and the scaling softens from $\lambda^{-1/2}$ to $\lambda^{-1/(2p)}$ — a flatter peak is wider and contributes more. In the complex-plane version this higher-order critical point is exactly what produces Airy-function (rather than Gaussian) behaviour near a coalescence of saddles.

Does the prefactor g(x) ever change the leading power of λ?

Only in the exceptional case that $g$ vanishes at the maximum. Generically $g(x_0)\neq0$ , so freezing $g$ at its peak value is exact to leading order and the amplitude sets only the overall constant, not the power of $\lambda$ . But if $g(x_0)=0$ the leading Gaussian contribution disappears and the first surviving term comes from $g''(x_0)$ times the second moment $\langle t^2\rangle=1/(\lambda|f''|)$ , costing an extra factor of $1/\lambda$ ; a nonzero $g'(x_0)$ contributes only an odd moment that vanishes, so you still land on the $g''$ term. The rule of thumb: a nonzero amplitude shifts the constant, a vanishing one demotes the leading power.