The Stirling Approximation

Write out $23!$ and you get a 23-digit number; write out $(6\times10^{23})!$ — the factorial that shows up the moment you count the ways a mole of gas can arrange itself — and the number of digits alone would fill a library. Statistical mechanics runs on exactly these monsters, because the multiplicity $\Omega$ of a macrostate is almost always a ratio of factorials, and entropy is its logarithm, $S=k_B\ln\Omega$ . To do any physics you must tame $\ln N!$ for $N\sim10^{23}$ , and you cannot do it by multiplying $N$ terms. Stirling's approximation is the escape hatch. It trades the jagged, discrete $\ln N!$ for a smooth, differentiable function of a continuous variable: $\ln N!\approx N\ln N-N$ , with a refinement $\ln N!\approx N\ln N-N+\tfrac12\ln(2\pi N)$ that adds the celebrated $\sqrt{2\pi N}$ prefactor. The first form is all thermodynamics usually needs, because entropy is extensive — it grows like $N$ — so the leading $N\ln N-N$ dominates and the correction is a vanishing crumb per particle. The second form is what you keep when you care about fluctuations, normalization constants, or a comparison against an exact small-system answer. The formula is not a lucky guess. Two independent routes deliver it: an elementary one that replaces the sum $\sum_{k=1}^N\ln k$ by the area under $\ln x$ , and a rigorous one that evaluates the Gamma integral $N!=\int_0^\infty x^N e^{-x}\,dx$ by Laplace's saddle-point method, which manufactures the $\sqrt{2\pi N}$ out of a Gaussian peak of width $\sqrt N$ . Both agree, and both expose why it works: a product of many terms, viewed through a logarithm, becomes a well-behaved integral. The reward is that entropy becomes an analytic function you can differentiate — the single move that produces the Boltzmann, Fermi–Dirac, and Bose–Einstein distributions.

Key takeaways

Stirling's approximation replaces the intractable $\ln N!$ with the smooth function $N\ln N-N$ , valid for large $N$ — the workhorse behind every entropy and distribution derivation.
The sharper form $N!\sim\sqrt{2\pi N}\,(N/e)^N$ adds the subleading $\tfrac12\ln(2\pi N)$ correction to the logarithm; keep it only for prefactors and fluctuations.
It follows two ways: approximate $\sum_{k=1}^N\ln k$ by $\int_1^N\ln x\,dx$ , or evaluate the Gamma integral $N!=\int_0^\infty x^N e^{-x}dx$ by Laplace's method.
Because entropy is extensive ( $S\sim N$ ), the leading term survives the thermodynamic limit and the $\ln N$ correction is negligible per particle.
The relative error of the leading form falls as $(\ln N)/(2N)$ ; restoring the $\tfrac12\ln(2\pi N)$ term drops it one whole power of $N$ faster, to $O(1/N)$ .
Key idea: taking a logarithm turns a product of $N$ terms into an integral you can actually do — $\ln N!\approx N\ln N-N$ — making entropy a smooth, differentiable function of particle number.

The core insight

Replace the discrete, intractable $\ln N!$ by the smooth function $N\ln N-N$ (plus a $\tfrac12\ln 2\pi N$ correction); the relative error vanishes as $1/N$ , so in the thermodynamic limit entropy becomes an analytic, differentiable function of particle and occupation numbers.

Where this lives today

Stirling's approximation is not a museum piece; it is load-bearing infrastructure in every corner of modern statistical and information science. In machine learning, the entropy and cross-entropy of categorical distributions — the loss functions that train classifiers — are Stirling-smoothed factorial counts in disguise; the same $-\sum x\ln x$ that Stirling produces from a binomial is the Shannon entropy minimized in every softmax layer. In computational chemistry and biophysics, partition-function estimates for polymers and lattice models lean on Stirling to convert configurational counts into free energies. And in numerical analysis the Lanczos and Spouge approximations to the Gamma function — the routines your calculator and every scientific library call to evaluate $\Gamma(z)$ — are engineered refinements of exactly the de Moivre–Stirling asymptotics derived here. The formula also anchors the theory of large deviations, the modern language for rare events: the rate functions that quantify the exponentially small probability of an atypical fluctuation are Stirling-derived entropies, and they underlie everything from risk modeling to the fluctuation theorems of non-equilibrium thermodynamics.

