Does there exist a polynomial (,) which detects all non-squares? I asked this question in MathStackExchange back in April, and it received more than 30 upvotes, but no answer was offered even after a bounty. I am reposting it here in hopes that someone can answer it.
Problem. Does there exist a two-variable polynomial $P(x, y)$ with integer coefficients such that a positive integer $n$ is not a perfect square if and only if there is a pair $(x, y)$ of positive integers such that $P(x, y)=n$?
Context. The answer is positive for polynomials in 3 variables! This appeared as a problem in USA Team Selection Test in 2013. It turns out that the polynomial $P(x, y, z)=z^2\cdot (x^2-zy^2-1)^2+z$  enjoys the following property: a positive integer $n$ is a not a perfect square if and only if $P(x, y, z)=n$ has a solution in positive integers $(x, y, z)\in \mathbb{N}^{3}$. This construction works nicely due to Pell's equation. If $n$ is not a perfect square, then Pell's equation $x^2-ny^2=1$ has a solution in positive integers $(x_0, y_0)$, and so we get $P(x_0, y_0, n) = n$. Conversely, if $P(x, y, z)=n$, then one can show that $n$ cannot be a perfect square because $n=z^2(x^2-zy^2-1)^2+z$ can be squeezed between two consecutive perfect squares:
$$
 (z(x^2-zy^2-1))^2 < n < (z(|x^2-zy^2-1|+1)^2
$$
Remark. It is clear that there is no single-variable polynomial $P(x)$ which could achieve the desired property. Indeed, there are arbitrary number of consecutive non-squares, and a polynomial $P(x)$ of degree $n>1$ cannot output a consecutive list of $n+1$ numbers. This last claim itself is a nice problem; for a solution, see Example 2.24 in page 11 of Number Theory: Concepts and Problems by Andreescu, Dospinescu and Mushkarov.
 A: This is just a long comment that might be helpful:
Treat $a$ as a parameter, and treat $x$ as a variable.  The Diophantine expression
$$
\exists x\ ((x^2<a) \land (a<(x+1)^2))
$$
defines $a$ as a nonsquare.  Adding two more variables, we can remove the inequalities.  To that end, let
$$
D(a,x,y,z):=(x^2+y-a)^2+(a+z-(x+1)^2)^2.
$$
Then $a$ is a nonsquare exactly when there are positive integers $x,y,z$ such that $D(a,x,y,z)=0$.  Now, for the new polynomial
$$
E(w,x,y,z):=w(1-D(w,x,y,z)^2),
$$
its positive outputs (on positive inputs) are exactly the positive nonsquares (although there are possibly nonpositive outputs).
We can reduce the number of variables needed to define $D$ from three to two by using an idea of Carl Schildkraut (taken from a comment to the MSE version of the question).  Use
$$
D(a,x,y):=((x+y-1)^2+x-a)\cdot ((x+y-1)^2+(x+y-1)+x-a)
$$
instead.  There is a corresponding $3$-variable $E$.
To answer the question fully it would suffice to find a $1$-variable version of $D$.
The hope here is that we could somehow translate Peter Taylor's near miss to this situation, where fractional considerations become less relevant.
