The first explicit example of nonnegative polynomial that is not a sum of squares was found by Motzkin in 1967. By the Arithmetic-geometric mean, $$\frac{x^4y^2 + x^2y^4 + 1}{3} \ge \sqrt{x^4y^2 \cdot x^2y^4 \cdot 1} = x^2y^2$$ hence $$x^4y^2 + x^2y^4 + 1 - 3x^2y^2 \ge 0.$$ The code belows construct the Motzkin polynomial using DynamicPolynomials.

In :
using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2

Out:
$$x^{4}y^{2} + x^{2}y^{4} - 3x^{2}y^{2} + 1$$

The Motzkin polynomial is nonnegative but is not a sum of squares as we can verify numerically as follows. We first need to pick an SDP solver, see here for a list of the available choices.

In :
using SumOfSquares, CSDP
factory = with_optimizer(CSDP.Optimizer);

In [ ]:
using SumOfSquares, MosekTools
factory = with_optimizer(Mosek.Optimizer);

In :
model = SOSModel(factory)
@constraint(model, motzkin >= 0) # We constraint motzkin to be a sum of squares

Out:
$(1)x^{4}y^{2} + (1)x^{2}y^{4} + (-3)x^{2}y^{2} + (1) \text{ is SOS}$
In :
optimize!(model)

CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 8.21e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  3.2410827e+01
Iter:  2 Ap: 9.49e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  3.0226979e+01
Iter:  3 Ap: 8.52e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  5.8043991e+00
Iter:  4 Ap: 4.43e-01 Pobj:  0.0000000e+00 Ad: 9.70e-01 Dobj: -2.0131592e+01
Declaring primal infeasibility.


We see that the problem is detected as infeasible...

In :
termination_status(model)

Out:
INFEASIBLE::TerminationStatusCode = 2

... and that the dual solution is a certificate of the infeasibility of the problem.

In :
dual_status(model)

Out:
INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4

Even if the Motzkin polynomial is not a sum of squares, it can still be certified to be nonnegative using sums of squares. Indeed a polynomial is certified to be nonnegative if it is equal to a fraction of sums of squares. The Motzkin polynomial is equal to a fraction of sums of squares whose denominator is $x^2 + y^2$. This can be verified numerically as follows:

In :
model = SOSModel(factory)
@constraint(model, (x^2 + y^2) * motzkin >= 0) # We constraint the (x^2 + y^2) * motzkin to be a sum of squares

Out:
$(1)x^{6}y^{2} + (2)x^{4}y^{4} + (1)x^{2}y^{6} + (-3)x^{4}y^{2} + (-3)x^{2}y^{4} + (1)x^{2} + (1)y^{2} \text{ is SOS}$
In :
optimize!(model)

CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 7.84e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  6.1906404e+01
Iter:  2 Ap: 9.45e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  6.0656534e+01
Iter:  3 Ap: 9.05e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  2.7707890e+01
Iter:  4 Ap: 7.75e-01 Pobj:  0.0000000e+00 Ad: 8.23e-01 Dobj:  5.8982728e+00
Iter:  5 Ap: 8.15e-01 Pobj:  0.0000000e+00 Ad: 9.46e-01 Dobj:  1.3960890e+00
Iter:  6 Ap: 7.61e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  4.0730107e-01
Iter:  7 Ap: 7.54e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  8.8693718e-02
Iter:  8 Ap: 7.91e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  4.7142486e-02
Iter:  9 Ap: 7.58e-01 Pobj:  0.0000000e+00 Ad: 9.68e-01 Dobj:  8.4582270e-03
Iter: 10 Ap: 7.89e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  4.0179083e-03
Iter: 11 Ap: 7.66e-01 Pobj:  0.0000000e+00 Ad: 9.92e-01 Dobj:  7.5288196e-04
Iter: 12 Ap: 7.96e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  3.7224432e-04
Iter: 13 Ap: 7.62e-01 Pobj:  0.0000000e+00 Ad: 9.73e-01 Dobj:  6.9932371e-05
Iter: 14 Ap: 6.52e-01 Pobj:  0.0000000e+00 Ad: 6.07e-01 Dobj:  2.9520639e-05


