Let
$$ \mathbf{u}=\sigma_x (a+bi)+\sigma_y (c+di)+\sigma_z (e+hi) $$
Squaring the vector (geometric product, Clifford algebra), we get:
$$ \begin{align} \mathbf{u}^2&=(a+bi)^2+(c+di)^2+(e+hi)^2\\ &=a^2-b^2+2iab+c^2-d^2+2icd+e^2-h^2+2ieh\\ &=(a^2-b^2+c^2-d^2+e^2-h^2)+2i(ab+cd+eh) \end{align} $$
I am now interested in a linear transformation $T$ which leaves $\mathbf{u}^2$ invariant.
$$ (T\mathbf{u})^2=\mathbf{u}^2 $$
Since $\mathbf{u}$ is a complex number, then we get two equations that must independently remain invariant by the transformation:
$$ \begin{align} k_1 &:=(a^2-b^2+c^2-d^2+e^2-h^2)\\ k_2 &:=2i(ab+cd+eh) \end{align} $$
I have previously asked the question [What symmetry corresponds to the square of a complex number, as opposed to the norm? ] for a single complex number $z=a+bi$, and $z^2=a^2-b^2+2iab$. I was told in the comments that the real part represents the set of Lorentz transformations in the plane, but with the imaginary constraint, it restricts it to mere reflexions O(1).
So with the added freedom of 2 or 3 orthogonal complex numbers, what is $T$?
If I had to take a guess, I would say its related to the $U(1)$ group in each plane: $\sigma_x,\sigma_y,\sigma_z$
Your expression for $\mathbf u^2$ is the product of complex matrices $u^Tu$, with $$u=\begin{bmatrix}u_1\\ u_2\\ u_3\end{bmatrix}=\begin{bmatrix}a+bi\\ c+di\\ e+hi\end{bmatrix}$$ $$u^T=\begin{bmatrix}u_1&u_2&u_3\end{bmatrix}.$$ To say that a transformation matrix $A$ preserves $u^Tu$, $$u^Tu=(Au)^T(Au)=u^TA^TAu$$ (for all $u$), is to say that $A^TA=I$. The set of all such matrices is the complex orthogonal group.
($A^T$ is the ordinary transpose, not the conjugate-transpose as in the unitary group.)