$\left(a^2+ab+b^2\right)\left(c^2+cd+d^2\right)=\left(ac-bd\right)^2+\left(ac-bd\right)\left(ad+bc+bd\right)+\left(ad+bc+bd\right)^2$

I claim that all primes of the form $3k+1$ can be written in the form $a^2+ab+b^2$ for some positive integers $a$ and $b$. It is well known that there exists infinitely many primes of this form. Lemma: Given any prime of the form $3k+1$, there exists a solution to the equation $x^2+x+1=0\left(\operatorname{mod}p\right)$ Proof: This is equivalent to the quadratic congruence $\left(2x+1\right)^2=-3\left(\operatorname{mod}p\right)$. So, using Legrende notation, we need to show that $\left(\frac{-3}{p}\right) = 1$. We also have $\left(\frac{-3}{p}\right)=\left(\frac{3}{p}\right)\left(\frac{-1}{p}\right)$. It is known that $\left(\frac{-1}{p}\right)=\left(-1\right)^{\frac{p-1}{2}}$, and by quadratic reciprocity, $\left(\frac{3}{p}\right)\left(\frac{p}{3}\right)=\left(-1\right)^{\frac{p-1}{2}\cdot\frac{3-1}{2}}=\left(-1\right)^{\frac{p-1}{2}}$. But as $p=3k+1$, $\left(\frac{p}{3}\right)=1$ and thus $\left(\frac{3}{p}\right)=\left(-1\right)^{\frac{p-1}{2}}$ and so finally $\left(\frac{-3}{p}\right)=\left(-1\right)^{p-1}=1$ So by our lemma, for any prime of the form $3k+1$, there exists an $x$ such that $x^2+x+1=mp$ for some integer $m$. So we have integers $a, b$ such that $a^2+ab+b^2=mp$ for some integer $m$. Note that as $x\le p-1$, $m < p$ and thus is coprime to $p$. I will now show that either $m=1$, and so we are done, or we can always find a smaller set of solutions to $a^2+ab+b^2=mp$ for some integer $m$.Suppose $m$ is strictly greater than $1$ and let $a=a'\left(\operatorname{mod}m\right)$ and $b=b'\left(\operatorname{mod}m\right)$ where $\frac{-m}{2}\le a',\ b'\le\frac{m}{2}$ Now consider the following: $\frac{\left(a^2+ab+b^2\right)\left(b'^2+a'b'+a'^2\right)}{m^2}=\left(\frac{ab'-ba'}{m}\right)^2+\left(\frac{ab'-ba'}{m}\right)\left(\frac{aa'+bb'+ba'}{m}\right)+\left(\frac{aa'+bb'+ba'}{m}\right)^2$ From our definitions of $a'$ and $b'$, and from the knowledge $a^2+ab+b^2=0\left(\operatorname{mod}m\right)$, it is clear that our new solution does consist of integers. Also, the result is still a multiple of $p$, as $m$ and $p$ are coprime. Not only that, but $a'^2+a'b'+b'^2\le\frac{3}{4}m^2$ and thus if we were to write this new solution as $np$ for some integer $n$, it is clear than $n<m$. Hence we repeat this process until $m=1$ and hence we are done.