$p = x^2 + xy + y^2$ if and only if $p \equiv 1 \text{ mod }3$?

$p=x^2+xy+y^2=(x-y)^2+3xy$, so this is about the behaviour of perfect squares. If $z$ is an integer not divisible by $3$, then $z^2\equiv 1\mod 3.$ Since $p\ne 3$ is prime, it is not divisible by $3$. Hence $x-y$ is not divisible by $3$, and it follows that $p\equiv 1\mod 3.$

Let $K = \mathbb Q(\sqrt{-3})$. Consider the integer ring $\mathcal O_K = \mathbb Z[\alpha]$, where $\alpha = \frac{1+\sqrt{-3}}{2}$. Recall that the field norm of an element $x+y\alpha\in\mathcal O_K$ is $$N_{K/\mathbb Q}(x+y\alpha) =(x+y\alpha)(x+y\overline\alpha) = x^2+xy+y^2.$$

You've correctly remarked that $\mathcal O_K$ is a PID. It follows that

A prime $p$ can be expressed as $p=x^2+xy+y^2$ if and only if $p$ is the norm of some element of $\mathcal O_K$.

If $p\ne 3$ (so that $p$ does not ramify in $\mathcal O_K$) and $p = N_{K/\mathbb Q}(x+y\alpha)$, then the ideal $(x+y\alpha)$ is a prime ideal of $\mathcal O_K$ lying above $p$ with norm $p$. So $p$ splits completely in $\mathcal O_K$.

Conversely, if $p = \mathfrak{p_1p_2}$, then since $\mathcal O_K$ is a PID, $\mathfrak p_1$ is generated by an element of $\mathcal O_K$, which must have norm $p$. Hence,

A prime $p\ne 3$ can be expressed as $p=x^2+xy+y^2$ if and only if $p$ splits completely in $\mathcal O_K$.

To examine when this can happen, we can use the following version of the Kummer-Dedekind Theorem:

Theorem: Let $p$ be a prime, and let $\beta\in \mathcal O_K$ be such that $K=\mathbb Q(\beta)$ and $p\nmid (\mathcal O_K:\mathbb Z[\beta])$. Let $f(X)$ be the minimal polynomial of $\beta$ over $\mathbb Q$. Suppose that $$f(X) \equiv f_1(X)^{e_1}\cdots f_m(X)^{e_m}\pmod p. $$ Then $p$ splits as $p\mathcal O_K = \mathfrak{p_1^{e_1}\cdots p_m^{e_m}}$ in $\mathcal O_K$.

In particular, taking $K$ as above and $\beta = \sqrt{-3}$, since $(\mathcal O_K, \mathbb Z[\sqrt{-3}]) = 2$, the above theorem applies to all primes $p>2$. The case $p=2$ can be checked manually.

So for $p\ge 5$ $$ \begin{align} p \text{ splits completely in }\mathcal O_K&\iff X^2+3\text{ splits into distinct linear factors mod } p\\&\iff X^2+3\text{ is reducible mod } p\\ &\iff\left(\frac{-3}p\right)=1\\ &\iff(-1)^{(p-1)/2}\left(\frac{3}p\right)=1\\ &\iff(-1)^{\frac{p-1}2}(-1)^{\frac{p-1}2\frac{3-1}2}\left(\frac p3\right) = 1\text{ by quadratic reciprocity}\\ &\iff \left(\frac p3\right) = 1\\&\iff p\equiv 1 \pmod 3. \end{align} $$ The result follows.

Here is an answer, maybe not the most illuminating for you and probably not that different from Mathmo123's answer, but involving Jacobi sums.

Let $p\equiv 1 ($mod $3$), $\lambda$ a generator of the group of characters over $\mathbb{F}_p^{\times}$ and set $\chi=\lambda^{\frac{p-1}{3}}$. So the character $\chi$ takes it values in $V=\{0,1,\alpha,\alpha^2\}$ where $\alpha=\frac{-1+\sqrt(-3)}{2}$. Given that the set $V$ is stable under multiplication, the Jacobi sum associated to $\chi$, $J(\chi,\chi)=\sum_{a+b=1} \chi(a)\chi(b)$ is in $\mathbb{Z}+\alpha\mathbb{Z}+\alpha^2\mathbb{Z}=\mathbb{Z}+\alpha\mathbb{Z}$. So we can write $J(\chi,\chi)=x+y\alpha$ with $x$ and $y$ in $\mathbb{Z}$. But we know $p=|J(\chi,\chi)|^2$; so $p=x^2-xy+y^2$.

Of course the same strategy works for $p\equiv 1$ mod $4$ $\Longleftrightarrow p=x^2+y^2$.