Patterns in solutions to $a^2 + b^2 + c^2 = n$ I have plotted solutions $(a,b,c)$ to $a^2 + b^2 + c^2 = n$ for $12000 \leq n < 12100$, rescaled to $S^2$ and projected onto the first two coordinates. (these are read from the lower left, across and upwards.  Sorry.)

While, over all $ a \leq  n \leq b$ there is clear tendency towards equidistribution, for a fixed $n$, there are marked patterns (or even no solution at all).  I've been struggling to find language for the types of patterns that are observing.  There are solution-free regions, there are linear patterns, circular patterns and other sub-varieties.  
I suspect for fixed $n$, the solutions are clustering along the intersection of two vareties, $S^2 \cap V $ and I am trying to characterize the equations of $V$.
This is my theory of why equidistribution might be so hard to prove; it's because there are in fact patterns. 
 A: There is no clustering of the solutions of $a^2+b^2+c^2=n$, even for individual $n$'s, assuming the number of solutions is large (e.g. when $n\equiv 1,2,3,5,6\pmod{8}$ and $n$ is large). This was proved by William Duke (Invent. Math. 92 (1988), 73-90), his paper is available (for free) here. 
For more recent results, e.g. what happens beyond equidistribution, see the work of Bourgain-Sarnak-Rudnick here and here.
A: There is a $48$-element symmetry group of the space of solutions. This by itself will create the appearance of the patterns when the number of solutions divided by $48$ is small. Imagine choosing $k$ random points in the fundamental domain $a \geq b \geq c \geq 0$ and reflecting them around the sphere. Patterns will appear, simply because it is unlikely for the $k$ points to be uniform in the triangle, and any nonuniformities will be magnified by some or all of the symmetries.
In your case $k \approx 2 \pi (12000)^{1/2} / 48 \approx 14$ is pretty small, and this might partially or entirely explain all the effect you see.
A: Here I plotted some solutions around $n= 120,000 \pm 50$ and hoping a parallel algorithm could accelerate the process for $n \approx 1.2 \times 10^6$.  Here we notice again, distinctive, but fainter patterns.

At around $n \asymp 10^6$ these solution sets look pretty homogeneous. I'm a little bit nervious about floating point artifacts or other programming issues.    The solution sets look like random spherical harmonics.  This might even have a name in the physics literature.  
 
There might be approximate algorithms for $n \asymp 10^8$ .
