The set of integers $n$ expressible as $n=x^2+xy+y^2$ 
Let $S$ be the set of integers $n$, such there exist integers $x,y$ with $$n=x^2+xy+y^2$$


Is the implication $$a,b\in S\implies ab\in S$$ true? If yes, how can I prove it?

I worked out $$n\in S\iff 4n\in S$$ and $$n\in S\iff 3n\in S$$
I tried two approaches. The first is to express $$(a^2+ab+b^2)(c^2+cd+d^2)$$ in the form $$f^2+fg+g^2$$ with polynomials $f,g$ with integer coefficients. I however could not find suitable $f$ and $g$.
The second approach is based on $$x^2+xy+y^2=\frac{(2x+y)^2+3y^2}{4}$$ If we have $n=x^2+xy+y^2$ , we have $u^2+3v^2=4n$ for some integers $u,v$ with equal parity. The main problem of this approach is to consider the equal parity.
Any ideas ?
 A: Denote $a^2+ab+b^2$ by $L(a, b)$. Numbers of the form $L$ are known as Loeschian numbers. What you are trying to prove is that the set $S$ of Loeschian numbers is multiplicative.
Here's a direct proof that the product of two Loeschian numbers is Loeschian.
\begin{align*}
L(a, b)L(c, d) &=(a^2+ab+b^2)(c^2+cd+d^2)\\
&= a^2c^2+a^2cd+a^2d^2+abc^2+abcd+abd^2+b^2c^2+b^2cd+b^2d^2\\
&=[a^2c^2+2abc^2+b^2c^2+2abcd+2b^2cd+b^2d^2]+{}\\
&\qquad +[a^2cd-abc^2+abcd-b^2c^2+abd^2-b^2cd]+{}\\
&\qquad +[a^2d^2-2abcd+b^2c^2]\\
&=(ac+bc+bd)^2+(ac+bc+bd)(ad-bc)+(ad-bc)^2\\
&=L(ac+bc+bd, ad-bc)
\end{align*}
In addition,
\begin{align*}
L(a, b)L(c, d)
&= (a^2+ab+b^2)(c^2+cd+d^2)\\
&= a^2c^2+a^2cd+a^2d^2+abc^2+abcd+abd^2+b^2c^2+b^2cd+b^2d^2\\
&=[a^2c^2-2abcd+b^2d^2]+{}\\
&\qquad +[a^2cd+abc^2+abcd-abd^2-b^2cd-b^2d^2]+{}\\
&\qquad +[a^2d^2+2abcd+b^2c^2+2abd^2+2b^2cd+b^2d^2]\\
&=(ac-bd)^2+(ac-bd)(ad+bc+bd)+(ad+bc+bd)^2\\
&=L(ac-bd, ad+bc+bd)
\end{align*}
It can also be proved that $n$ has the form $a^2+ab+b^2$ iff $n$ has the form $u^2+3v^2$. (Sketch of a proof: Where $a+b+c=0$, $L(a, b)=L(b, c)=L(c, a)$. Moreover, at least one of $a, b$ and $c$ is even. $L(a, 2v)=a^2+2av+4v^2=u^2+3v^2$ where $u=a+v$. Conversely, $u^2+3v^2=L(a, b)$ where $b=2v$ and $a=u-v$. Thus $L$ and $u^2+3v^2$ represent the same set of integers.)
Thus the set of numbers of the form $u^2+3v^2$ is multiplicative. Again there is a direct proof, and it is slightly simpler than the proof above.
A: The answer is yes.  If $m=u^2+uv+v^2$ and $n=x^2+xy+y^2$, then
$$m=(u+\omega v)(u+\bar{\omega}v)\text{ and }n=(x+\omega y)(x+\bar{\omega}y)\,,$$
where $\omega:=\frac{1+\sqrt{-3}}{2}$ and $\bar{\omega}:=\frac{1-\sqrt{-3}}{2}$.  Now,
$$\begin{align}(u+\omega v)(x+\omega y)&=ux+\omega (uy+vx)+\omega^2 vy=ux+\omega(uy+vx)+(\omega -1)vy
\\
&=(ux-vy)+\omega (uy+vx+vy)\,.\end{align}$$
since $\omega^2-\omega +1=0$.  Thus,
$$mn=(f+\omega g)(f+\bar{\omega}g)=f^2+fg+g^2\,,$$
where $f:=ux-vy$ and $g:=uy+vx+vy$.

In fact, $S$ consists of all natural numbers of the form $$n:=3^\alpha \prod_{i=1}^r\,p_i^{\beta_i}\,\prod_{j=1}^s\,q_j^{2\gamma_j}\,.$$ where $\alpha,r,s,\beta_1,\beta_2,\ldots,\beta_r,\gamma_1,\gamma_2,\ldots,\gamma_s\in\mathbb{Z}_{\geq 0}$, $p_1<p_2<\ldots<p_r$ are prime natural numbers congruent to $1$ modulo $3$, and $q_1<q_2<\ldots<q_s$ are prime natural numbers congruent to $2$ modulo $3$.  Note that, for such $n$, there are precisely $N_n$ pairs $(x,y)\in\mathbb{Z}\times \mathbb{Z}$ such that $n=x^2+xy+y^2$, where
$$N_n:=6\,\prod_{i=1}^r\,\left(\beta_i+1\right)\,.$$
For example, $N_3=6$ as $$(x,y)=\pm(1,1), \pm (2,-1), \pm(-1,2)$$ are all the integral solutions to $3=x^2+xy+y^2$ . Another example is $N_7=12$, since $$(x,y)=\pm(2,1),\pm(1,2),\pm(3,-1),\pm(-1,3),\pm(3,-2),\pm(-2,3)$$
are the integral solutions to $7=x^2+xy+y^2$.  Also, there are $$\left\lceil\frac{N_n}{12}\right\rceil=\left\lceil\frac{1}{2}\prod_{i=1}^r\,\left(\beta_i+1\right)\right\rceil$$
pairs $(x,y)\in\mathbb{Z}_{\geq 0}\times\mathbb{Z}_{\geq 0}$ with $x\geq y$ such that $n=x^2+xy+y^2$.
A: Note that
\begin{align*}
&\qquad16(a^2+ab+b^2)(c^2+cd+d^2)\\&=\left((2a+b)^2+3b^2\right)\left((2c+d)^2+3d^2\right)\\&=\left(\frac{(2a+b)(2c+d)+3bd}{2}\right)^2+3\left(\frac{d(2a+b)-b(2c+d)}{2}\right)^2
\end{align*}