The physical picture

Here is the whole idea in one image. The logarithm of a factorial is a sum of $\ln k$ from $k=1$ to $N$ — a staircase of bars. The function $\ln x$ is smooth and slowly varying, so the total area of those bars is very nearly the area under the curve $\ln x$ from $1$ to $N$ . That area is $N\ln N-N+1$ , and the leading two terms are Stirling. The approximation is good precisely because $\ln$ changes gently from one integer to the next: the bars hug the curve.

Figure

The sum of ln k tracks the area under ln x — the origin of N ln N − N

The second picture explains the mysterious $\sqrt{2\pi N}$ and shows the approximation's quality with numbers on the plot. The Gamma integral $N!=\int_0^\infty x^N e^{-x}\,dx$ has an integrand $x^N e^{-x}$ that is a sharp bump centered at $x=N$ with width $\sqrt N$ . As $N$ grows the bump becomes needle-thin relative to its position, so replacing it by a Gaussian captures essentially all the area — Laplace's method — and the Gaussian's normalization is precisely $\sqrt{2\pi N}$ . The payoff figure is the punchline of this whole article: the relative error of the two Stirling forms plotted against $N$ . The leading form's error slides down a shallow slope ( $\sim(\ln N)/2N$ ); the moment you restore the $\tfrac12\ln(2\pi N)$ term, the error drops onto a steeper line one full power of $N$ lower.

Figure

Punchline: adding the ½ln(2πN) term buys one whole extra power of N in accuracy

Why physicists usually throw the $\sqrt{2\pi N}$ away: entropy per particle is $s=S/N=k_B\big(\ln N-1+\tfrac{1}{2N}\ln 2\pi N+\dots\big)$ . The correction term $\tfrac{1}{2N}\ln 2\pi N$ marches to zero as $N\to\infty$ , so it is genuinely irrelevant to bulk thermodynamics — it resurfaces only when finite size or fluctuations are the point.

The core equations

The headline result is the leading Stirling formula. It is the approximation you reach for by reflex whenever a factorial appears inside a logarithm, and it is exact in the sense that the ratio of the two sides tends to $1$ .

\ln N! \approx N\ln N - N \qquad (1)

where: $N$ — a large positive integer (typically a particle or occupation number); $N!=N\cdot(N-1)\cdots 2\cdot 1$ ; $\ln$ — natural logarithm. Read it as: the log of a factorial grows like $N\ln N$ , faster than linearly but far slower than the factorial itself. Keeping one more term restores the prefactor and turns the approximation from "good" into "startlingly good even at modest $N$ ". This is the de Moivre–Stirling form, the one to quote when you need $N!$ itself rather than just its logarithm.

N! \sim \sqrt{2\pi N}\,\left(\frac{N}{e}\right)^{N} \qquad (2)

where: $\sim$ — asymptotic equality (ratio $\to1$ as $N\to\infty$ ); $e=2.71828\ldots$ ; the $\sqrt{2\pi N}$ prefactor is the width of the Gaussian peak of the Gamma integrand. Taking the logarithm gives $\ln N!=N\ln N-N+\tfrac12\ln(2\pi N)$ . The full story is an asymptotic series in $1/N$ , whose successive terms are governed by the Bernoulli numbers. It is not a convergent series — its terms eventually grow — but truncated at its smallest term it is exquisitely accurate.

\ln N! = N\ln N - N + \tfrac{1}{2}\ln(2\pi N) + \frac{1}{12N} - \frac{1}{360 N^{3}} + \cdots \qquad (3)

