Finding primes so that $x^p+y^p=z^p$ is unsolvable in the $p$-adic units

Question

On my number theory exam yesterday, we had the following interesting problem related to Fermat's last theorem:

Suppose $p>2$ is a prime. Show that $x^p+y^p=z^p$ has a solution in $\mathbb{Z}_p^{\times}$ if and only if there exists an integer $a$ such that $p\not\mid a(a+1)$ and $$(a+1)^p=a^p+1\pmod{p^2}.$$ (I wil not show this here, it's a good exercise in Hensel's lemma)

Now I'm interested for which primes this last property holds true. This is what I found so far:

Proposition: for every prime $p\equiv 1\pmod{3}$ there exists an $a\in\mathbb{Z}$ such that $p\not\mid a(a+1)$ and $p^2\mid(a+1)^p-a^p-1$.
Sketch of proof: We will often use the fact that if $p\mid a-b$ then $p^2\mid a^p-b^p$ for $a,b,p\in\mathbb{Z}$. There exists some $n\in\mathbb{Z}$ such that $p \mid n^2+n+1$. We will show that we can take $a=n$ (it is clear that $p\not \mid n(n+1)$). If $n\equiv 1\pmod{p}$, we are done because $(n+1)^p\equiv 2\equiv n^p+1\pmod{p^2}$. So we may assume that the order of $n$ mod $p$ is 3.

We see that $n^{3p}\equiv 1\pmod{p^2}$, hence $p^2\mid \frac{n^{3p}-1}{n^p-1}=n^{2p}+n^p+1$ and therefore $(n+1)^p\equiv -n^{2p}\equiv n^p+1\pmod{p^2}$. $\blacksquare$

What can we say about primes $p\equiv 2\pmod{3}$? Say a prime is bad if it satisfies the above property. It is easy to see that for $p$ to be good we only have to check that for all $1\le a \le (p-1)/2$ we have $p^2\not\mid (a+1)^p-a^p-1$.
Numerical data suggests that infinitely many primes $p\equiv 2\pmod{3}$ are good and infinitely many are bad. I'd like to find a proof of this.
My numerical data furthermore suggests that there is no easy divisibility criterion to judge badness for primes, and that that there are a lot more good primes than bad ones, in terms of asymptotics.

I welcome all comments/ideas/references.

---

Edit: Here are some observations. I've listed the first bad primes, together with the set of$a\in\mathbb{Z}$ such that $1\le a \le (p-1)/2$ and $(a+1)^p\equiv a^p+1\pmod{p^2}$:

$59\qquad [3, 4, 11, 14, 15, 20] \\ 83\qquad [8, 30, 36] \\ 179\qquad [2, 59, 88] \\ 227\qquad [36, 82, 92] \\ 419\qquad[80, 110, 150] \\ 443\qquad [108, 125, 201] \\ 701\qquad [23, 61, 132, 146, 153, 252] \\ 857\qquad [6, 143, 244] \\ 887\qquad [20, 132, 168] \\ 911\qquad [14, 64, 242] \\ 929 \qquad[234, 253, 266] \\ 971 \qquad[9, 97, 108] \\ 977 \qquad[102, 182, 331] \\ 1091\qquad [64, 234, 358] \\ 1109 \qquad[176, 400, 523] \\ 1193\qquad [224, 227, 473] \\ 1217\qquad [186, 228, 410] \\ 1223\qquad [304, 349, 409] \\ 1259\qquad [19, 62, 264] \\ 1283\qquad [19, 134, 449] \\ 1289\qquad [412, 535, 593] \\ 1439\qquad [93, 199, 294] \\ 1487\qquad [167, 416, 649] \\ 1493 \qquad[141, 179, 367] \\ 1613 \qquad[227, 473, 739] \\ 1637 \qquad[76, 279, 574] \\ 1811 \qquad[39, 123, 498, 628, 691, 743] \\ 1847\qquad [172, 362, 698] \\ 1901\qquad [3, 475, 634] \\ 1997\qquad [125, 671, 840] \\ 2003 \qquad[31, 312, 840] \\ 2087 \qquad[31, 202, 586] \\ 2243 \qquad[252, 593, 1077] \\ 2423\qquad [353, 704, 857] \\ 2477\qquad [17, 688, 1020] \\ 2579 \qquad[294, 576, 693] \\ 2591 \qquad[322, 345, 368] \\ 2729\qquad [1024, 1049, 1089] \\ 2777\qquad [488, 1078, 1269] \\ 2969\qquad [341, 789, 1158] \\ 3089\qquad [677, 1015, 1053]\\ 3137\qquad [1227, 1308, 1427]\\ 3191\qquad [776, 916, 1383]\\ 3203\qquad [1214, 1305, 1587]\\ 3251\qquad [164, 337, 1260]$

I also found this sequence on OEIS, along with more numerical data, where it is claimed (but not proved) that the density of bad primes is roughly $\frac 16$.

My numerical data also suggests that the number of $a\in\mathbb{Z}$ satisfying $1\le a \le (p-1)/2$ and $p^2\mid (a+1)^p-a^p-1$ is always either $3$ or $6$, but I think that's rather a bold statement. (edit: indeed, 20123 is a counterexample with $9$ solutions, found by Michael Stocker. The weaker statement that the number of such solutions is divisible by $3$ is true, however: see here.)

---

Edit: Someone showed me these slides, where the above criterion for solvability of $x^p+y^p=z^p$ in the $p$-adic units is stated. (There is a remark that it can be proven with eisenstein reciprocity that for all non-Wieferich primes the first case of FLT holds true - maybe there is some way to use this for some (partial) result on our problem? I'm really curious.)

Answer 1

Generating solutions for Fermat's Last Theorem modulo $p^2$ for prime $59$:

Assume:

$$x^p=y^p+z^p$$ $$p \not \| xyz$$ $$gcd(x,y,z)=1$$ $$x^{p-1} \equiv y^{p-1} \equiv z^{p-1} \equiv 1 \pmod{p^2}$$

We generate solutions for:

$$(y+z)^p \equiv y^p+z^p \pmod{p^2}$$

Note:

$$\frac{x^{p+1}-y^{p+1}}{x-y} \equiv (p+1)y^p + a(x-y) \equiv x+y \pmod{p^2}$$ $$\frac{x^{p+1}-z^{p+1}}{x-z} \equiv (p+1)z^p + b(x-z) \equiv x+z \pmod{p^2}$$

$ \implies $

$$az \equiv x-py \pmod{p^2}$$ $$by \equiv x-pz \pmod{p^2}$$

$ \implies $

$$z(a-1) \equiv y(1-p) \pmod{p^2}$$ $$y(b-1) \equiv z(1-p) \pmod{p^2}$$ $$z(a-p) \equiv y(b-p) \equiv x(1-p) \equiv (y+z)(1-p) \pmod{p^2}$$

$ \implies $

$$(a-1)^{p-1} \equiv (b-1)^{p-1} \equiv (a-p)^{p-1} \equiv (b-p)^{p-1} \equiv 1+p \pmod{p^2}$$

Note:

$$(a-1)(b-1) \equiv (1-p)^2 \equiv 1-2p \pmod{p^2}$$

$ \implies $

$$ab \equiv a+b-2p \pmod{p^2}$$

To generate our solutions modulo 59^2 we first find solutions for $a$ in:

$$(a-1)^{p-1} \equiv (a-p)^{p-1} \equiv 1+p \pmod{59^2}$$

$$\{64,122, 1428,1491,1496,1549,1992,2045,2050,2113,3419,3477\}$$

For out exapmle we take $a=64$.

We find $b$ by calculating:

$$ab \equiv a+b-2p \pmod{p^2}$$

$$a=64$$ $$b=1491$$

We take $r=2$ in $x \equiv r^{59} \equiv 2^{59} \equiv 946 \pmod{p^2}$

We now use $a,b,x$ to solve $y,z$ in:

$$z(a-p) \equiv y(b-p) \equiv x(1-p) \pmod{p^2}$$

$$y \equiv 2869 \pmod{p^2}$$ $$z \equiv 1558 \pmod{p^2}$$

We now have $x,y,z$ with:

$$x^{p-1} \equiv y^{p-1} \equiv z^{p-1} \equiv 1 \pmod{p^2}$$

and:

$$x^p \equiv (y+z)^p \equiv y^p + z^p \equiv (2869 + 1558)^{59} \equiv 4427^{59} \equiv 2869^{59} + 1558^{59} \pmod{59^2}$$

Running $r$ through $1..58$ generates $58$ solutions.

Note:

$$\frac{\frac{x^{p+1}-y^{p+1}}{x-y} - (p+1)y^p}{x-y} \equiv \frac{\frac{946^{59}-2869^{59}}{946-2869}-60(2869)^{59}}{946-2869} \equiv a \equiv 64 \pmod{p^2}$$ $$\frac{\frac{x^{p+1}-z^{p+1}}{x-z} - (p+1)z^p}{x-z} \equiv \frac{\frac{946^{59}-1558^{59}}{946-1558}-60(1558)^{59}}{946-1558} \equiv b \equiv 1491 \pmod{p^2}$$