Lagrange proved that every nonnegative integer is a sum of 4 squares.
Gauss proved that every nonnegative integer is a sum of 3 triangular numbers.
Is there a 2-variable polynomial $f(x,y) \in \mathbf{Q}[x,y]$ such that $f(\mathbf{Z} \times \mathbf{Z})=\mathbf{N}$?
This is a cute problem! I toyed with it and didn't really get anywhere - I got the strong impression that it requires fields of mathematics that I am not expert in.
Indeed, given that the problem seems related to that of counting integer solutions to the equation $f(x,y) = c$, one may need to use arithmetic geometry tools (e.g. Faltings' theorem). In particular if we could reduce to the case when the genus is just 0 or 1 then presumably one could kill off the problem. (One appealing feature of this approach is that arithmetic geometry quantities such as the genus are automatically invariant (I think) with respect to invertible polynomial changes of variable such as $(x,y) \mapsto (x,y+P(x))$ or $(x,y) \mapsto (x+Q(y),y)$ and so seem to be well adapted to the problem at hand, whereas arguments based on the raw degree of the polynomial might not be.)
Of course, Faltings' theorem is ineffective, and so might not be directly usable, but perhaps some variant of it (particularly concerning the dependence on c) could be helpful. [Also, it is overkill - it controls rational solutions, and we only care here about integer ones.] This is far outside of my own area of expertise, though...
The other thing that occurred to me is that for fixed c and large x, y, one can invert the equation $f(x,y) = c$ to obtain a Puiseux series expansion for y in terms of x or vice versa (this seems related to resolution of singularities at infinity, though again I am not an expert on that topic; certainly Newton polytopes seem to be involved). In some cases (if the exponents in this series expansion are favourable) one could then use Archimedean counting arguments to show that f cannot cover all the natural numbers (this is a generalisation of the easy counting argument that shows that a 1D polynomial of degree 2 or more cannot cover a positive density set of integers), but this does not seem to work in all cases, and one may also have to use some p-adic machinery to handle the other cases. One argument against this approach though is that it does not seem to behave well with respect to invertible polynomial changes of variable, unless one works a lot with geometrical invariants.
Anyway, to summarise, it seems to me that one has to break out the arithmetic geometry and algebraic geometry tools. (Real algebraic geometry may also be needed, in order to fully exploit the positivity, though it is also possible that positivity is largely a red herring, needed to finish off the low genus case, but not necessary for high genus, except perhaps to ensure that certain key exponents are even.)
EDIT: It occurred to me that the polynomial $f(x,y)-c$ might not be irreducible, so there may be multiple components to the associated algebraic curve, each with a different genus, but presumably this is something one can deal with. Also, the geometry of this curve may degenerate for special c, but is presumably stable for "generic" c (or maybe even all but finitely many c).
It also occurs to me that one use of real algebraic geometry here is to try to express f as something like a sum of squares. If there are at least two nontrivial squares in such a representation, then f is only small when both of the square factors are small, which is a 0-dimensional set and so one may then be able to use counting arguments to conclude that one does not have enough space to cover all the natural numbers (provided that the factors are sufficiently "nonlinear"; if for instance $f(x,y)=x^2+y^2$ then the counting arguments barely fail to provide an obstruction, one has to use mod p arguments or something to finish it off...)
EDIT, FOUR YEARS LATER: OK, now I know a bit more arithmetic geometry and can add to some of my previous statements. Firstly, it's not Faltings' theorem that is the most relevant, but rather Siegel's theorem on integer points on curves - the enemy appears to be those points $(x,y)$ where $x,y$ are far larger than $f(x,y)$, and Siegel's theorem is one of the few tools available to exclude this case. The known proofs of this theorem are based on two families of results in Diophantine geometry: one is the Thue-Siegel-Roth theorem and its variants (particularly the subspace theorem), and the other is the Mordell-Weil theorem and its variants (particularly the Chevalley-Weil theorem). A big problem here is that all of these theorems have a lot of ineffectivity in them. Even for the very concrete case of Hall's conjecture on lower bounding $|x^2-y^3|$ for integers $x,y$ with $x^2 \neq y^3$, Siegel's theorem implies that this bound goes to infinity as $x,y \to \infty$, but provides no rate; as I understand it, the only known lower bounds are logarithmic and come from variants of Baker's method.
As such, a polynomial such as $f(x,y) = (x^2 - y^3 - y)^4 - y + C$ for some large constant C already looks very tough to analyse. (I've shifted $y^3$ here by $y$ to avoid the degenerate solutions to $x^2=y^3$, and to avoid some cheap way to deal with this polynomial from the abc conjecture or something.) The analogue of Hall's conjecture for $|x^2-y^3-y|$ suggests that $f(x,y)$ goes to $+\infty$ as $x,y \to \infty$ (restricting $x,y$ to be integers), but we have no known growth rate here due to all the ineffectivity. As such, we can't unconditionally rule out the possibility of an infinite number of very large pairs $(x,y)$ for which $x^2-y^3-y$ happens to be so close to $y^{1/4}$ that we manage to hit every positive integer value in $f(x,y)$ without hitting any negative ones. However, one may be able to get a conditional result assuming some sufficiently strong variant of the abc conjecture. One should also be able to exclude large classes of polynomials $f$ from working; for instance, if the curve $f(x,y)=0$ meets the line at infinity at a lot of points in a transverse manner, then it seems that the subspace theorem may be able to get polynomial bounds on solutions $(x,y)$ to $f(x,y)=c$ in terms of $c$, at which point a lot of other tools (e.g. equidistribution theory) become available.
Another minor addendum to my previous remarks: the generic irreducibility of $f(x,y)-c$ follows from Bertini's second theorem, as one may easily reduce to the case when $f$ is non-composite (not the composition of two polynomials of lower degree).
The search turned up a 1981 paper by John S.Lew (in the Unsolved problems section)
which discusses related problems, and ends up stating this one. The author's problems are:
His main conjecture is that the only solutions to A are Cantor's $x+ \frac12(x+y-1)(x+y-2)$, which apparently goes to the time of Polya. Lew states C independently from empirical observations.
What can be said about the following stronger question? Let $f(x,y)\in \mathbf{Q}[x,y]$ such that $f(\mathbf{Z}\times \mathbf{Z})$ is a subset of $\mathbf{N}$. Let $g(n)$ be the number of elements of $f(\mathbf{Z}\times \mathbf{Z})\cap \lbrace 0,1,\dots,n\rbrace$. How fast can $g(n)$ grow? Is it always true that $g(n) =O(n/\sqrt{\log(n)})$? If true this is best possible since if $f(x,y)=x^2+y^2$ then $g(n)\sim cn/\sqrt{\log(n)}$.
After thinking about this problem for a bit using rather a naive approach, looking at regions where f grows faster than quadratically (as mentioned in Qiaochu's attempt), it certainly appears that obtaining a negative answer to this problem is very difficult. To obtain a positive answer might be easier since we only need to exhibit a single polynomial with $f(\mathbb{Z}\times\mathbb{Z})=\mathbb{N}$ although, in all likelihood, there do not exist such examples. However, even using the naïve approach, it quickly becomes clear what kind of behaviour could lead to a polynomial $f$ having the required properties. Polynomials of degree less than four can be quickly dismissed. Then, assuming $f$ has degree $2n$, look at the leading order terms $f_{2n}$, which is a homogeneous polynomial of degree $2n$. Away from the zeros of $f_{2n}$, it grows at rate $R^{2n}$ ($R=\sqrt{x^2+y^2}$) and dominates the lower order terms, so $f$ cannot cover a strictly positive density of the integers. The difficulties occur when we look close to certain curves which are asymptotic to the zeros of $f_{2n}$ (being homogeneous, the zeros of $f_{2n}$ lie on a finite set of lines radiating out from the origin). On such curves, $f_{2n}$ will grow at a rate less than $O(R^{2n})$ and cancellation with lower order terms can occur. Looking at how much cancellation can occur seems to lead to difficult problems of Diophantine approximation.
The kind of polynomial which are plausible candidates for a polynomial mapping onto the positive integers are as follows. $$ f(x,y)=a\prod_{i=1}^d\left(x-\alpha_iy\right)^{2n} - q(x,y)\tag{1}\label{47389_1} $$ where $\alpha$ is a real algebraic integer with minimal polynomial of degree $d > 2$ over the rationals and conjugates $\alpha_1,\dotsc,\alpha_d$, $a$ is a positive integer, and $q(x,y)$ is a polynomial of degree $2n(d-2)$. By Dirichlet's theorem, we know that there are infinitely many integer x,y such that $\lvert x/y-\alpha_i\rvert < y^{-2}$, so the leading order term of \eqref{47389_1} is less than some fixed multiple of $R^{2n(d-2)}$ infinitely often. So there will be some cancellation with q. On the other hand, by the Thue–Siegel–Roth theorem, we know that $\lvert x/y-\alpha\rvert > y^{-2-\epsilon}$ for all large x,y, so the leading order terms of \eqref{47389_1} grow at least as fast as $O(R^{2n(d-2)-\epsilon})$ which, at least, means that \eqref{47389_1} cannot go negative very quickly. The question then, is there an algebraic number $\alpha$ such that $\lvert x/y-\alpha\rvert\ge cy^{-2}$ for some positive constant $c$ and all integers $x$, $y$? In that case, the leading order term of \eqref{47389_1} would be bounded below by a multiple $R^{2n(d-2)}$ and $q$ could be chosen such that $f(\mathbb{Z}\times\mathbb{Z})\subseteq\mathbb{N}$. It is then possible, but still unlikely, that \eqref{47389_1} gives a polynomial mapping onto the positive integers. Looking in my copy of Hindry & Silverman (Diophantine Geometry, An Introduction) it mentions that it is an open problem whether there exist such algebraic numbers (and there are no known examples or counterexamples with $d > 2$) but it is conjectured that there aren't any. So, polynomials such as \eqref{47389_1} appear unlikely to do what we want, but proving this seems to be very difficult. Of course, that $f$ actually covers $\mathbb{N}$ would be a much stronger statement than $\lvert x/y-\alpha\rvert \ge cy^{-2}$, so maybe an expert in this area could actually rule out such examples, but it still looks like a very tricky problem.
We can also try polynomials such as $$ f(x,y)=a\prod_{i=1}^d\left((x-\alpha_iy)^r-p(\alpha_i)y^s\right)^{2n} - q(x,y)\tag{2}\label{47389_2} $$ where, now, $p$ is a polynomial with integer coefficients, $r$, $s$ are positive integers and $q$ has degree less than $2n$. This is even more difficult than \eqref{47389_1} to deal with and whether or not such polynomials can provide what the question is asking for depends on how small $(x/y-\alpha)-p(\alpha)^{1/r}y^{s/r-1}$ can be. This also looks like a very difficult problem in Diophantine approximation.
So, although I expect that the answer to this is no, there are no such $f$, any method of proving this has to cope with possibilities such as \eqref{47389_1} and \eqref{47389_2}. Just these two cases look extremely difficult to handle. Maybe it is possible though, and an expert on such areas would be able to say something more about them than I can.
I don't think so (but I haven't checked this argument very thoroughly):
First we claim that any such $f$ has degree two. Clearly the leading term of $f$ cannot be odd, so suppose by contradiction that $f$ has degree at least four. Pick a constant $R$ large enough so that in the region $D_1$ consisting of points satisfying $|x|, |y| \ge R$, there exists a constant $c$ such that $f(x, y) \ge c \text{ min}(|x|, |y|)^4$. (Edit: This isn't always possible, but I think it can be salvaged.) Then it's not hard to see that $\sum_{D_1} \frac{1}{f(x, y)}$ converges. But $D_2 = \mathbb{Z}^2 - D_1$ can be partitioned into $4R - 2$ not necessarily disjoint lines in which one of $x$ or $y$ is fixed. On any of these regions $f$ cannot be linear, so it either grows at least quadratically or is a constant; we can ignore the lines on which $f$ is constant. It follows that $\sum_L \frac{1}{f(x, y)}$ converges for any line $L$ on which $f$ is nonconstant, hence $\sum \frac{1}{f(x, y)}$ converges if we sum over every point of $\mathbb{Z}^2$ except the lines on which $f$ is constant. Since we have only thrown out finitely many of the values of $f$ in this sum, those values cannot contain every positive integer.
But if $f$ is quadratic, it is a constant plus the sum of squares of two polynomials with rational coefficients and there are many integers not representable as the sum of squares of two rational numbers.
How about an alternative question: does there exist a polynomial $f\in\mathbb{Q}[x,y]$ with integer values at lattice points, and of degree at least two on each variable, such that for any prime $p$, the map $f:\mathbb{Z}\times\mathbb{Z}\to\mathbb{Z}$ composing with $\mathbb{Z}\to\mathbb{F}_p$ is surjective? Or more specifically, does a degree two such polynomial exist? The last part shouldn't be too hard, but I don't know how to solve it.
If we can make a (single variable) polynomial function $g(x)$ from $\mathbb{Z}$ onto $\mathbb{N}$, we could compose it with the Cantor pairing function, but such a $g(x)$ seems implausible for some reason...
In this pre-print, I showed that the existence of such a polynomial would contradict Vojta's conjecture. The rough outline of the argument is as follows.
Suppose such an $F \in \mathbb{Q}[x, y]$ exists. The surfaces $X_n : z^n = F(x, y)$ will then have infinitely many integer points (since every positive $n$th power is represented by $F$).
If the integral points on $X_n$ are not (Zariski) dense for some odd $n$, then they lie on finitely many curves. By increasing $n$, we can reduce this to rational curves with a single puncture. However, for such curves either the integer values of $z$ represented can either be both positive and negative (contradicting that the image of $F$ is contained in $\mathbb{N}$) or must all be perfect squares (contradicting surjectivity).
If instead, for infinitely many $n$, $X_n$ is not of log general type, then this shows the surface $\mathbb{A}^2 \setminus \{F(x, y) = 0\}$ is also not of log general type. There is a classification of such $F$ (e.g. this paper of Aoki), that $F(x, y) = c$ is a rational curve for almost all $c$, and either the normalisation of this curve is missing at most 2 points, or there exists an auxiliary polynomial with a similar set of restrictions. They also give a canonical form (over $\overline{\mathbb{Q}}$) for such functions. As $F(x, y) = c$ is rational and has a rational point for all $c \in \mathbb{N}$ forces $F(x, y) = T$ to be a rational curve over $\mathbb{Q}(T)$ (as was pointed out to me by Samir Siksek), and so $F : \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ is surjective. By considering each case in the classification via a weaker version of the canonical form, I show that the set $\{(x, y) \in \mathbb{R}^2 : F(x, y) < 0\}$ is large enough to contain integral points, and so $F(\mathbb{Z}, \mathbb{Z})$ is not contained in $\mathbb{N}$.
We therefore conclude that $X_n$ must be of log general type and has dense integral points. This contradicts Vojta's conjecture (here, taken as an integral generalisation of the Bombieri-Lang conjecture).