where: $\tfrac{1}{12N}$ — first Bernoulli correction; the alternating higher terms carry Bernoulli numbers $B_{2j}$ . For $N\sim10^{23}$ everything past $N\ln N-N$ is noise. The whole edifice rests on the integral representation of the factorial — the Gamma function — which is what makes Laplace's method available and extends $N!$ to non-integer arguments.

\Gamma(N+1) = N! = \int_{0}^{\infty} x^{N} e^{-x}\,\mathrm{d}x \qquad (4)

where: $\Gamma$ — Euler's Gamma function; $x$ — a dummy variable over the positive reals. The integrand $x^Ne^{-x}$ peaks at $x=N$ ; that peak is the seed of the saddle-point evaluation. Finally, the single feature that makes Stirling indispensable rather than merely convenient: its derivative is clean. Treating $N$ as continuous,

\frac{\mathrm{d}}{\mathrm{d}N}\ln N! \approx \ln N \qquad (5)

where: the chain-rule term $N\cdot\tfrac1N=1$ exactly cancels the derivative of $-N$ . This cancellation is why maximizing $\ln\Omega=\ln\frac{N!}{\prod_i n_i!}$ over occupation numbers $n_i$ produces $\ln n_i$ terms — and hence the exponential equilibrium distributions.

Derivation AThe Gamma integral by Laplace's method

The goal here is to get the full Stirling formula, prefactor and all, from a single controlled approximation rather than from hand-waving about sums and integrals. The strategy is Laplace's method: when an integral is dominated by a sharp peak, expand the exponent to second order about its maximum and do the resulting Gaussian integral exactly. Applied to the Gamma integral $N!=\int_0^\infty x^Ne^{-x}dx$ , the peak sits at $x=N$ , its curvature sets a width $\sqrt N$ , and the Gaussian normalization is exactly the $\sqrt{2\pi N}$ we are chasing. Watch the $2\pi$ appear from nothing but the width of a bump.

Result.

N! \sim \sqrt{2\pi N}\,\left(\frac{N}{e}\right)^{N}

Step 1. Write the integrand as one exponential: $x^Ne^{-x}=e^{f(x)}$ with $f(x)=N\ln x-x$ . Why: Laplace's method needs the whole integrand in the form $e^{f}$ so that a large parameter multiplies a smooth exponent — here $N$ is large and $f$ scales with it.

Step 2. Locate the maximum: $f'(x)=\tfrac{N}{x}-1=0$ gives $x_0=N$ . Why: for large $N$ the integrand is overwhelmingly concentrated at its single interior maximum; everything else contributes negligibly.

Step 3. Measure the curvature: $f''(x)=-\tfrac{N}{x^2}$ , so $f''(N)=-\tfrac1N$ . Why: the second derivative fixes how fast the peak falls off — a curvature $-1/N$ means a Gaussian of variance $N$ , i.e. width $\sqrt N$ .

Step 4. Expand $f$ to second order about $x_0=N$ :

f(x)\approx f(N)-\frac{1}{2N}(x-N)^2=(N\ln N-N)-\frac{(x-N)^2}{2N}

Why: near a smooth maximum the exponent is parabolic; the constant piece $N\ln N-N$ pulls out of the integral, leaving a pure Gaussian.

Step 5. Do the Gaussian integral, extending the lower limit to $-\infty$ (legitimate because the peak is sharp):

N!\approx e^{N\ln N-N}\int_{-\infty}^{\infty}e^{-(x-N)^2/2N}\,\mathrm{d}x=\left(\frac{N}{e}\right)^{N}\sqrt{2\pi N}

