Rotations preserve cross products?

In relation to this question here, I'm trying to prove that $\bf R(a \times b) = R(a) \times R(b)$ for any given rotation $\bf R$ in $\mathbb{R}^3$.

Edit: If one can prove that this holds for the unit vectors (e.g. $\bf R(\hat j) \times R(\hat k) = R(\hat i)$) then that's all there is to it, given the distributivity of cross product and linearity of $\bf R$.

Any rotation in $\Bbb{R}^3$ can be represented by $R_x(\alpha)R_y(\beta)R_y(\gamma)$, where $$\begin{alignat}{1} R_x(\theta) &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta \\[3pt] 0 & \sin \theta & \cos \theta \\[3pt] \end{bmatrix} \\[6pt] R_y(\theta) &= \begin{bmatrix} \cos \theta & 0 & \sin \theta \\[3pt] 0 & 1 & 0 \\[3pt] -\sin \theta & 0 & \cos \theta \\ \end{bmatrix} \\[6pt] R_z(\theta) &= \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\[3pt] \sin \theta & \cos \theta & 0\\[3pt] 0 & 0 & 1\\ \end{bmatrix}. \end{alignat} $$ If we show that a rotation around one axis preserves the cross product then we know that the rotation around the others preserves it as well because the arithemetic is analogous and we can apply, for example, $R_x$ to all elements in the equality $$R_y(\beta)R_y(\gamma)(a \times b) = R_y(\beta)R_y(\gamma)a \times R_y(\beta)R_y(\gamma) b$$

So we have $$\begin{bmatrix} \cos \theta & -\sin \theta & 0 \\[3pt] \sin \theta & \cos \theta & 0\\[3pt] 0 & 0 & 1\\ \end{bmatrix}\begin{bmatrix} a_1 \\[3pt] a_2\\[3pt] a_3\\ \end{bmatrix} \times \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\[3pt] \sin \theta & \cos \theta & 0\\[3pt] 0 & 0 & 1\\ \end{bmatrix}\begin{bmatrix} b_1 \\[3pt] b_2\\[3pt] b_3\\ \end{bmatrix}=\\ \begin{bmatrix} -a_3b_2\cos\theta +a_2b_3\cos\theta - a_3b_1\sin\theta + a_1b_3\sin\theta\\[3pt] a_3b_1\cos\theta -a_1b_3\cos\theta - a_3b_2\sin\theta + a_2b_3\sin\theta\\[3pt] -a2b1\cos^2\theta + a_1b_2\cos^2\theta-a_2b_1\sin^2\theta +a_1b_2\sin^2\theta\\ \end{bmatrix}=\\ \begin{bmatrix} (-a_3b_2+a_2b_3)\cos\theta+(-a_3b_1+a_1b_3)\sin\theta\\[3pt] (a_3b_1-a_1b_3)\cos\theta+(-a_3b_2+a_2b_3)\sin\theta\\[3pt] a_1b_2-a_2b_1\\[3pt] \end{bmatrix} $$ Which you will see is equal to $$ \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\[3pt] \sin \theta & \cos \theta & 0\\[3pt] 0 & 0 & 1\\ \end{bmatrix}\left(\begin{bmatrix} a_1 \\[3pt] a_2\\[3pt] a_3\\ \end{bmatrix}\times\begin{bmatrix} b_1 \\[3pt] b_2\\[3pt] b_3\\ \end{bmatrix}\right) $$

We can think of $\Bbb{R}^3$ as the space of imaginary quaternions $\{ai+bj+ck|i,j,k \in \Bbb{R}\}$. Then the following facts are standard (the first one is in the linked wikipedia article, and the other two are easy to verify given a quaternion multiplication table):

• Any rotation on $\Bbb{R}^3$ may be realized as the conjugate action of some unit quaternion: that is, given any rotation $\bf{R}$, there is some unit quaternion $q({\bf R})$ such that ${\bf R}x=q({\bf R})x q({\bf R})^{-1}$.
• Conjugating by a unit quaternion commutes with taking its imaginary part: if $q$ is a unit quaternion and $r$ an arbitrary quaternion, then $$ q\Im(r)q^{-1}=\Im(qrq^{-1}) $$
• If $x,y \in \Bbb{R}^3$, their quaternion product is given in terms of the ordinary vector dot product and cross product: $$ xy=-x \cdot y+x \times y \, ; $$ in particular, $x \times y = \Im(xy)$.

If you believe all three of these facts, the proof you want is immediate. Setting $q=q({\bf R})$, we have

\begin{eqnarray} {\bf R}(x \times y) &=& q(x \times y)q^{-1}\\ & = &q\Im(xy)q^{-1}\\ &=&\Im(qxyq^{-1})\\ &=&\Im((qxq^{-1})(qyq^{-1}))\\ &=&qxq^{-1} \times qyq^{-1} \\ &=&{\bf R}x \times {\bf R}y \end{eqnarray} as desired.