Does there exist any straight-forward theory that can explain what happens to colored noise that passes through a static non-linearity? That is, if you have colored noise $$v = H u \, , \; u \sim N(0,\sigma)$$ generated by driving an LTI filter $H$ by white noise, and pass it though a static non-linear function $y = f(v)$, e.g. $f(v) = v^2$; what is the spectrum of $y$?
If your nonlinearity can be expressed as a polynomial (i.e., in terms of addition and multiplication), you can make use of:
The linearity of the Fourier transform, i.e., if $f$ and $g$ are (benign) functions, $a$ and $b$ are numbers and $ℱ$ denotes the Fourier transform, then:
$$ℱ(a·f+b·g) = a·ℱ(f) + b·ℱ(g)$$
The convolution theorem, which states that multiplication in the time domain corresponds to convolution in the frequency domain (and vice versa):
$$ℱ(f·g) = ℱ(f) ∗ ℱ(g)$$
So, in your specific example, the spectrum of $y$ is the spectrum of $v$ convoluted with itself:
$$ℱ(y) = ℱ(v·v) = ℱ(v) ∗ ℱ(v)$$
It is impossible to derive an analytic expression that describes what happens to the spectrum.
From experience however, passing your signal through that $f$, it will increase the high frequencies. This is because your new mapping ($f$) is increasing the slope of the signal (dependant on amplitude). If that signal was to be fed to a piece of electronic equipment, then the waveform would initially distort and ultimately clip.
That is, it would hit the maximum voltage (or current) it is designed for and stay there. In short, sinusoids turn to square waves. This is the basis of electric guitar distortion.
The ...pass it though a static non-linear function... is called "Function Composition". In some cases, deriving the impact of some composition to the frequency spectrum might be "easy". That is, there might be analytical expressions that characterise the Fourier Transform's coefficients. But in the general case, what you are likely to end up with after composition, is a new function whose Fourier sum (or integral) you are going to have to evaluate, in order to "see what does it look like".
In terms of "...a straightforward theory..." there are basically two things that are relatively close to what you are looking for. The first thing is the Taylor series expansion and the second thing is nonlinear or dynamic convolution.
Taylor series expansion is used for the linearisation of a nonlinear function around a given operating point. It would basically tell you how much would your original function's slope change around that point (within a given range of course, preferably small, to be accurate). But if the signal "crossed" that band of observation, you would have to pick a new point and linearise the function again and so on. Again, for simple functions this might be "easy".
Related to this is Dynamic Convolution (and also, this paper for what you are after) which uses not one $h$, but one $h$ per amplitude range. Therefore, the resulting signal that is passing through a system with a static nonlinearity (because if it is not static, then you have one $h$ per amplitude per unit of time) is the sum of a set of conditional convolutions depending on the input waveform's amplitude.
These could be used to represent the effect of the non-linearity on the input waveform (or function) which might provide some insight for specific waveforms.
But it is impossible to derive an analytic expression of what is happening to the spectrum of some $u$ as it passes through some non-linearity like $f(x)=x^2$.
In fact, if you look at it from a practical point of view, non-linearities tend to create harmonics. That is, they tend to add stuff to the spectrum of $u$. Linear systems do not do that because a simple change in linear slope is basically a multiplication. Say for instance your new mapping is some $f(x,a)=a \cdot x$ (instead of your $f(x)=x^2$). Clearly, you change the $a$ you change the slope of $x$ but due to the linearity property of the Fourier Transform, the impact this has on the coefficients is to simply enlarge them but not move them around or generate new ones.
A sinusoid through a non-linearity at low amplitudes, may be "riding" its linear part and therefore, one spectral line enters the system, one spectral line exits the system. BUT!, as the sinusoid grows louder and approaches clipping, it starts to look more and more like a square wave and as a square wave it is composed of a sum of odd-integer harmonics. That is, one spectral line enters the system, many spectral lines exit the system.
It is this enrichment that it is impossible to describe analytically.
Hope this helps.