Why: $\int e^{-u^2/2N}du=\sqrt{2\pi N}$ is the standard Gaussian result — the $2\pi$ is its normalization and the $N$ is the variance set in Step 3. The prefactor was never mysterious: it is the area of a bell curve. *Verified: taking the logarithm reproduces $\ln N!=N\ln N-N+\tfrac12\ln(2\pi N)$ ; a numpy check of $\sqrt{2\pi N}(N/e)^N/N!$ gives $0.983,\,0.996,\,0.999$ at $N=5,20,100$ , matching Abramowitz & Stegun 6.1.41.*

Go deeper · Derivation B · 1 minThe error orders: what the prefactor buys you

Derivation BThe error orders: what the prefactor buys you

This derivation is the article's owned result: a precise accounting of how good each Stirling form is, expressed as the power of $N$ at which the relative error decays. It matters because "Stirling is accurate for large $N$ " is too vague to use — you need to know whether dropping the $\tfrac12\ln(2\pi N)$ term costs you a factor of $N$ or a factor of $\ln N$ . The answer, confirmed by a machine table over six decades of $N$ , is sharp: the leading form errs at order $(\ln N)/N$ , and restoring the prefactor buys a clean extra power of $N$ , taking the error to $O(1/N)$ relative to $\ln N!$ .

Result.

\frac{|\,N\ln N-N-\ln N!\,|}{\ln N!}\sim\frac{\tfrac12\ln(2\pi N)}{N\ln N},\qquad \frac{|\,N\ln N-N+\tfrac12\ln(2\pi N)-\ln N!\,|}{\ln N!}\sim\frac{1}{12N\ln N}

Step 1. The absolute error of the leading form is exactly the dropped terms: $\ln N!-(N\ln N-N)=\tfrac12\ln(2\pi N)+\tfrac{1}{12N}+\cdots\approx\tfrac12\ln(2\pi N)$ . Why: subtract (1) from the full series (3); the first surviving term is the prefactor's logarithm.

Step 2. Divide by $\ln N!\approx N\ln N$ to get the relative error $\approx\tfrac{\ln(2\pi N)}{2N\ln N}\sim\tfrac{1}{2N}$ (since $\ln 2\pi N\approx\ln N$ for large $N$ ). Why: relative error is what governs how many correct digits you keep; the $1/2N$ scaling means one extra correct digit per factor $\sim3$ in $N$ .

Step 3. The absolute error of the full form is the next dropped term, $\tfrac{1}{12N}$ ; dividing by $\ln N!\approx N\ln N$ gives relative error $\sim\tfrac{1}{12N^2\ln N}$ — down by a factor $N$ versus Step 2. Why: the prefactor term was the entire $O(1/N)$ -relative deficit; removing it exposes the much smaller $O(1/N^2)$ -relative Bernoulli term.

Step 4. Confirm on a table. At $N=10,10^2,\dots,10^6$ the leading relative error is $1.38\times10^{-1},\,8.86\times10^{-3},\,7.40\times10^{-4},\dots$ (ratio to $\tfrac12\ln(2\pi N)/\ln N!$ equals $1.00$ to three figures), while the full-form relative error is $5.5\times10^{-4},\,2.3\times10^{-6},\,1.4\times10^{-8},\dots$ — each one power of $N$ smaller. Why: the numbers march down the two predicted slopes exactly, pinning the error orders and vindicating the "drop the prefactor in bulk, keep it for finite size" rule. *Verified: `verify_aq-sm-stirling-approximation.py` reproduces both error-order columns and the ratio-to-prediction $\to1$ across six decades of $N$ .*

Limiting cases

The right truncation depends entirely on how large $N$ is and what you are computing:

Case	Result	Meaning
Thermodynamic limit	$\ln N!\to N\ln N-N$	Leading term suffices; relative error $\sim\tfrac{\ln N}{2N}\to0$ . The default in deriving entropy and distributions.
Finite-size / fluctuations	$\ln N!\approx N\ln N-N+\tfrac12\ln(2\pi N)$	Keep the prefactor for normalization constants, Gaussian fluctuation amplitudes, and comparisons to exact small- $N$ answers.
High precision	full series $+\tfrac{1}{12N}-\tfrac{1}{360N^3}+\cdots$	Divergent (Gevrey-1) but, truncated at its smallest term, exponentially accurate.
Small $N\lesssim1$ — breakdown	use exact factorials	$0!=1!=1$ , yet $N\ln N-N$ gives $0$ and $-1$ ; the approximation is simply wrong here.

