Find three polynomials whose squares sum up to $x^4 + y^4 + x^2 + y^2$

Prove that $$p(x,y) = x^4 + y^4 + x^2 + y^2$$ can be written as a sum of squares of three polynomials over $x,y$ for real numbers.

$(\sqrt{2\sqrt{2}-2}y^2+x)^2 + (-\sqrt{2\sqrt{2}-2}xy + y)^2 + (x^2 + (1-\sqrt{2})y^2)^2 = x^4+y^4+x^2+y^2$.

In fact, there are many other solutions by using brute force, as suggested by Will Jagy in the comments.

Step 1 Let's suppose that $x^4+y^4+x^2+y^2 = \sum_{i=1}^3 (a_i x^2 + b_i xy + c_i y^2 + d_i x + e_i y)^2$. If we write $\vec{a} = (a_1,a_2,a_3)$ and so forth, we obtain the following identities by comparing coefficients:

Part (a): $x^2: \| \vec{d} \|^2 = 1$, $y^2: \| \vec{e} \|^2 = 1$, $xy: \langle \vec{d}, \vec{e} \rangle = 0$.

Part (b): $x^3: \langle \vec{a}, \vec{d} \rangle = 0$, $y^3: \langle \vec{c}, \vec{e} \rangle = 0$, $x^2y: \langle \vec{b}, \vec{d} \rangle = - \langle \vec{a}, \vec{e} \rangle$, $xy^2: \langle \vec{b} ,\vec{e} \rangle = - \langle \vec{c}, \vec{d} \rangle$.

Part (c): $x^4: \| \vec{a} \|^2 = 1$, $y^4: \| \vec{c} \|^2 = 1$, $x^2y^2: \| \vec{b} \|^2 + 2 \langle \vec{a}, \vec{c} \rangle = 0$, $x^3y: \langle \vec{b}, \vec{a} \rangle = 0$, $xy^3: \langle \vec{b}, \vec{c} \rangle = 0$.

(All inner products are standard Euclidean inner products)

Step 2 We consider an orthonormal basis of $\mathbb{R}^3$: $\vec{d}, \vec{e}, \vec{f}$. This used part (a).

Using part (b), we get

$\vec{a} = x_1 \vec{e} + y_1 \vec{f}$,

$\vec{b} = -x_1 \vec{d} - x_2 \vec{e} + y_3 \vec{f}$,

$\vec{c} = x_2 \vec{d} + y_2 \vec{f}$

for some reals $x_1,x_2,y_1,y_2,y_3$.

Step 3 Part (c) tells us that

$x_1^2 + y_1^2 = 1$, $x_2^2 + y_2^2 = 1$. (From norm of $\vec{a}, \vec{c}$ being 1.)

$y_1y_3 = x_1x_2$, $y_2y_3 = x_1x_2$ (From $\vec{b}$ perpendicular to $\vec{a}, \vec{c}$.)

$x_1^2+x_2^2+y_3^2 + 2y_1y_2 = 0$.

In the second equation, if $y_3$ is not 0, then $y_1 = y_2$, which forces $x_1 = x_2 = y_1 = y_2 = y_3 = 0$ in the third equation, contradicting the first equation.

So $y_3 = 0$, which implies $x_1x_2 = 0$. WLOG assume $x_1 = 0$, then $y_1 = \pm 1$. The third equation becomes

$0 = x_2^2 + 2y_1y_2 = 1 - y_2^2 + 2y_1y_2 = 2 - (y_2 - y_1)^2$.

So $y_2 - y_1 = \pm \sqrt{2}$. The first equation tells us that $y_1,y_2$ has norms less than 1, which implies that $y_2 = \pm (1 - \sqrt{2})$. This gives $x_2 = \pm \sqrt{2\sqrt{2} - 2}$.

To summarize, up to symmetry of $(x_1,y_1) \leftrightarrow (x_2,y_2)$, we have $x_1 = 0$, $y_1 = \pm 1$, $y_2 = \pm (1 - \sqrt{2})$ (same sign as $y_1$), $x_2 = \pm \sqrt{2\sqrt{2} - 2}$ and $y_3 = 0$.

Step 4 By picking arbitrary orthonormal basis $\vec{d}, \vec{e}, \vec{f}$, we can substitute the values for $x_1,x_2,y_1,y_2,y_3$ above to get solutions for $\vec{a}, \vec{b}, \vec{c}$. All these would give the required sum of squares decomposition. The decomposition I took above comes from choosing the standard orthonormal basis of $\mathbb{R}^3$.