Is the representation of any prime of the form $6n+1$ as $a^2+3b^2$ essentially unique? Two well-known results in number theory are:

Fermat's $4n+1$ theorem: Every prime of the form $4n+1$ can be represented as $a^2+b^2 (a,b \in\mathbb{N})$.


Euler's $6n+1$ theorem: Every prime of the form $6n+1$ can be represented as $a^2+3b^2 (a,b \in\mathbb{N})$.

Looking at the Mathworld entries on these theorems here and here, I notice that representation of primes of the form $4n+1$ is stated to be unique (up to order), but that there is no mention of uniqueness in respect of representation of primes of the form $6n+1$.  Uniqueness does however seem to hold at least for small primes of this form.
Question: Is the representation of any prime of the form $6n+1$ as $a^2+3b^2$ essentially unique?
If this is the case then a reference to a proof would be appreciated.
 A: This follows from very old results on representations of integers by quadratic forms. In particular it is a special case of a result of Euler which shows that two essentially distinct representations of $\,m\,$ imply  $\,m\,$ is composite $ $ (the proof constructs a proper factor of $\,m\,$ via a quick gcd computation). $\, $
Appended below is a classic elementary proof of Euler's result that requires no knowledge of ideal theory of quadratic number fields. It is excerpted from Wieb Bosma's thesis (1990) pp. 14-16 (which has a nice concise historical introduction to primality testing). It deserves strong emphasis that the arithmetical essence of this proof is much clearer when it is translated into the language of quadratic number fields and their ideal theory - as is often the case for such results. The use of ideals essentially simplifies such (nonlinear) quadratic (form) arithmetic by linearizing it into arithmetic of ideals (modules) - making available analogs of the powerful tools of linear algebra.
As hinted in the final few paragraphs below, this result was part of Euler's research on idoneal numbers for primality testing. For more on such see Ernst Kani's paper.




[10] Z.I. Borevich, I. R. Shafarevich, Number Theory, Orlando: Academic Press 1966.
[159] A. Weil, Number theory, an approach through history, Boston: Birkhauser 1984.
A: This problem has a completely elementary proof, using just high-school calculus and basic properties of divisibility.

Proposition: Let $n \geq 2$ be an arbitrary natural number, and let $p$ be a prime number such that $p \nmid n$. If the prime number $p$ has the form $x^2+ny^2$, for some natural numbers $a$ and $b$, then this representation is necessarily unique. (We can eliminate the condition $p \nmid n$.) (In the case of your problem you can replace it with $n=3$.)


In other words, for $n \geq 2$, a prime number either has no representation of the form $x^2+ny^2$, or if it does, this representation is unique, up to sign. (In the case of $n=1$, a prime number either has no representation of the form $x^2+y^2$, or if it does, this representation is unique, up to sign and permutations.)


Proof: Suppose on contrary that for positive integers $(a, b)\neq (c, d) \in \mathbb{N}^2$ we have:
$$a^2+nb^2=p=c^2+nd^2.$$
Then clearly we have
$$a^2\cong -nb^2 \mod p,$$
$$c^2\cong -nd^2 \mod p.$$
Multiplying respectively and vice versa together, we get that:
$$(ac)^2\cong (nbd)^2 \mod p,$$
$$-n(ad)^2\cong -n(bc)^2 \mod p.$$
Since $p \nmid n$ we get that:
$$p \mid (ac-nbd)(ac+nbd),$$
$$p \mid (ad+bc)(ad-bc).$$

