Nonnegativity conditions for a polynomial in two variables?

There is Stengle's positivensatz,, otherwise, if your polynomial were homogeneous, you could use Polya's theorem, and if you could make your domain compact, you can use Handelman's theorem, but otherwise, Stengle is your man.

For Polya and Handelman see Powers and Reznick:

\bib{MR1854339}{article}{
  author={Powers, Victoria},
  author={Reznick, Bruce},
  title={A new bound for P\'olya's theorem with applications to polynomials
  positive on polyhedra},
  note={Effective methods in algebraic geometry (Bath, 2000)},
  journal={J. Pure Appl. Algebra},
  volume={164},
  date={2001},
  number={1-2},
  pages={221--229},
  issn={0022-4049},
  review={\MR{1854339 (2002g:14087)}},
  doi={10.1016/S0022-4049(00)00155-9},

}

You can rewrite your polynomial in squared variables (so it will be of degree 8) and homogenise it, obtaining a ternary form of degree 8.

There is a result of Hilbert in Acta Math. 17(1893) that applies here; he shows that a real ternary form of degree $2d$ satisfies $f(x,y,z)\geq 0$ for all $x,y,z$ iff it is representable as the rational function $p(x,y,z)/q(x,y,z)$, with $p$ and $q$ homogeneous globally nonnegative of degrees $4d-4$, resp., $2d-4$, and $p$ being a sum of squares of forms. (the result also says how many terms there are, etc)

Here, as $d=4$, we can assume that $q$ is a sum of squares of forms, as it is of degree 4 (by much more famous Hilbert's result from 1888).

We showed that this can be made into an algorithm, cf. http://www.optimization-online.org/DB_HTML/2001/11/399.html (published here: http://www.sciencedirect.com/science/article/pii/S0377221703005162)

Here's a trick (which sometimes works), but it might be contained in the comments above. As stated, the Newton polyhedron (convex hull of the exponents, the latter viewed as lattice points in ${\bf R}^2$), $K$, of $P$, is contained in the double unit square in ${\bf R}^2$. First, an obvious necessary condition for nonnegativity: the coefficients of the extreme points of $K$ must be positive (as opposed to negative).

Next, for each edge of $K$, there is a corresponding polynomial (just take the exponents appearing in the edge, and use the coefficients from $P$); if $P$ is nonnegative, each of these edge polynomials must be nonnegative---but after a change of variables, the edge polynomials are really just polynomials in one variable, for which testing for nonnegativity is just the quadratic formula in your example).

This gives some necessary conditions.

Now a sufficient condition for positivity: If $P(x,y) >$ on the interior of the positive orthant, and each of the edge polynomials is positive (which implies the coefficients at the extreme points of $K$ are all positive), then for any $Q$ with the same supporting monomials as $P$ and $Q$ having no nonnegative coefficients, there exists an integer $M$ such that the product, $Q^M P$, has no negative coefficients. (This is a generalization of Polya and Meissner's theorems; when $K$ is the double unit square, this is Meissner's theorem.)

The simplest (but probably not the best) choice for $Q$ is $\sum x^w$ where $w$ varies over the lattice points in $K$ (where $x^w = x_1^{w(1)}x_2^{w(2)}$ is the usual monomial, with $w = (w(1),w(2))$. Another choice is to replace the negative coefficients in $P$ by their absolute value. In any event, you can try testing $Q^m P$ for large (whatever that means in this context) choices of $m$.

And a variation: if you suspect your thing is only nonnegative, add $\epsilon Q = \epsilon \cdot \sum_{w \in K} x^w$ to $P$, where $0 <\epsilon$ is much smaller than all the absolute values of the nonzero coefficients of $P$, and test $Q^m \cdot (P + \epsilon Q)$ for nonnegative coefficients. If succcessful, reduce $\epsilon$. If unsuccessful, increase $\epsilon$; this might give estimates of lower bounds.

The problem, as usual, is that it can't give a negative answer, although it might suggest where to look.