I've been wondering about this one for a while; I find it a little weird how abruptly it happens. Basically, why do we need just three uniforms for $Z_n$ to smooth out like it does? And why does the smoothing-out happen so relatively quickly?
$Z_2$:
$Z_3$:
(images shamelessly stolen from John D. Cook's blog: http://www.johndcook.com/blog/2009/02/12/sums-of-uniform-random-values/)
Why doesn't it take, say, four uniforms? Or five? Or...?
We can take various approaches to this, any of which may seem intuitive to some people and less than intuitive to others. To accommodate such variation, this answer surveys several such approaches, covering the major divisions of mathematical thought--analysis (the infinite and the infinitesimal), geometry/topology (spatial relationships), and algebra (formal patterns of symbolic manipulation)--as well as probability itself. It culminates in an observation that unifies all four approaches, demonstrates there is a genuine question to be answered here, and shows exactly what the issue is. Each approach provides, in its own way, deeper insight into the nature of the shapes of the probability distribution functions of sums of independent uniform variables.
The Uniform $[0,1]$ distribution has several basic descriptions. When $X$ has such a distribution,
The chance that $X$ lies in a measurable set $A$ is just the measure (length) of $A \cap [0,1]$, written $|A \cap [0,1]|$.
From this it is immediate that the cumulative distribution function (CDF) is
$$F_X(x) = \Pr(X \le x) = |(-\infty, x] \cap [0,1]| = |[0,\min(x,1)]| = \begin{array}{ll} \left\{ \begin{array}{ll} 0 & x\lt 0 \\ x & 0\leq x\leq 1 \\ 1 & x\gt 1. \end{array}\right. \end{array} $$
The probability density function (PDF), which is the derivative of the CDF, is $f_X(x) = 1$ for $0 \le x \le 1$ and $f_X(x)=0$ otherwise. (It is undefined at $0$ and $1$.)
The characteristic function (CF) of any random variable $X$ is the expectation of $\exp(i t X)$ (where $i$ is the imaginary unit, $i^2=-1$). Using the PDF of a uniform distribution we can compute
$$\phi_X(t) = \int_{-\infty}^\infty \exp(i t x) f_X(x) dx = \int_0^1 \exp(i t x) dx = \left. \frac{\exp(itx)}{it} \right|_{x=0}^{x=1} = \frac{\exp(it)-1}{it}.$$
The CF is a (version of the) Fourier transform of the PDF, $\phi(t) = \hat{f}(t)$. The most basic theorems about Fourier transforms are:
The CF of a sum of independent variables $X+Y$ is the product of their CFs.
When the original PDF $f$ is continuous and $X$ is bounded, $f$ can be recovered from the CF $\phi$ by a closely related version of the Fourier transform,
$$f(x) = \check{\phi}(x) = \frac{1}{2\pi} \int_{-\infty}^\infty \exp(-i x t) \phi(t) dt.$$
When $f$ is differentiable, its derivative can be computed under the integral sign:
$$f'(x) = \frac{d}{dx} \frac{1}{2\pi} \int_{-\infty}^\infty \exp(-i x t) \phi(t) dt = \frac{-i}{2\pi} \int_{-\infty}^\infty t \exp(-i x t) \phi(t) dt.$$
For this to be well-defined, the last integral must converge absolutely; that is,
$$\int_{-\infty}^\infty |t \exp(-i x t) \phi(t)| dt = \int_{-\infty}^\infty |t| |\phi(t)| dt$$
must converge to a finite value. Conversely, when it does converge, the derivative exists everywhere by virtue of these inversion formulas.
It is now clear exactly how differentiable the PDF for a sum of $n$ uniform variables is: from the first bullet, the CF of the sum of iid variables is the CF of one of them raised to the $n^\text{th}$ power, here equal to $(\exp(i t) - 1)^n / (i t)^n$. The numerator is bounded (it consists of sine waves) while the denominator is $O(t^{n})$. We can multiply such an integrand by $t^{s}$ and it will still converge absolutely when $s \lt n-1$ and converge conditionally when $s = n-1$. Thus, repeated application of the third bullet shows that the PDF for the sum of $n$ uniform variates will be continuously $n-2$ times differentiable and, in most places, it will be $n-1$ times differentiable.
The blue shaded curve is a log-log plot of the absolute value of the real part of the CF of the sum of $n=10$ iid uniform variates. The dashed red line is an asymptote; its slope is $-10$, showing that the PDF is $10 - 2 = 8$ times differentiable. For reference, the gray curve plots the real part of the CF for a similarly shaped Gaussian function (a normal PDF).
Let $Y$ and $X$ be independent random variables where $X$ has a Uniform $[0,1]$ distribution. Consider a narrow interval $(t, t+dt]$. We decompose the chance that $X+Y \in (t, t+dt]$ into the chance that $Y$ is sufficiently close to this interval times the chance that $X$ is just the right size to place $X+Y$ in this interval, given that $Y$ is close enough:
$$\eqalign{ f_{X+Y}(t) dt = &\Pr(X+Y\in (t,t+dt])\\ & = \Pr(X+Y\in (t,t+dt] | Y \in (t-1, t+dt]) \Pr(Y \in (t-1, t+dt]) \\ & = \Pr(X \in (t-Y, t-Y+dt] | Y \in (t-1, t+dt]) \left(F_Y(t+dt) - F_Y(t-1)\right) \\ & = 1 dt \left(F_Y(t+dt) - F_Y(t-1)\right). }$$
The final equality comes from the expression for the PDF of $X$. Dividing both sides by $dt$ and taking the limit as $dt\to 0$ gives
$$f_{X+Y}(t) = F_Y(t) - F_Y(t-1).$$
In other words, adding a Uniform $[0,1]$ variable $X$ to any variable $Y$ changes the pdf $f_Y$ into a differenced CDF $F_Y(t) - F_Y(t-1)$. Because the PDF is the derivative of the CDF, this implies that each time we add an independent uniform variable to $Y$, the resulting PDF is one time more differentiable than before.
Let's apply this insight, starting with a uniform variable $Y$. The original PDF is not differentiable at $0$ or $1$: it is discontinuous there. The PDF of $Y+X$ is not differentiable at $0$, $1$, or $2$, but it must be continuous at those points, because it is the difference of integrals of the PDF of $Y$. Add another independent uniform variable $X_2$: the PDF of $Y+X+X_2$ is differentiable at $0$,$1$,$2$, and $3$--but it does not necessarily have second derivatives at those points. And so on.
The CDF at $t$ of a sum of $n$ iid uniform variates equals the volume of the unit hypercube $[0,1]^n$ lying within the half-space $x_1+x_2+\cdots+x_n \le t$. The situation for $n=3$ variates is shown here, with $t$ set at $1/2$, $3/2$, and then $5/2$.
As $t$ progresses from $0$ through $n$, the hyperplane $H_n(t): x_1+x_2+\cdots+x_n=t$ crosses vertices at $t=0$, $t=1, \ldots, t=n$. At each time the shape of the cross section changes: in the figure it first is a triangle (a $2$-simplex), then a hexagon, then a triangle again. Why doesn't the PDF have sharp bends at these values of $t$?
To understand this, first consider small values of $t$. Here, the hyperplane $H_n(t)$ cuts off an $n-1$-simplex. All $n-1$ dimensions of the simplex are directly proportional to $t$, whence its "area" is proportional to $t^{n-1}$. Some notation for this will come in handy later. Let $\theta$ be the "unit step function,"
$$\theta(x) = \begin{array}{ll} \left\{ \begin{array}{ll} 0 & x \lt 0 \\ 1 & x\ge 0. \end{array}\right. \end{array} $$
If it were not for the presence of the other corners of the hypercube, this scaling would continue indefinitely. A plot of the area of the $n-1$-simplex would look like the solid blue curve below: it is zero at negative values and equals $t^{n-1}/(n-1)!$ at the positive one, conveniently written $\theta(t) t^{n-1}/(n-1)!$. It has a "kink" of order $n-2$ at the origin, in the sense that all derivatives through order $n-3$ exist and are continuous, but that left and right derivatives of order $n-2$ exist but do not agree at the origin.
(The other curves shown in this figure are $-3\theta(t-1) (t-1)^{2}/2!$ (red), $3\theta(t-2) (t-2)^{2}/2!$ (gold), and $-\theta(t-3) (t-3)^{2}/2!$ (black). Their roles in the case $n=3$ are discussed further below.)
To understand what happens when $t$ crosses $1$, let's examine in detail the case $n=2$, where all the geometry happens in a plane. We may view the unit "cube" (now just a square) as a linear combination of quadrants, as shown here:
The first quadrant appears in the lower left panel, in gray. The value of $t$ is $1.5$, determining the diagonal line shown in all five panels. The CDF equals the yellow area shown at right. This yellow area is comprised of:
The triangular gray area in the lower left panel,
minus the triangular green area in the upper left panel,
minus the triangular red area in the low middle panel,
plus any blue area in the upper middle panel (but there isn't any such area, nor will there be until $t$ exceeds $2$).
Every one of these $2^n=4$ areas is the area of a triangle. The first one scales like $t^n=t^2$, the next two are zero for $t\lt 1$ and otherwise scale like $(t-1)^n = (t-1)^2$, and the last is zero for $t\lt 2$ and otherwise scales like $(t-2)^n$. This geometric analysis has established that the CDF is proportional to $\theta(t)t^2 - \theta(t-1)(t-1)^2 - \theta(t-1)(t-1)^2 + \theta(t-2)(t-2)^2$ = $\theta(t)t^2 - 2 \theta(t-1)(t-1)^2 + \theta(t-2)(t-2)^2$; equivalently, the PDF is proportional to the sum of the three functions $\theta(t)t$, $-2\theta(t-1)(t-1)$, and $\theta(t-2)(t-2)$ (each of them scaling linearly when $n=2$). The left panel of this figure shows their graphs: evidently, they are all versions of the original graph $\theta(t)t$, but (a) shifted by $0$, $1$, and $2$ units to the right and (b) rescaled by $1$, $-2$, and $1$, respectively.
The right panel shows the sum of these graphs (the solid black curve, normalized to have unit area: this is precisely the angular-looking PDF shown in the original question.
Now we can understand the nature of the "kinks" in the PDF of any sum of iid uniform variables. They are all exactly like the "kink" that occurs at $0$ in the function $\theta(t)t^{n-1}$, possibly rescaled, and shifted to the integers $1,2,\ldots, n$ corresponding to where the hyperplane $H_n(t)$ crosses the vertices of the hypercube. For $n=2$, this is a visible change in direction: the right derivative of $\theta(t)t$ at $0$ is $0$ while its left derivative is $1$. For $n=3$, this is a continuous change in direction, but a sudden (discontinuous) change in second derivative. For general $n$, there will be continuous derivatives through order $n-2$ but a discontinuity in the $n-1^\text{st}$ derivative.
The integration to compute the CF, the form of the conditional probability in the probabilistic analysis, and the synthesis of a hypercube as a linear combination of quadrants all suggest returning to the original uniform distribution and re-expressing it as a linear combination of simpler things. Indeed, its PDF can be written
$$f_X(x) = \theta(x) - \theta(x-1).$$
Let us introduce the shift operator $\Delta$: it acts on any function $f$ by shifting its graph one unit to the right:
$$(\Delta f)(x) = f(x-1).$$
Formally, then, for the PDF of a uniform variable $X$ we may write
$$f_X = (1 - \Delta)\theta.$$
The PDF of a sum of $n$ iid uniforms is the convolution of $f_X$ with itself $n$ times. This follows from the definition of a sum of random variables: the convolution of two functions $f$ and $g$ is the function
$$(f \star g)(x) = \int_{-\infty}^{\infty} f(x-y)g(y) dy.$$
It is easy to verify that convolution commutes with $\Delta$. Just change the variable of integration from $y$ to $y+1$:
$$\eqalign{ (f \star (\Delta g)) &= \int_{-\infty}^{\infty} f(x-y)(\Delta g)(y) dy \\ &= \int_{-\infty}^{\infty} f(x-y)g(y-1) dy \\ &= \int_{-\infty}^{\infty} f((x-1)-y)g(y) dy \\ &= (\Delta (f \star g))(x). }$$
For the PDF of the sum of $n$ iid uniforms, we may now proceed algebraically to write
$$f = f_X^{\star n} = ((1 - \Delta)\theta)^{\star n} = (1-\Delta)^n \theta^{\star n}$$
(where the $\star n$ "power" denotes repeated convolution, not pointwise multiplication!). Now $\theta^{\star n}$ is a direct, elementary integration, giving
$$\theta^{\star n}(x) = \theta(x) \frac{x^{n-1}}{{n-1}!}.$$
The rest is algebra, because the Binomial Theorem applies (as it does in any commutative algebra over the reals):
$$f = (1-\Delta)^n \theta^{\star n} = \sum_{i=0}^{n} (-1)^i \binom{n}{i} \Delta^i \theta^{\star n}.$$
Because $\Delta^i$ merely shifts its argument by $i$, this exhibits the PDF $f$ as a linear combination of shifted versions of $\theta(x) x^{n-1}$, exactly as we deduced geometrically:
$$f(x) = \frac{1}{(n-1)!}\sum_{i=0}^{n} (-1)^i \binom{n}{i} (x-i)^{n-1}\theta(x-i).$$
(John Cook quotes this formula later in his blog post, using the notation $(x-i)^{n-1}_+$ for $(x-i)^{n-1}\theta(x-i)$.)
Accordingly, because $x^{n-1}$ is a smooth function everywhere, any singular behavior of the PDF will occur only at places where $\theta(x)$ is singular (obviously just $0$) and at those places shifted to the right by $1, 2, \ldots, n$. The nature of that singular behavior--the degree of smoothness--will therefore be the same at all $n+1$ locations.
Illustrating this is the picture for $n=8$, showing (in the left panel) the individual terms in the sum and (in the right panel) the partial sums, culminating in the sum itself (solid black curve):
It is useful to note that this last approach has finally yielded a compact, practical expression for computing the PDF of a sum of $n$ iid uniform variables. (A formula for the CDF is similarly obtained.)
The Central Limit Theorem has little to say here. After all, a sum of iid Binomial variables converges to a Normal distribution, but that sum is always discrete: it never even has a PDF at all! We should not hope for any intuition about "kinks" or other measures of differentiability of a PDF to come from the CLT.
I think the more surprising thing is that you get the sharp peak for $n=2$.
The Central Limit Theorem says that for large enough sample sizes the distribution of the mean (and the sum is just the mean times $n$, a fixed constant for each graph) will be approximately normal. It turns out that the uniform distribution is really well behaved with respect to the CLT (symmetric, no heavy tails (well not much of any tails), no possibility of outliers), so for the uniform the sample size needed to be "large enough" is not very big (around 5 or 6 for a good approximation), you are already seeing the OK approximation at $n=3$.
You could argue that the probability density function of a uniform random variable is finite,
so its integral the cumulative density function of a uniform random variable is continuous,
so the probability density function of the sum of two uniform random variables is continuous,
so its integral the cumulative density function of the sum of two uniform random variables is smooth (continuously differentiable),
so the probability density function of the sum of three uniform random variables is smooth.
