Stirling's approximation by the saddle-point method

Derivation · the atom

Stirling's approximation by the saddle-point method

D-019 stirling-saddle-point Home PU-205 Threads chance Depends on — verified

Statement

N! grows as N^N e^−N √(2πN), obtained by recognising that the integrand of the gamma function is overwhelmingly concentrated near its maximum.

Why it matters

Statistical mechanics is built on ln N! with N ~ 10²³, so this is not an optional refinement. More importantly, the method — expand about the maximum, keep the quadratic term, do the Gaussian integral — is the same technique that later gives the WKB approximation, the classical limit of the path integral, and mean-field theory. Teaching it once, properly, pays four times.

Assumptions

N ≫ 1.The relative error is of order 1/12N. At N = 10 this is under 1%; at N = 10²³ it is unmeasurable.

The integrand has a single sharp maximum.Multiple comparable maxima require summing contributions from each. This is not academic — it is exactly the situation at a first-order phase transition.

The Gaussian tail may be extended to ±∞.The true lower limit is x = 0, but the integrand there is smaller by a factor of order e^−N, so the extension costs nothing.

Derivation

N! = Γ(N+1) = ∫₀^∞ x^N e^−x dx

The integral representation. Verify by parts: it satisfies the same recursion and initial value as the factorial. Converting a discrete product into an integral is what makes calculus available. A

= ∫₀^∞ e^f(x) dx, f(x) = N ln x − x

Write the integrand as a single exponential. Necessary because the two factors individually behave badly — one grows without bound, the other decays — while their combination has a clean maximum. A

f′(x) = N/x − 1 = 0 ⟹ x₀ = N, f(N) = N ln N − N

Locate the maximum. It sits at x = N, and the peak height already gives the leading behaviour. B

f″(x) = −N/x², f″(N) = −1/N

Curvature at the maximum. Negative, confirming it is a maximum, and small — the peak is broad in absolute terms but narrow relative to N, its width scaling as √N against a position of N. B

f(x) ≈ f(N) − (x−N)² / 2N

Taylor expand to second order. The linear term is absent because we are at a stationary point — that is the whole reason for locating it first. B

N! ≈ e^{N ln N − N} ∫_−∞^∞ e^−u²/2N du = N^N e^−N √(2πN)

Substitute u = x − N, extend the limits, and evaluate the Gaussian integral ∫ e^−au² du = √(π/a) with a = 1/2N. C

Result

N! ≈ N^N e^−N √(2πN) ln N! ≈ N ln N − N + ½ ln(2πN)

Reading. In statistical mechanics only the logarithm is needed, and only the first two terms of it: the √(2πN) contributes ½ ln(2πN), which is of order ln N against terms of order N and is discarded without hesitation. But it must be derived before it can be discarded honestly.

Units check. N! and the argument x of the gamma integral are pure numbers, so every term — N^N, e^−N, √(2πN) — is dimensionless, as a count must be. The logarithm takes a dimensionless argument, as required.

Limiting cases

Large N: relative error falls as 1/12N — the next term of the asymptotic series. At N = 100 the factorial is accurate to 0.083%.
Keeping only the peak height gives ln N! ≈ N ln N − N, the form used throughout statistical mechanics; the √ term is a controlled correction.
For the logarithm — the quantity actually used — the accuracy is far better than for the factorial itself.

Breaks when

N is small. At N = 1 the approximation gives 0.922 against a true value of 1.
The exponent has several comparable maxima — each must be summed, and interference between them is the physics of a first-order transition.
The maximum lies at or near a boundary of the integration range, where the symmetric Gaussian extension is illegitimate.

Failure modes

Using ln N! ≈ N ln N and dropping the −N. That term is the same order as the first and dropping it is a gross error, not a refinement.
Keeping √(2πN) in a statistical-mechanics calculation where it is negligible, while elsewhere dropping terms of order N. Inconsistent truncation is worse than either choice made consistently.
Applying the approximation to the small factorials that appear in combinatorial prefactors.

Worked number

At N = 10: exact 10! = 3,628,800; Stirling gives 3,598,696 — a relative error of 0.83%. At N = 5 the error is 1.7% (120 vs 118.02); at N = 100 it is 0.083% (9.333 × 10¹⁵⁷ vs 9.325 × 10¹⁵⁷). The error tracks 1/12N at each point.

Run the check →

Tier 3 · watch out

The same saddle-point argument reappears as WKB (PU-401) and the classical limit of the path integral (PU-401/L07). What is a broad-but-narrow peak here becomes a stationary-phase point there; the Gaussian around it is the fluctuation determinant. Learn the bookkeeping once and the later derivations are the same moves in a different alphabet.

Discussion

