May $p^3$ divide $(a+b)^p-a^p-b^p$?

We explain the pattern observed by Joe Silverman, deducing the existence of infinitely many solutions, some of which even have $p^5 | (a+b)^p - a^p - b^p$.

Lemma. If $n \equiv 1 \bmod 3$ then the homogeneous polynomial $P_n(a,b) := (a+b)^n - a^n - b^n \in {\bf Z}[a,b]$ has a factor $(a^2+ab+b^2)^2$.

Proof : Either
i) evaluate $P_n(x,1)$ and its derivative at a cube root of unity $\rho$, or
ii) use the $S_3$ symmetry of $-P_n(a,b) = a^n + b^n + c^n$ where $a+b+c=0$: double roots at $(a:b:c) = (1:\rho:\rho^2)$ and $(1:\rho^2:\rho)$ are the only ways to get the total number of roots to be $1 \bmod 3$. $\Box$

[These ${\bf Q}(\rho)$-rational points on the $n$-th Fermat curve are well known.]

Corollary. If $p$ is a prime congruent to $1 \bmod 3$, and $a,b$ are integers such that $p^k | a^2+ab+b^2$, then $p^{2k+1} | P_n(a,b)$.

Proof : Observe that $P_p \in p {\bf Z}[a,b]$ and use the Lemma.

Now $k=1$, and thus $2k+1=3$, is easy to obtain: choose $a$ arbitrarily, and let $b\equiv ra \bmod p$ where $r$ is a cube root of unity mod $p$. We can even get a factor of $p^5$ by choosing $a,b$ such that $a^2 + ab + b^2 = p^2$ (that is, so that $a,b,p$ are sides of a triangle with a $120^\circ$ angle); for example, $$ (3+5)^7 - 3^7 - 5^7 = 120 \cdot 7^5. $$

There are apparently lots of examples. The smallest is $a=1$, $b=2$, and $p=7$.

For what it's worth, consider $b=p-1$. It seems (experimentally at least, I haven't tried to prove it) that for every $p\equiv1\pmod6$ there exists an $a$ in the range $2\le a\le \frac{1}{2}p$ such that $$(a+p-1)^p-a^p-(p-1)^p\equiv0\pmod{p^3}.\qquad (*)$$ I checked for $p<500$, and in this range, for each $p$ there is a unique such $a$. Further, again just for $p<500$, if $p\not\equiv1\pmod6$, then there are no solutions to $(*)$.