Worked examples

1 · Watch the gap close: ln(100!) by hand-scale reasoning (confidence builder)

Compute $\ln(100!)$ exactly, then with each Stirling form, and see how much the prefactor buys.

Step 1. Exact value (from $\ln\Gamma(101)$ ): $\ln(100!)=\mathbf{363.739}$ . This is the truth we are approximating.

Step 2. Leading Stirling: $100\ln 100-100=100(4.60517)-100=460.517-100=\mathbf{360.517}$ . The error is $363.739-360.517=3.222$ , a relative error of $0.886\%$ .

Step 3. Add the prefactor: $\tfrac12\ln(2\pi\cdot100)=\tfrac12\ln(628.3)=\tfrac12(6.443)=3.222$ . So the full form gives $360.517+3.222=\mathbf{363.739}$ — matching the exact value to three decimals.

Step 4. Read the lesson: the entire leading-form deficit at $N=100$ is $3.222$ , and it is exactly the prefactor term $\tfrac12\ln(2\pi N)$ . That is the owned error-order result made concrete — the correction you dropped is the whole error, to leading order. Contrast with $N=10$ : exact $\ln 10!=15.104$ , leading form $10\ln10-10=13.026$ (error $2.08$ , a full $14\%$ ), and adding $\tfrac12\ln(2\pi\cdot10)=2.07$ recovers $15.10$ . Small $N$ is where the prefactor earns its keep.

2 · The entropy of a spin-½ paramagnet from a binomial (exam level)

A crystal has $N=6.0\times10^{23}$ spin-½ atoms. Using Stirling, find the entropy when a fraction $x=N_\uparrow/N=1/4$ point up, and compare it to the maximum.

Step 1. The multiplicity is the binomial $\Omega=\binom{N}{N_\uparrow}$ — the number of ways to choose which atoms are up. Its log is $\ln\Omega=\ln N!-\ln N_\uparrow!-\ln(N-N_\uparrow)!$ .

Step 2. Apply the leading Stirling to each factorial. The $-N$ pieces cancel (since $N_\uparrow+(N-N_\uparrow)=N$ ), leaving the Shannon form $\ln\Omega=-N\big[x\ln x+(1-x)\ln(1-x)\big]$ .

Step 3. Plug in $x=1/4$ : $x\ln x=0.25\ln0.25=0.25(-1.386)=-0.3466$ and $(1-x)\ln(1-x)=0.75\ln0.75=0.75(-0.2877)=-0.2158$ . The bracket is $-0.5623$ , so $\ln\Omega=0.5623\,N$ and $S=k_B\ln\Omega=\mathbf{0.562\,Nk_B}$ .

Step 4. The maximum entropy is at $x=1/2$ : $S_{\max}=Nk_B\ln2=\mathbf{0.693\,Nk_B}$ . So the quarter-up state carries $0.562/0.693=81\%$ of the maximum entropy — still enormous, because $\ln\Omega$ is flat near its peak. Only Stirling makes this a one-line calculation: without it, differentiating $\binom{N}{N_\uparrow}$ to find the peak would be hopeless.

Step 5. Sanity check the leading-order use: at $N\sim10^{23}$ the dropped prefactor $\tfrac12\ln[2\pi x(1-x)N]$ is about $\tfrac12\ln(10^{23})\approx27$ , utterly negligible against $\ln\Omega\sim10^{23}$ . The extensive entropy is safe.

Origins · 1 minJames Stirling and how this was discovered

James Stirling

1692–1770 · Scotland and London