The strength of Stirling's formula is that every symbol in it carries a separate physical meaning, and the saddle-point derivation hands each one to you in turn. The factor N^Ne^−N is nothing but the height of the integrand at its peak, e^f(N) with f(N) = N ln N − N. The factor √(2πN) is the width of that peak — the area under a Gaussian of curvature |f″(N)| = 1/N, so its width scales as √N. A factorial that spans hundreds of orders of magnitude is thus reconstructed from just two numbers: how tall the dominant term is, and how many terms of comparable size cluster around it.

This is why the result is structured as (peak height) × (peak width) and why the two pieces are treated so differently in practice. In statistical mechanics one works with ln N! for N ~ 10²³, and there the split becomes a hierarchy: N ln N − N is of order 10²³, while ½ ln(2πN) is of order 50. Keeping the width term in an entropy is like measuring the distance to the Moon and worrying about the thickness of a coat of paint. But — and this is the discipline the page insists on — you may only discard it after you have derived it and confirmed its size. Dropping the −N, by contrast, is not a refinement but a factor-of-e^N catastrophe, because that term is the same order as the one you keep.

The deeper point is that this derivation is a template, not a one-off. It is Laplace's method: any integral of the form ∫e^{N g(x)} dx, with a large parameter N multiplying a smooth function with a single interior maximum, is asymptotically √(2π/N|g″|) e^{N g(x₀)}. Rotate the contour into the complex plane and the same argument becomes the saddle-point (steepest-descent) method, where the maximum is a saddle of the analytic integrand; that is the version that yields WKB, the semiclassical limit ℏ → 0 of the Feynman path integral (the classical path is the stationary point of the action), and the leading asymptotics of Bessel and Airy functions. The 1/12N correction quoted on the page is simply the next term of Watson's lemma — the contribution of the cubic and quartic pieces of f that we truncated at second order.

There is also a probabilistic reading that unifies the two examples below. The gamma integrand x^Ne^−x, once normalised, is the probability density of a sum of N independent exponential waiting times; its mean is N and its variance is N, so by the central limit theorem it approaches a Gaussian of width √N centred on N — which is exactly what the quadratic expansion produces. Stirling's formula and the CLT are two faces of the same concentration-of-measure phenomenon, and the √(2πN) is the normalisation of that limiting Gaussian.

Common misconceptions.

“Stirling says N! ≈ N ln N.” No — that confuses the formula for N! with the leading term of ln N!, and even then omits the essential −N.
“The √(2πN) is negligible, so it never matters.” It is negligible inside a logarithm of order N. In a bare prefactor — the peak height of a binomial, a partition-function ratio — it is often the entire answer (Example 2).
“Bigger N always means safe.” Not if the exponent develops a second comparable maximum (a first-order transition) or if the maximum sits at the boundary of the range, where only half the Gaussian lies inside and an extra factor of ½ appears (Problem 5).
“Apply it everywhere factorials appear.” The small factorials in combinatorial prefactors (a 3! symmetry factor) must be evaluated exactly; Stirling is an asymptotic statement.

Worked examples

Example 1 — Entropy of a two-state paramagnet (the width term is rightly dropped). A crystal has N independent spins, each ↑ or ↓. Find the configurational entropy when a fraction x = 0.75 point up, for one mole of spins.

W = N!n! (N−n)!, S = k ln W, n = xN

Number of ways to place n up-spins among N sites, and the Boltzmann entropy.

ln W ≈ (N ln N − N) − (n ln n − n) − ((N−n) ln(N−n) − (N−n))

Apply ln m! ≈ m ln m − m to each factorial. Every ½ln(2πm) term is dropped — consistently, since all are of order ln N ≪ N.

the −N, +n, +(N−n) cancel; put n = xN: ln W = −N[x ln x + (1−x) ln(1−x)]

The extensive term N ln N cancels identically because n + (N−n) = N; what survives is the Shannon mixing entropy per spin.

S = −R[x ln x + (1−x) ln(1−x)] = −(8.314)[0.75 ln 0.75 + 0.25 ln 0.25]

Use N = N_A so N_Ak = R = 8.314 J K⁻¹ mol⁻¹. Substitute numbers only now.

= −(8.314)[0.75(−0.2877) + 0.25(−1.3863)] = −(8.314)(−0.5623)

0.75 ln 0.75 = −0.2158, 0.25 ln 0.25 = −0.3466; sum −0.5623.

S ≈ 4.68 J K⁻¹ mol⁻¹

Answer. The configurational entropy at 75% alignment is 4.68 J K⁻¹ mol⁻¹ (compare the maximum R ln 2 = 5.76 J K⁻¹ mol⁻¹ at x = ½). Retaining the √(2πN) factors would have changed this by ~10⁻²² of its value — correctly discarded.

Example 2 — Peak of a binomial distribution (the width term is the whole answer). Toss a fair coin N = 10⁴ times. What is the probability of exactly N/2 heads?

P = 12^N N!((N/2)!)²

