Why audio clipping produces harmonics?
Clipping is a non-linear operation. Consider a system that clips all input whose absolute value is larger than 1: $$ y(t) = \begin{cases}x(t)\text{ if }|x(t)|\leq 1 \\ 1 \text{ if } x(t) &gt; 1 \\ -1 \text{ if } x(t) &lt; -1 \end{cases}$$ While this system is linear for small inputs, it becomes non-linear as soon as the input becomes large enough.

Non-linear systems with no memory can be modeled by a power series: \begin{align}y(t) &amp;= \sum_{n=0}^\infty a_n x^n \\ &amp;= a_0 + a_1x(t) + a_2x^2(t) + a_3 x^3(t) + \cdots \end{align} Usually the coefficients $a_n$ become small as $n$ increases. Let's focus on $n=2$. If the input is a sine wave, then the clipped output contains the square of the sine wave, which has a DC term plus a term of twice the frequency of the input. It is the first distortion-causing harmonic.

If you play a bit with the equations you'll see that the terms for $n=3,4,\ldots$ each introduces new harmonics. When the input is more complex, such as a sum of sine waves of different frequencies, the harmonics will be more interesting.

Another, maybe more intuitive way to see it is as follows. A clipped sine wave looks a bit like a square wave -- the more clipping there is, the more similar it is. It is natural, then, that its Fourier series will resemble that of a square wave, which has plenty of harmonics.

If you have access to a spectrum analyzer, you can easily experiment with this by producing clipped sound in your computer, and seeing the output in the spectrum analyzer.