Why does a non-linear system lead to interaction and frequency mixing between inputs? When we have a system that is nonlinear and we apply a sum of two different frequency sine waves as an input, we see the output of this system has components that are at the sum frequency of the two sine waves, the difference frequency of the two sine waves, and multiple harmonics for each of the sine waves. Why is this? Why does a nonlinear system allow for interaction/mixing while a linear system does not?
 A: Really, the remarkable thing is that a linear system does not ever lead to any interaction between frequencies. You put two sine waves waves in, and you're guaranteed to get exactly two sine waves out again. Why would that be? Why is it not possible for anything else to come out?
Well, it is essentially the definition of linearity: a function $f:V\to W$ is linear if for any $u,v \in V$ and $\mu,\nu \in \mathbb{R}$,
$$
  f(\mu\cdot u + \nu\cdot v) = \mu\cdot f\,u + \nu\cdot f\,v.
$$
Specifically, $f(u + v) = f\,u + f\,v$, so if both $f\,u$ and $f\,v$ are just single sine waves than $f(u + v)$ will be simply the superposition of two sine waves, nothing else.
But in the general, i.e. nonlinear case, there is no reason whatsoever to expect that this should be the case. $f(u + v)$ could be something completely and utterly different from $f\,u + f\,v$.
In fact, even the notion of “frequency mixing” and so on again assumes that the function is of a very simple form, namely dominated by linear and quadratic terms. It is the quadratic contribution that gives you sum- and difference frequencies. The reason being
$$
  (\sin (\omega_0\cdot t) + \sin(\omega_1\cdot t))^2
    = \sin(\omega_0\cdot t)^2
        + 2\cdot \sin (\omega_0\cdot t) \cdot\sin(\omega_1\cdot t)
        + \sin(\omega_1\cdot t)^2
$$
The $\sin(\omega_i\cdot t)^2$ terms† are often ignored. The reason is that these can be written in terms of $\sin(2\cdot\omega_i\cdot t)$, and double-frequency is often already present in the signal anyway (real-world signals can very well be periodic at frequency $\tfrac{\omega_i}{2\pi}$, but they won't be exact sinusoidals, meaning they can be interpreted as a Fourier series of integer-multiple frequencies). But
$$
  2\cdot \sin (\omega_0\cdot t) \cdot\sin(\omega_1\cdot t)
   = \cos ((\omega_0-\omega_1)\cdot t) - \cos ((\omega_0+\omega_1)\cdot t)
$$
...and those sum-and difference frequencies were quite definitely not at all in either of the individual signals.
For a sufficiently pathological nonlinear function, the quadratic approximation will be no good either – you'll get a whole mess of frequencies across the spectrum out of only two sines. But often, what we're interested in (or try to build) are systems that are to good approximation linear, and then the remainder is small and can be modelled by a quadratic.

†I reject the commonly used, but inconsistent $\sin^2 x$ notation. $\sin(x)^2 = (\sin x)^2$, whereas $\sin^2x$ should actually mean $\sin(\sin x)$, at least if $\tan^{-1}$ denotes the arctangent.