The central binomial coefficient divided by the total number of sequences 2^N.

N! ≈ N^Ne^−N√(2πN), (N/2)! ≈ (N/2)^N/2e^−N/2√(πN)

Now the full formula is mandatory: this is a bare ratio, not a logarithm. Note √(2π·N/2) = √(πN).

N!((N/2)!)² ≈ N^Ne^−N√(2πN)(N/2)^Ne^−N(πN) = 2^N √(2πN)πN = 2^N√2πN

The peak-height parts cancel: N^N/(N/2)^N = 2^N and the e^−N cancels. Only the width factors remain.

P = 2^−N · 2^N√2πN = √2πN = √2π·10⁴

The exponentially large pieces have vanished completely; the answer is purely the Gaussian-normalisation prefactor.

P(N/2) ≈ √(2/πN) = 7.98 × 10⁻³

Answer. About 0.80%. This equals 1/(σ√(2π)) with standard deviation σ = √N/2 = 50 — the peak of the Gaussian that the binomial becomes. Had we used only N ln N − N we would have obtained the nonsense P = 1, since the leading terms cancel exactly. Here the √(2πN) is not a correction; it is the physics.

Problems

(Warm-up.) Estimate ln(10!) using the three-term formula ln N! ≈ N ln N − N + ½ln(2πN) and compare with the exact value. What does the two-term formula alone give?
Solution
Exact: 10! = 3 628 800, so ln(10!) = 15.104. Two-term: 10 ln 10 − 10 = 23.026 − 10 = 13.026 — an error of 2.08, i.e. the missing ½ln(2πN). Three-term: add ½ln(2π·10) = ½ln(62.83) = ½(4.140) = 2.070, giving 15.096. The residual error is 15.104 − 15.096 = 0.008, in excellent agreement with the next asymptotic term 1/(12N) = 1/120 = 0.0083.
(Order of magnitude.) How many decimal digits does 100! have? Use log₁₀(100!) = ln(100!)/ln 10.
Solution
ln(100!) ≈ 100 ln 100 − 100 + ½ln(200π) = 460.517 − 100 + ½ln(628.3) = 360.517 + 3.221 = 363.739. Divide by ln 10 = 2.302585: log₁₀(100!) = 157.97. The number of digits is ⌊157.97⌋ + 1 = 158, which matches the exact 100! exactly.
(Testing the error law.) Predict the relative error of the leading Stirling formula N! ≈ N^Ne^−N√(2πN) at N = 20, then verify against the exact 20! = 2.4329 × 10¹⁸.
Solution
The next term of the asymptotic series multiplies the formula by (1 + 1/(12N) + …), so the leading formula undershoots by about 1/(12·20) = 1/240 = 0.417%. Check: 20²⁰ = 1.04858 × 10²⁶, e⁻²⁰ = 2.06115 × 10⁻⁹, √(40π) = 11.210. Product = 2.4225 × 10¹⁸. Relative error = (2.4329 − 2.4225)/2.4329 = 0.427% — as predicted.
(Entropy of mixing.) An ideal gas of N molecules, half of species A and half of species B, is arranged on N sites. Using Stirling, show that the mixing entropy is S = Nk ln 2 and evaluate it for one mole.
Solution
The number of arrangements is W = N!/((N/2)!)². With the two-term formula, ln W = (N ln N − N) − 2(½N ln(½N) − ½N) = N ln N − N ln(N/2) = N ln 2. Hence S = k ln W = Nk ln 2. For one mole, S = R ln 2 = 8.314 × 0.6931 = 5.76 J K⁻¹ mol⁻¹. (This is the discrete-lattice route to the −Nk ∑ x_i ln x_i mixing entropy.)
(Method reuse — a boundary saddle.) Find the leading large-N asymptotics of I(N) = ∫₀^π e^{N cos θ} dθ by Laplace's method. Explain why the answer carries a factor ½ relative to a naive interior-maximum estimate, and connect it to the modified Bessel function I₀(N).
Solution
The exponent g(θ) = cos θ is maximal at θ = 0, which lies on the boundary of the range. Expand cos θ ≈ 1 − θ²/2, so g(0) = 1, g″(0) = −1. Then I(N) ≈ e^N ∫₀^∞ e^−Nθ²/2 dθ. Because the maximum sits at the endpoint, only half of the symmetric Gaussian lies inside the range, giving ½√(2π/N):
I(N) ≈ e^N · ½√(2π/N) = e^N√(π/2N)
An interior maximum would have kept the full Gaussian and doubled this — the ½ is precisely the “maximum at a boundary” caveat from the derivation's assumptions. Since I₀(N) = (1/π)∫₀^π e^{N cos θ} dθ, this reproduces the standard result I₀(N) ≈ e^N/√(2πN), the same saddle-point structure that gave Stirling.