Stirling was a Scottish mathematician — nicknamed "the Venetian" for years spent in Italy — whose 1730 book Methodus Differentialis studied the asymptotics of series and factorials. The formula that bears his name, $\ln N!\approx N\ln N-N+\tfrac12\ln(2\pi N)$ , has a tangled paternity: Abraham de Moivre had already found that $N!$ grows like $C\,\sqrt N\,(N/e)^N$ while working on the normal approximation to the binomial, but it was Stirling who identified the constant as $C=\sqrt{2\pi}$ . The two corresponded and effectively share the result — which is why the sharp form is often called the de Moivre–Stirling formula. It is a rare case of a piece of pure eighteenth-century asymptotics becoming, two centuries later, the single most-used approximation in statistical physics: every time an entropy is differentiated, Stirling is silently at work.

History

1730

— Stirling publishes Methodus Differentialis, giving the asymptotic expansion of ln N! and fixing the prefactor constant as √(2π).

1730

— de Moivre, studying the normal approximation to the binomial, independently establishes the √N (N/e)^N growth; the sharp formula is jointly the de Moivre–Stirling result.

1774

— Laplace develops the saddle-point ("Laplace's method") that rigorously produces the full Stirling series from the Gamma integral.

1877

— Boltzmann links entropy to microstate counting, S = k ln W, making ln N! the central object of statistical mechanics and Stirling its indispensable tool.

1886

— Poincaré founds the theory of asymptotic (divergent-but-useful) series, opening it with Stirling's series and its Euler–Maclaurin/Bernoulli-number terms.

Going further

