Question

One way to generate integer solutions is as follows: Let $ v = (p + qi + rj + sk)$, where $p,q,r,s$ are rational integers and ${1,i,j,k}$ is the usual $\mathbb{R}$-basis for the algebra of quaternions. Let $v^{\prime}= (p -qi -rj - sk)$ and $w = vv^{\prime} = (p^{2}+q^{2}+r^{2}+s^{2}).$ If we expand $v^{2}$ in the form $x + yi + zj +tk$ for integers $x,y,z,t$, then we do have $x^{2} + y^{2} + z^{2} + t^{2} = w^{2}$. This is somehwat analogous to generating the Pythagorean triple $(x^{2}-y^{2})^{2}$ + $(2xy)^{2}$ = $(x^{2}+y^{2})^{2}$ by taking $({\rm Re}(x+iy)^{2})^{2} + ({\rm Im}(x+iy)^{2})^{2} = (x^{2}+y^{2})^{2}$. It may be better to work with Hurwitz quaternions for this question. Answer extended following the answer of Alex qubeat below, and slight rephrasing of the original question to place more emphasis on polynomials: At least in the context of polynomials in $\mathbb{R}[u]$ (it is easy to run out of letters in this game, so $u$ denotes an indeterminate here), more solutions may be manufactured using the fact that the solutions to $x(u)^{2} +y(u)^{2} + z(u)^{2} + t(u)^{2} = w(u)^{2}$ with $x(u),y(u),z(u),t(u),w(u) \in \mathbb{R}[u]$ have a semigroup structure. Alex's answer combines a fixed solution with an essentially ``constant" solution, to produce other solutions, but solutions can be combined in other ways. Let $\mathbb{H}$ denote the algebra of real quaternions. Then the map $N: \mathbb{H}[u] \to \mathbb{R}[u]$ with $N(x(u) +iy(u)+jz(u) + kt(u)) =x(u)^{2} +y(u)^{2} + z(u)^{2} +t(u)^{2}$ (for real polynomials $x(u),y(u),z(u),t(u))$ is multiplicative. The polynomials $p(u) \in \mathbb{H}[u]$ such that $N(p(u))$ is a non-zero square in $\mathbb{R}[u]$ are closed under multiplication. As in the integral case, one way to ensure that $N(x(u) + iy(u) + jz(u) + kt(u))$ is a square is to take $x(u) + iy(u) +jz(u)+kt(u)$ to be of the form $(a(u)+ib(u)+jc(u)+kd(u))^{2}$ for real polynomials $a(u),b(u),c(u),d(u)$, but it is also clear from this discussion that $N(p_{1}(u)p_{2}(u) \ldots p_{n}(u))$ is a square in $\mathbb{R}[u]$ as long as any given $p_{i}(u)\in \mathbb{H}[u]$ occurs an even number of times in the (non-commuting) product. In fact, it is permissible to count the total number of occurences of $\mathbb{H}$-conjugates of any given $p_{i}(u)$, that is elements of the form $v^{-1}p_{i}(u)v$, where $v$ is a non-zero element of $\mathbb{H}$, since $N$ takes a constant value on the $\mathbb{H}^{\times}$-conjugacy class of $p_{i}(u)$.

Further remarks. May 6: The polynomial ring $\mathbb{H}[u]$ has a natural division ring of fractions, $\mathbb{H}(u)$, which is isomorphic to a certain ring of $2 \times 2$ matrices over $\mathbb{C}(u)$ (I can't get the latex right for matrices, but the ring should be clear- in particular, the determinants of elements in the ring are elements of $\mathbb{R}(u)$). There is a natural ring homomorphism $\sigma$ from $\mathbb{H}(u)$ to ${\rm SO}(3,\mathbb{R}(u))$, obtained by letting $\mathbb{H}(u) \backslash {0}$ act by conjugation of the $\mathbb{R}(u)$-span of ${i,j,k}.$ The question about polynomials amounts to determining $S\sigma$, where $S$ is the set of elements of $\mathbb{H}(u)$ whose determinants are squares in $\mathbb{R}(u)$. This is the same as $T\sigma$, where $T$ is the set of elements of determinant $1$, since the non-zero elements of $\mathbb{R}(u)$ acts trivially by conjugation on $\mathbb{H}(u)$. Now it is clear that $T\sigma$ is a normal subgroup of ${\rm SO}(3,\mathbb{R})$, and that the factor group is an elementary Abelian $2$-group (that is, an Abelian group of exponent $2$). I do not see at present how to calculate the size of this group, but record this line of thinking in case anyone else can exploit it further.

More remarks, May 13: I think that $\mathbb{H}(u)$ has an involutory automorphism $\sigma$ which fixes $\mathbb{R}(u)$ elementwise, and which has $i \sigma = -i$,$j\sigma = k$ and $k\sigma =j.$ Then $N(p(u)\sigma) = N(p(u))$ for all $p(u) \in \mathbb{H}(u)$. Hence, for any $p(u) \in \mathbb{H}(u)$, the element $q(u) = p(u).p(u)\sigma$ will have $N(q(u)) =w(u)^{2}$ for some $w(u) \in \mathbb{R}(u)$, so this will yield a new type of solution (at least in the context of this answer).