Let X1,X2,… be independent, identically distributed random variables with finite mean μ=E[X1] and finite positive variance σ2=Var(X1)>0. For the standardized sum Zn=σn∑i=1nXi−nμ=σn(Xˉn−μ), we have limn→∞P(Zn≤z)=Φ(z)=2π1∫−∞ze−t2/2dt for every z∈R; that is, Zn converges in distribution to N(0,1).
Why is it true?
When you add up many small, independent random effects, the quirks of their individual distributions wash out, and the rescaled fluctuation around the mean universally settles into the bell-shaped Gaussian curve N(0,1) — which is why measurement errors, test scores, and thermal noise all look approximately normal.
Proof sketch
Set Yi=σXi−μ, so E[Yi]=0, Var(Yi)=1, and Zn=n1∑i=1nYi. Let φY(t)=E[eitY1] be the characteristic function of Y1. Since E[Y12]=1<∞, a second-order Taylor expansion at t=0 gives φY(t)=1−2t2+o(t2) as t→0. By independence, the characteristic function of Zn is φZn(t)=[φY(nt)]n=(1−2nt2+o(n1))n→e−t2/2 as n→∞ for each fixed t∈R. Because e−t2/2 is the characteristic function of N(0,1), Lévy's continuity theorem implies ZndN(0,1).