Let $\alpha_+=ac-nbd$, $\alpha_-=ac+nbd$, $\beta_+=ad+bc$, $\beta_-=ad-bc$. What we get is translated in this way: At least one of $\alpha_{\pm}$ is divisible by $p$ (Some information in the congruences are extra, but I think it is better not to remove them to see the symetries.)
Now look at these two relations(Brahmagupta's identity):
$$p^2=p.p=(a^2+nb^2)(c^2+nd^2)=(ac-nbd)^2+n(ad+bc)^2=\alpha_+^2+n\beta_+^2,$$
$$p^2=p.p=(a^2+nb^2)(c^2+nd^2)=(ac+nbd)^2+n(ad-bc)^2=\alpha_-^2+n\beta_-^2.$$
We know that $p$ divides at least one of $\alpha_{\pm}$.

*

*Claim: If $p\mid \alpha_+$, then $p \mid \beta_+$.
Proof: If $p\mid \alpha_+$, then we have $p\mid \alpha_+^2$. Also we know that $p\mid p^2$, so we can conclude that $p \mid (p^2-\alpha_+^2)$. By the first identity ($p^2=\alpha_+^2+n\beta_+^2,$), we can conlude that $p\mid n\beta_+^2$. Since $p\nmid n$, by Euclid's lemma we can conclude that $p \mid \beta_+^2$. Finaly notice that $p$ is a prime, so we can conclude that $p \mid \beta_+$.


*Claim: If $p\mid \alpha_-$, then $p \mid \beta_-$.
Proof: If $p\mid \alpha_-$, then we have $p\mid \alpha_-^2$. Also we know that $p\mid p^2$, so we can conclude that $p \mid (p^2-\alpha_-^2)$. By the second identity ($p^2=\alpha_-^2+n\beta_-^2,$), we can conlude that $p\mid n\beta_-^2$. Since $p\nmid n$, by Euclid's lemma we can conclude that $p \mid \beta_-^2$. Finaly notice that $p$ is a prime, so we can conclude that $p \mid \beta_-$.
On the other hand we know that $p\mid \alpha_{\pm}$. If we consider these two claims, then we can conclud that one of the pairs $(\dfrac{\alpha_+}{p},\dfrac{\beta_+}{p})$ or $(\dfrac{\alpha_-}{p},\dfrac{\beta_-}{p})$ is a pair of integers.

Case (I): Suppose that $(\dfrac{\alpha_+}{p},\dfrac{\beta_+}{p}) \in \mathbb{Z^2}$, then form the first Brahmagupta's identity we get that:
$$1=\dfrac{p^2}{p^2}=(\dfrac{\alpha_+}{p})^2+n(\dfrac{\beta_+}{p})^2,$$
but notice that this equation has just only these solutions: $(\dfrac{\alpha_+}{p},\dfrac{\beta_+}{p})=(\pm 1, 0)$, which is equivalent to $(ac-nbd, ad+bc)=(\pm p, 0)$. Replace $c=-a\dfrac{d}{b}$ in the relation $ac-nbd= \pm p$, which gives you:
$$\pm p =ac-nbd=-a.a\dfrac{d}{b}-nbd\dfrac{b}{b}=-\dfrac{d}{b}(a^2+nb^2)=-\dfrac{d}{b}p.$$
Now look at these three ralation: $(a,b)\neq(c,d) \in \mathbb{N}^2$, $ad+bc=0$, and $\dfrac{d}{b}=\pm 1$. It is impossible. (Notice that the relation $ad+bc=0$ is equivalent to $\dfrac{c}{a}=-\dfrac{d}{b}$, so we can conclude that $a=\pm c$ and $b=\pm d$, ... to get a contradiction.)
(In this case whence we realize that $ad+bc=0$, we were able to get a contradiction, but for the sake of symmetry I write it in that form.)

Case (II): Suppose that $(\dfrac{\alpha_-}{p},\dfrac{\beta_-}{p}) \in \mathbb{Z^2}$, then form the second Brahmagupta's identity we get that:
$$1=\dfrac{p^2}{p^2}=(\dfrac{\alpha_-}{p})^2+n(\dfrac{\beta_-}{p})^2,$$
but notice that this equation has just only these solutions: $(\dfrac{\alpha_-}{p},\dfrac{\beta_-}{p})=(\pm 1, 0)$, which is equivalent to $(ac+nbd, ad-bc)=(\pm p, 0)$. Replace $c=a\dfrac{d}{b}$ in the relation $ac+nbd=\pm p$, which gives you:
$$\pm p =ac+nbd=a.a\dfrac{d}{b}+nbd\dfrac{b}{b}=\dfrac{d}{b}(a^2+nb^2)=\dfrac{d}{b}p.$$
Now look at these three ralation: $(a,b)\neq(c,d)$, $ad-bc=0$, and $\dfrac{d}{b}=\pm 1$. It is impossible. (Notice that the relation $ad-bc=0$ is equivalent to $\dfrac{c}{a}=\dfrac{d}{b}$, so we can conclude that $a=\pm c$ and $b=\pm d$, ... to get a contradiction.)

Remark: The condition $p\nmid n$ is not a serious condition. Because if $p\mid n$ and $p=a^2+nb^2$, then necessarily we have $a=0$, and by simple argument we can show that $b=1, n=p$, and the representation is unique clearly.
The proof for the fact about $n=1$, is similar.
A: If you require $a, b \in \mathbb{N}$ then the representation is literally unique, and we can argue as Daniel Fischer does in the comments: the ring of Eisenstein integers $\mathbb{Z}[\omega]$, where $\omega = \frac{-1 + \sqrt{-3}}{2}$ is a primitive third root of unity, is a unique factorization domain, and the primes congruent to $1 \bmod 6$ admit a further factorization
$$p = (x - y \omega)(x - y \omega^2) = x^2 + xy + y^2$$
where $x - y \omega$ is a prime in $\mathbb{Z}[\omega]$ and $x - y \omega^2$ is its conjugate. (This can be proven by inspecting the quotient $\mathbb{Z}[\omega]/p \cong \mathbb{F}_p[\omega]/(\omega^2 + \omega + 1)$.) We can write $x - y \omega = x + \frac{y}{2} - \frac{y}{2} \sqrt{-3}$  which gives
$$p = \left( x + \frac{y}{2} \right)^2 + 3\left( \frac{y}{2} \right)^2$$
and this representation consists of integers as long as $y$ is even. Now, $x - y \omega$ in the factorization above is unique up to multiplication by units, and the units of $\mathbb{Z}[\omega]$ are $\pm 1, \pm \omega, \pm \omega^2$. Multiplying by $\pm 1$ doesn't affect the values of $|x|$ and $|y|$ so now we inspect only the results of multiplying by $\omega$ and $\omega^2$. This gives
$$(x - y \omega) \omega = x \omega - y \omega^2 = x \omega + y (\omega + 1) = y + (x + y) \omega$$
$$(x - y \omega) \omega^2 = x \omega^2 - y = - x (\omega + 1) - y = (- x - y) - x \omega.$$
This means that if $p = a^2 + 3b^2 = (a + \sqrt{-3} b)(a - \sqrt{-3} b)$ then $b$ can only take on one of the six values $\pm \frac{y}{2}, \pm \frac{x+y}{2}, \pm \frac{x}{2}$, and some casework gives that depending on the parities of $x$ and $y$, exactly one of $x, y, x + y$ is even (the case where $x$ and $y$ are both even can't occur because then $p$ would be divisible by $4$). So $b$ is (the absolute value of) this unique even value divided by $2$, which uniquely determines $a$.
