Points in a triangular lattice at the same distance from the origin and "breaking of symmetry" Introduction
I was trying to simulate what would happen to a certain physical system taking place in a triangular lattice (the physical details are not relevant to the discussion), when I came across a question that I could not answer and I am not familiar enough with the underlying maths to even know what to look for online.
Consider the ($2\mathrm{D}$) triangular lattice spaned by the two units vectors $e_1 = (1, 0)$ and $e_2 = (\frac{1}{2}, \frac{\sqrt{3}}{2})$. 
How I came across my problem (skip if only interested in the actual maths)
I was interested in a two-variable function $f(\overrightarrow{x_1}, \overrightarrow{x_2})$, with $\overrightarrow{x_1}$ and $\overrightarrow{x_2}$ belonging to the lattice. The function is translationally and rotationally invariant so that it only depends on the norm of the difference $f(\overrightarrow{x_1}, \overrightarrow{x_2}) = g(| \overrightarrow{x_1}- \overrightarrow{x_2}|)$. For simplicity, I will assume in the following that $\overrightarrow{x_1} = 0$ and note $\tilde{f}(\overrightarrow{x}) = f(\overrightarrow{0}, \overrightarrow{x}),$ with $\overrightarrow{x}$ in the lattice
My (poorly-written) code was giving the expected results for $\overrightarrow{x} = e_1, \overrightarrow{x} = e_2, \overrightarrow{x} = e_1 + e_2, \overrightarrow{x} = 2 e_1, \dots$ In fact, if you regroup together all vectors $\overrightarrow{x}$ within a same distance of the origin, and you sort the groups ($X_1, X_2, \dots, X_n, \dots)$ from the smallest to the largest distance (so $X_1$ would be the six sites nearest to the origin, $X_2$ the six next-nearest sites, and so on and so forth), I had the right results for all $X_1$, $X_2$, $X_3$ all the way up to $X_{19}$. But for some reason it was breaking down at $X_{20}$. I first thought that this was because my code was not good for large distances, but it worked again for $X_{21}$, $X_{22}, \dots$ so this was really strange
To help visualize this I plotted the different $X_i$'s for $i$ from $1$ to $30$ as seen below (sorry for the horrible color choice but it is to help distinguish groups that are close to each other):

The troublesome sites are the ones that are circled in red. After a while trying to figure out what was wrong, I realized that my code (I told you it was poorly written) is probably making the implicit assumption that all sites within a group were "equivalent" in some sense (any of them can be transformed into any other from the same group by a rotation of angle $k \times 60^{\circ}$, $k \in \{1, 2, 3, 4, 5 \}$ around the origin, and/or a reflexion with respect to the X or Y axis). But it is easy to see that this is not the case for $X_{20}$, as this group contains $18$ vectors instead of $6$ or $12$. This comes from the fact that $3 e_1 + 5 e_2$ and $7 e_2$ have the same norm while not being related by a rotation of $k \times 60^{\circ}$ and/or reflexion. I do not know the next time it happens, but from looking at the first $30$ groups of vectors, this seems to be the only occurence.
Actual maths
I guess this is a well-known problem in group theory/number theory, but can you find all points $\overrightarrow{x_1}$, $\overrightarrow{x_2}$ in a triangular lattice that are at the same distance from the origin while breaking the symmetries of the original lattice? How does this expand to different type of lattices and/or different dimensions? The closest thing that I am familiar with is the equivalent problem for a square lattice, which is to find all non-trivial integer solutions to the equation:
$$ a^2 + b^2 = c^2 + d^2, \quad (a,b,c,d) \in \mathbb{N}^4$$
For instance, $(6, 8, 10, 0)$ is an acceptable solution. But I am looking here for solutions in the case of a triangular lattice, and possibly how to generalize it to different lattice types.
 A: This seems to be the prettiest one: for
$$  u^2 + uv + v^2 = x^2 + xy+y^2,  $$
take coprime integers $p,q,r,s$  and
$$  u = rs+ps- rq  $$
$$  v = pq-ps + rq  $$
$$  x = pq -rs + rq  $$
$$  y = -pq+rs+ps  $$
$$
\left(
\begin{array}{r}
      u  \\
      v    \\
      x    \\
      y    \\
\end{array}
\right) =
\left(
\begin{array}{rrrr}
      0   &   1   &   1 &    -1 \\
      1   &   0  &   -1 &     1 \\
      1   &  -1  &    0  &    1 \\
     -1   &   1  &    1   &   0 \\
\end{array}
\right)
\left(
\begin{array}{r}
      pq   \\
      rs    \\
      ps    \\
      rq    \\
\end{array}
\right)
$$

parisize = 4000000, primelimit = 500000
?   u = 0*p*q+r*s+p*s- r*q ;   v = p*q+ 0*r*s-p*s +r*q ;     x = p*q -r*s +0*p*s+ r*q ;    y = -p*q+r*s+p*s + 0* r*q  ;
? u
%2 = s*p + (-r*q + s*r)
? v
%3 = (q - s)*p + r*q
? x
%4 = q*p + (r*q - s*r)
? y
%5 = (-q + s)*p + s*r
? u^2 + u*v + v^2
%6 = (q^2 - s*q + s^2)*p^2 + (r*q^2 - s*r*q + s^2*r)*p + (r^2*q^2 - s*r^2*q + s^2*r^2)
? f = u^2 + u*v + v^2
%7 = (q^2 - s*q + s^2)*p^2 + (r*q^2 - s*r*q + s^2*r)*p + (r^2*q^2 - s*r^2*q + s^2*r^2)
? g = x^2 + x*y + y^2
%8 = (q^2 - s*q + s^2)*p^2 + (r*q^2 - s*r*q + s^2*r)*p + (r^2*q^2 - s*r^2*q + s^2*r^2)
? f - g
%9 = 0
? 
? 


A: From Dickson, 1929:
The number of representations of any positive $n$ by $q=x^2 +xy+y^2$ is $6E(n),$ where $E(n)$ is the excess of the number of divisors $3h+1$ over the number of divisors $3h+2.$ If $n=2^k m,$ then we have $E(n) = 0$ when $k$ is odd, but when $k$ is even we have $E(n) = E(m).$
I should add that the situation with the prime $2$ is the same as the situation with prime $5$ or $11$ or any $p \equiv 2 \pmod 3.$  If, in factoring the number $n,$ the exponent of such a prime $p$ is odd, there are no representations at all. If the exponent of such $p$ is even, dividing by $p^2$ does not change the number of representations, so we may keep dividing by such $p^2$ until the only prime factors are $3$ itself and any $q \equiv 1 \pmod 3. $ Multiplying by $3$ does not change the number of representations (needs a little proof) so we keep dividing by that as well. In the end, we are asking for the number of divisors of a number that is a product of primes $q_j \equiv 1 \pmod 3.$ Then multiply by $6$
