Already for three-term progressions it's somewhat surprising that
there are infinitely many solutions, because the usual probabilistic
guess for the expected number of solutions leads to a convergent sum:
a random number of size $M$ is a square with probability about
$M^{-1/2}$, so we're summing something like $1/(abc)$ over all three-term progressions
$(a,b,c)$, etc. To be sure such a guess cannot account for non-random
patterns arising from polynomial identities, but it does suggest that
past a certain point such identities will be the only source of solutions.

Now a mindless exhaustive search over progressions
$(x,x+y,x+2y)$ with $0 < x,y < 10^4$ finds only the first six examples
$$
(1,7),\phantom+
(4,26),\phantom+
(15,97),\phantom+
(56,362),\phantom+
(209,1351),\phantom+
(780,5042)
$$
of an infinite family associated with the solutions
$(2,1)$, $(7,4)$, $(26,15)$, $(97,56)$, $(362,209)$, $(1351,780)$, etc.
of the Pell equation $x^2-3y^2=1$. If it can be proved that these are
the only solutions then it will immediately follow that there are
no four-term arithmetic progressions with the same property.
But that seems like a very hard problem.

Here's the **gp** code; with a bound of $10^4$ it takes only
a few minutes. One can surely do better with a more intelligent
search procedure (e.g. start by finding all solutions of $ab+1=r^2$
by factoring $r^2-1$).

```
H = 10^4
progsq(x,y,n) = sum(i=0,n-2,sum(j=i+1,n-1,issquare((x+i*y)*(x+j*y)+1)))
for(x=1,H,for(y=1,H,if(progsq(x,y,3)==3,print([x,y]))))
```