Texts. Pathria & Beale, Statistical Mechanics, 3rd ed., Appendix B (Stirling's formula) · Kardar, Statistical Physics of Particles, Ch. 2–4 · Reif, Fundamentals of Statistical and Thermal Physics, Appendix A.6 · Sethna, Statistical Mechanics: Entropy, Order Parameters, and Complexity, Ch. 3 · Abramowitz & Stegun, Handbook of Mathematical Functions, 6.1.37–6.1.42.

Provenance. Leading and full Stirling formulas, the asymptotic series with Bernoulli corrections (Abramowitz & Stegun 6.1.37–6.1.42), the Laplace/saddle-point derivation ( $x_0=N$ , $f''=-1/N$ ), and the Euler–Maclaurin constant $\tfrac12\ln2\pi$ are standard. A verifier should reproduce: the ratio $\sqrt{2\pi N}(N/e)^N/N!\to1$ ( $0.983,0.996,0.999$ at $N=5,20,100$ ); the owned error-order table (leading relative error $\sim\tfrac12\ln(2\pi N)/\ln N!$ , full relative error one power of $N$ smaller, over $N=10$ – $10^6$ ); the worked numbers $\ln100!=363.739$ , leading $360.517$ , full $363.739$ , and $\ln10!=15.104$ ; and the paramagnet entropies $0.562\,Nk_B$ ( $x=1/4$ ) vs $0.693\,Nk_B$ ( $x=1/2$ ). Constant $k_B=1.380649\times10^{-23}\,\mathrm{J/K}$ (CODATA, exact since the 2019 SI). No non-standard numbers are asserted.

Prerequisites. Microstates and Macrostates · Boltzmann Entropy · Logarithms and integration. Leads to. Boltzmann Entropy · Statistical Ensembles · Accessible Phase Space and Energy Shells.

Is Stirling's approximation exact, or just "good enough"?

Neither, exactly. It is asymptotically exact: the ratio of the two sides of $N!\sim\sqrt{2\pi N}(N/e)^N$ tends to $1$ as $N\to\infty$ , so for the astronomical $N$ of statistical mechanics it is exact for every practical purpose. But the full series $N\ln N-N+\tfrac12\ln(2\pi N)+\tfrac{1}{12N}-\cdots$ is divergent — if you keep summing its terms they eventually blow up. The trick is that you truncate it at its smallest term, where the leftover error is exponentially tiny. So it is a controlled, quantifiable approximation whose error you can bound, not a vague heuristic.

Why do we take the logarithm before approximating, rather than approximating N! directly?

Two reasons. First, $N!$ itself is hyper-astronomical — for a mole it has more digits than atoms in your body — so there is no useful "smooth curve" to fit it to. Its logarithm, by contrast, is a tame $\sim N\ln N$ , an ordinary extensive quantity. Second, physics wants the logarithm anyway: entropy is $S=k_B\ln\Omega$ , and $\Omega$ is built from factorials. Taking the log turns the product-of-factorials multiplicity into a sum, and Stirling turns that sum into a differentiable function — exactly what you need to maximize $\ln\Omega$ and find equilibrium.

When can I safely drop the ½ln(2πN) term?

Whenever you are computing an extensive quantity in the thermodynamic limit — which is most of bulk thermodynamics. Entropy scales as $S\sim N$ , while the prefactor term contributes only $\sim\ln N$ , so per particle it is $\tfrac{\ln N}{2N}\to0$ . Keep it in three situations: (1) fluctuations and finite-size corrections, where the $\sqrt N$ width is the whole point; (2) normalization constants of probability distributions; (3) any comparison against an exact small- $N$ result. And whenever you keep it for one factorial, keep it for all of them in the same calculation.

Where does the 2π come from? It seems to appear out of nowhere.

It is the normalization of a Gaussian. The Gamma integrand $x^Ne^{-x}$ is a bump peaked at $x=N$ ; near the peak it looks like a bell curve $e^{-(x-N)^2/2N}$ of variance $N$ . The area under that bell curve is $\int e^{-u^2/2N}du=\sqrt{2\pi N}$ , and the $2\pi$ is simply the constant in the Gaussian integral $\int e^{-u^2/2}du=\sqrt{2\pi}$ . So the $2\pi$ that shows up in factorial asymptotics is the same $2\pi$ from the normal distribution — no coincidence, since de Moivre found the formula while studying exactly that distribution.

The derivative d/dN(ln N!) ≈ ln N looks too clean. Why does the "+1" disappear?

Differentiate $N\ln N-N$ : the product rule on $N\ln N$ gives $\ln N+N\cdot\tfrac1N=\ln N+1$ , and the derivative of $-N$ is $-1$ , so they cancel to leave exactly $\ln N$ . That cancellation is not an accident — it is the reason Stirling is so beloved for variational problems. When you maximize $\ln\Omega=\ln N!-\sum_i\ln n_i!$ over the occupation numbers $n_i$ , each derivative is just $-\ln n_i$ , and setting the Lagrange-constrained variation to zero gives $\ln n_i=\alpha-\beta\epsilon_i$ , i.e. $n_i\propto e^{-\beta\epsilon_i}$ — the Boltzmann distribution falls out in one line.

Is Stirling's formula actually due to Stirling?

Only partly. Abraham de Moivre, working on the normal approximation to the binomial around 1730, had already shown that $N!\approx C\sqrt N(N/e)^N$ for some constant $C$ — he just could not pin down $C$ . James Stirling supplied the missing piece, identifying $C=\sqrt{2\pi}$ . The two corresponded, and the sharp form is fairly called the de Moivre–Stirling formula. Stirling's own contribution appears in his 1730 Methodus Differentialis; the rigorous saddle-point derivation you see today came later, with Laplace's method.

How can a divergent series be useful? Doesn't divergence mean it's wrong?

A divergent asymptotic series is a different animal from a convergent one. For a fixed $N$ , its terms first shrink — dramatically — reach a minimum, and only then grow. If you stop at that minimum term, the remaining error is smaller than any power of $1/N$ : exponentially small. So for large $N$ you get spectacular accuracy from just two or three terms, even though summing all of them would diverge. This is generic in physics — perturbation series in quantum field theory behave the same way. The rule is: truncate at the smallest term, never "to convergence".