Now the problem is declared feasible by the solver...

In :
termination_status(model)

Out:
OPTIMAL::TerminationStatusCode = 1

... and the primal solution is a feasible point, hence it is a certificate of nonnegativity of the Motzkin polynomial.

In :
primal_status(model)

Out:
FEASIBLE_POINT::ResultStatusCode = 1

One may consider ourself lucky to have had the intuition that $x^2 + y^2$ would work as denominator. In fact, the search for the denominator can be carried out in parallel to the search of the numerator. In the example below, we search for a denominator with monomials of degrees from 0 to 2. If none is found, we can increase the maximum degree 2 to 4, 6, 8, ... This gives a hierarchy of programs to try in order to certify the nonnegativity of a polynomial by identifying it with a fraction of sum of squares polynomials. In the case of the Motzkin polynomial we now that degree 2 is enough since $x^2 + y^2$ works.

In :
model = SOSModel(factory)
X = monomials([x, y], 0:2)

Out:
6-element MonomialVector{true}:
x²
xy
y²
x
y
1 

We create a quadratic polynomial that is not necessarily a sum of squares since this is implied by the next constraint: deno >= 1.

In :
@variable(model, deno, Poly(X))

Out:
$$(noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname)$$

We want the denominator polynomial to be strictly positive, this prevents the trivial solution deno = 0 for instance.

In :
@constraint(model, deno >= 1)

Out:
$(noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname - 1) \text{ is SOS}$
In :
@constraint(model, deno * motzkin >= 0)

Out:
$(noname)x^{6}y^{2} + (noname)x^{5}y^{3} + (noname + noname)x^{4}y^{4} + (noname)x^{3}y^{5} + (noname)x^{2}y^{6} + (noname)x^{5}y^{2} + (noname)x^{4}y^{3} + (noname)x^{3}y^{4} + (noname)x^{2}y^{5} + (-3 noname + noname)x^{4}y^{2} + (-3 noname)x^{3}y^{3} + (-3 noname + noname)x^{2}y^{4} + (-3 noname)x^{3}y^{2} + (-3 noname)x^{2}y^{3} + (-3 noname)x^{2}y^{2} + (noname)x^{2} + (noname)xy + (noname)y^{2} + (noname)x + (noname)y + (noname) \text{ is SOS}$
In :
optimize!(model)

CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 7.95e-01 Pobj:  0.0000000e+00 Ad: 7.98e-01 Dobj: -1.4969375e+00
Iter:  2 Ap: 8.28e-01 Pobj:  0.0000000e+00 Ad: 7.89e-01 Dobj:  3.6682062e-01
Iter:  3 Ap: 8.25e-01 Pobj:  0.0000000e+00 Ad: 8.48e-01 Dobj: -4.2971015e-03
Iter:  4 Ap: 7.85e-01 Pobj:  0.0000000e+00 Ad: 8.48e-01 Dobj:  1.9519502e-04
Iter:  5 Ap: 7.02e-01 Pobj:  0.0000000e+00 Ad: 8.45e-01 Dobj: -1.7443054e-04
Iter:  6 Ap: 6.84e-01 Pobj:  0.0000000e+00 Ad: 7.45e-01 Dobj:  3.4354877e-05
Iter:  7 Ap: 6.59e-01 Pobj:  0.0000000e+00 Ad: 7.05e-01 Dobj: -1.2384593e-05
Iter:  8 Ap: 7.09e-01 Pobj:  0.0000000e+00 Ad: 5.96e-01 Dobj: -8.3434950e-06
Iter:  9 Ap: 6.92e-01 Pobj:  0.0000000e+00 Ad: 7.86e-01 Dobj:  1.3482976e-06
Iter: 10 Ap: 7.08e-01 Pobj:  0.0000000e+00 Ad: 6.19e-01 Dobj: -2.1547496e-07
Iter: 11 Ap: 6.29e-01 Pobj:  0.0000000e+00 Ad: 8.05e-01 Dobj:  1.1922145e-07
Iter: 12 Ap: 7.29e-01 Pobj:  0.0000000e+00 Ad: 6.16e-01 Dobj: -1.2524057e-08
Iter: 13 Ap: 8.49e-01 Pobj:  0.0000000e+00 Ad: 8.38e-01 Dobj:  8.7355197e-09
Iter: 14 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 7.32e-01 Dobj: -1.3938737e-09
Iter: 15 Ap: 8.46e-01 Pobj:  0.0000000e+00 Ad: 8.14e-01 Dobj:  5.3843217e-10
Iter: 16 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 7.60e-01 Dobj: -5.6655207e-11
Iter: 17 Ap: 9.05e-01 Pobj:  0.0000000e+00 Ad: 8.60e-01 Dobj:  1.1193794e-11
Iter: 18 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 7.25e-01 Dobj: -2.4333986e-12
Iter: 19 Ap: 6.67e-01 Pobj:  0.0000000e+00 Ad: 7.60e-01 Dobj:  9.9750431e-13
Iter: 20 Ap: 3.18e-02 Pobj:  0.0000000e+00 Ad: 5.86e-01 Dobj: -3.0148476e-13
Iter: 21 Ap: 8.80e-02 Pobj:  0.0000000e+00 Ad: 3.28e-01 Dobj:  3.1504853e-13
Iter: 22 Ap: 3.34e-03 Pobj:  0.0000000e+00 Ad: 3.27e-01 Dobj: -6.7360709e-14
Iter: 23 Ap: 1.21e-03 Pobj:  0.0000000e+00 Ad: 4.07e-01 Dobj:  1.4309873e-13
Iter: 24 Ap: 5.29e-04 Pobj:  0.0000000e+00 Ad: 1.98e-01 Dobj: -9.4979405e-14
Iter: 25 Ap: 2.24e-02 Pobj:  0.0000000e+00 Ad: 1.06e-01 Dobj:  1.0185252e-13
Iter: 26 Ap: 1.74e-02 Pobj:  0.0000000e+00 Ad: 1.20e-01 Dobj: -6.4880959e-14
Iter: 27 Ap: 2.36e-02 Pobj:  0.0000000e+00 Ad: 1.49e-01 Dobj:  8.2554804e-14
Iter: 28 Ap: 2.93e-03 Pobj:  0.0000000e+00 Ad: 1.65e-01 Dobj: -4.8411558e-14
Iter: 29 Ap: 8.78e-02 Pobj:  0.0000000e+00 Ad: 1.78e-01 Dobj:  4.5825457e-14
Iter: 30 Ap: 1.32e-03 Pobj:  0.0000000e+00 Ad: 3.17e-01 Dobj: -2.4389675e-14
Iter: 31 Ap: 4.83e-02 Pobj:  0.0000000e+00 Ad: 2.09e-01 Dobj:  3.0447835e-14
Iter: 32 Ap: 6.47e-03 Pobj:  0.0000000e+00 Ad: 2.76e-01 Dobj: -1.2564183e-14
Iter: 33 Ap: 2.78e-03 Pobj:  0.0000000e+00 Ad: 1.90e-01 Dobj:  4.5354455e-15
Iter: 34 Ap: 1.54e-05 Pobj:  0.0000000e+00 Ad: 7.77e-02 Dobj:  2.5022451e-16
Stuck at edge of primal feasibility, giving up.

In :
termination_status(model)

Out:
ALMOST_OPTIMAL::TerminationStatusCode = 7
In :
primal_status(model)

Out:
NEARLY_FEASIBLE_POINT::ResultStatusCode = 2

We can check the denominator found by the program using JuMP.getvalue

In :
value(deno)

Out:
$$11788.376578732263x^{2} + 11788.412205445507y^{2} + 20644.79368009344$$

Because a picture is worth a thousand words let's plot the beast. We can easily extend Plots by adding a recipe to plot bivariate polynomials.

In :
using RecipesBase
using MultivariatePolynomials
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
x, y, (x, y) -> p(variables(p) => [x, y])
end

In :
using Plots
plot(range(-2, stop=2, length=100), range(-2, stop=2, length=100), motzkin,
st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))

Out:
In [ ]: