Prime numbers which divide $n^3-3n+1$
Solution 1:
A "problem" with a fully elementary proof is that concepts are often pulled out of nowhere. This is already seen in your proof, since you do not explain why you consider writing $n \equiv a + 1/a \bmod p$. How did you decide to try to solve such a congruence? Why are you interested in a proof making no use of groups or fields?
Motivation for the congruence $n \equiv a + 1/a \bmod p$ comes from looking in $\mathbf R$. The polynomial $x^3 - 3x + 1$ has three real roots: $2\cos(2\pi/9)$, $2\cos(4\pi/9)$, and $2\cos(8\pi/9)$. In terms of the primitive 9th root of unity $z = e^{2\pi i/9}$, there are 6 primitive 9th roots of unity $z$, $z^2$, $z^4$, $z^5 = 1/z^4$, $z^7 = 1/z^2$, and $z^8 = 1/z$, and the three roots of $x^3 - 3x + 1$ are the sum of a primitive 9th root of unity and its complex conjugate: $2\cos(2\pi/9) = z + \overline{z} = z + 1/z$, $2\cos(4\pi/9) = z^2 + \overline{z^2} = z^2 + 1/z^2$, and $2\cos(8\pi/9) = z^4 + \overline{z^4} = z^4 + 1/z^4$. This suggests that when $x^3 - 3x + 1$ has a root $n \bmod p$, $n \bmod p$ should have the form $a + 1/a$ for a 9th root of unity $a$ in characteristic $p$. That is analogous to $2\cos(2\pi/9) = z + 1/z$ above. The intuition behind $a^6 + a^3 + 1 \equiv 0 \bmod p$ comes from $a \bmod p$ being a primitive 9th root of unity in characteristic $p$ and $x^9 - 1$ factoring like this: $$ x^9 - 1 = (x-1)(x^2+x+1)(x^6 + x^3 + 1). $$ In this factorization, the root of $x - 1$ is $1$ and the roots of $x^2+x+1$ are the nontrivial cube roots of unity, so the roots of $x^6 + x^3 + 1$ are the primitive 9th roots of unity.
The condition $p \equiv \pm 1 \bmod 9$ that you are interested in is the same as $p^2 \equiv 1 \bmod 9$, and that very strongly suggests looking at the field $\mathbf F_{p^2}$ and showing its group of nonzero elements contains an element of order $9$ (primitive 9th root of unity in characteristic $p$). Then $9 \mid (p^2-1)$ by Lagrange's theorem. But you want to avoid group theory. I don't see how to derive $p^2 \equiv 1 \bmod 9$ in a nice way by working only in $\mathbf F_p^\times$ instead of in $\mathbf F_{p^2}^\times$.
If $r$ is one of the roots of $x^3-3x+1$ in $\mathbf R$ then the other two roots are $r^2-2$ and $-r^2-r+2$. This has an analogue in $\mathbf F_p$ if $p \not= 3$ and $n^3 - 3n + 1 \equiv 0 \bmod p$ for an integer $n$: $x^3 - 3x + 1$ has three different roots in $\mathbf F_p$, namely $n \bmod p$, $n^2 - 2 \bmod p$, and $-n^2 - n + 2 \bmod p$. (When $p = 3$ these roots are all equal.) Maybe you can get something useful by working with all three of the roots in $\mathbf F_p$ and not just one of the roots mod $p$.
[Update: since you know some field theory already, here is how I would solve this problem using field theory, rather than the method using more machinery below. The key insight, as mentioned above, is that $x^3 - 3x + 1$ has roots of the form $z + 1/z$ as $z$ runs over primitive 9th roots of unity. This is true in $\mathbf C$ and in fact it's true in every field where there are 9 ninth roots of unity: if $z^9 = 1$ and $z^3 \not= 1$ then $z + 1/z$ is a root of $x^3 -3x + 1$. For a prime $p \not= 3$, the polynomial $x^9 - 1$ has nine distinct roots in characteristic $p$ and only 3 of those roots are cube roots of unity, so $x^9-1$ has 6 roots in characteristic $p$ that are primitive 9th roots of unity. Pairing the 6 roots together as $z$ and $1/z$, the three sums $z + 1/z$ are the three roots of $x^3 - 3x + 1$. All 3 roots can be expressed in terms of any one root $r$: the roots are $r$, $r^2 - 2$, and $-r^2 - r + 2$, so a field containing one root of $x^3 - 3x + 1$ contains all three of its roots. Let $\zeta$ be a primitive 9th root of unity in characteristic $p$. It does not have to lie in $\mathbf F_p$, but $\mathbf F_p(\zeta)$ is a finite field since $\zeta$ is a root of $x^6 + x^3 + 1$. Since $\zeta + 1/\zeta$ is a root of $x^3 - 3x + 1$ and a field containing one root of this polynomial contains all three roots, your question is about showing, for $p \not= 3$, that $\zeta + 1/\zeta$ lies in $\mathbf F_p$ if and only if $p \equiv \pm 1 \bmod 9$.
("If") In a field $K$ of characteristic $p$, the elements of its subfield $\mathbf F_p$ are the solutions in $K$ to $a^p = a$: that polynomial equation has at most $p$ solutions in $K$ and all $p$ elements of $\mathbf F_p$ are solutions of that equation (the $p-1$ nonzero elements of $\mathbf F_p$ satisfy $a^{p-1} = 1$, so $a^p = a$, and $a = 0$ also fits $a^p = a$). If $p \equiv \pm 1 \bmod 9$ then $\zeta + 1/\zeta \in \mathbf F_p$ since $\zeta + 1/\zeta$ equals its own $p$th power: $$ (\zeta + 1/\zeta)^p = \zeta^p + (1/\zeta)^p = \zeta^p + 1/\zeta^p = \zeta + 1/\zeta $$ since $\zeta^p$ only depends on $p \bmod 9$ and thus $\zeta^p$ is $\zeta$ or $1/\zeta$ from $p \equiv \pm 1 \bmod 9$.
("Only if") The field extension $\mathbf F_p(\zeta)/\mathbf F_p(\zeta + 1/\zeta)$ has degree at most $2$ since $\zeta$ is a root of $x^2 - (\zeta + 1/\zeta)x + 1$. If $\zeta + 1/\zeta$ lies in $\mathbf F_p$, then $\mathbf F_p(\zeta)/\mathbf F_p$ has degree at most 2, so $\zeta \in \mathbf F_{p^2}$ (there is exactly one quadratic extension field of $\mathbf F_p$ and it has size $p^2$). That implies $\zeta^{p^2-1} = 1$, since $\mathbf F_{p^2}^\times$ is a multiplicative group of order $p^2-1$. Since $\zeta$ has multiplicative order 9 (it is a primitive 9th root of unity!), from $\zeta^{p^2-1} = 1$ we must have $9 \mid (p^2-1)$, so $p^2 \equiv 1 \bmod 9$ and thus $p \equiv \pm 1 \bmod 9$.
This concludes my update to this answer.]
I do not see a fully elementary proof for what you seek, but I do see a natural proof if the reader knows algebraic number theory. First of all, the motivation for knowing that all the roots of $x^3-3x+1$ are expressible in terms of one of them comes from Galois theory since $x^3 - 3x+1$ is irreducible over $\mathbf Q$ with a square discriminant $81$. That the discriminant is $81$ tells us by algebraic number theory that the only ramifying prime in $\mathbf Q(r)$ is 3. The Kronecker-Weber theorem tells us the abelian Galois extension $\mathbf Q(r)/\mathbf Q$ has to lie in some cyclotomic field, and we are led to check how $\mathbf Q(r)$ is related to the 9th cyclotomic field because the only primes at which the $m$th cyclotomic field ramifies are factors of $m$. So $\mathbf Q(r)$ should lie in a 3-power cyclotomic field. In fact, $\mathbf Q(r)$ is the real subfield of $\mathbf Q(\zeta_9)$, where $\zeta_9 = e^{2\pi i/9}$. The connection between prime splitting in Galois extensions of $\mathbf Q$ and Frobenius elements in Galois groups over $\mathbf Q$ tells us $p$ splits in $\mathbf Q(r)$ if and only if $p \equiv \pm 1 \bmod 9$ because under the natural isomorphism of ${\rm Gal}(\mathbf Q(\zeta_9)/\mathbf Q)$ with $(\mathbf Z/9\mathbf Z)^\times$, the real subfield $\mathbf Q(r) = \mathbf Q(\zeta_9 + 1/\zeta_9)$ corresponds to the subgroup $\pm 1 \bmod 9$. This, to me, explains why $x^3 - 3x + 1$ has a root mod $p$ if and only if $p=3$ or $p \equiv \pm 1 \bmod 9$ in a very conceptual way. Of course it does not at all qualify as an "elementary proof" making no use of abstract algebra. Sorry.