Does there exist a polynomial 𝑃(𝑥,𝑦) which detects all non-squares?

Question

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.

Answer 1

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.