How to get the left and right rotation quaternions between two $\mathbb{R}^4$ unit vectors? This question is somewhat similar to this, but in higher dimension.
I have two non-collinear unit vectors in $\mathbb{R}^4$, $\bf a$ and $\bf b$. I know from Wikipedia that there is at least one pair of quaternions $Q_L$ and $Q_R$ that is able to rotate $\bf a$ to $\bf b$, i.e.:
$$
{\bf b} = Q_L {\bf a} Q_R
$$
(where multiplication is the Hamilton product)
Questions:


*

*What is an easy to understand and fast to compute way to find $Q_L$ and $Q_R$?

*What are the rotation quaternions that moves $\bf a$ towards $\bf b$, but by a given angle $\theta$ ?


Edit: the comments made me realize I need further constraints to get the rotation I need. So, $\bf a$ and $\bf b$ spans a subspace. The rotation I need is closed in this subspace. I think this restriction reduces the number of possible solutions from infinite to two, (as there are two pairs of quaternions that describes the same rotation).
 A: Let $a,b$ be two unit quaternions, for simplicity $a,b\ne\pm1$ and $b\ne\pm a$.
There are many pairs $(p,q)$ of quaternions for which $b=paq$, e.g. $(ba^{-1},1)$ or $(1,a^{-1}b)$.
You can even parametrize them all by solving for $p$ in terms of $q$ or vice-versa.
There is a unique pair (up to $\pm$) for which $F(x)=pxq$ acts as a 2D rotation in the plane spanned by $a,b$ but fixes all points in the plane orthogonal to it, and that is $(\sqrt{ba^{-1}},\sqrt{a^{-1}b})$.
The square roots are defined by halving the convex angle of polar form, i.e. $\sqrt{e^{\theta\mathbf{u}}}=e^{(\theta/2)\mathbf{u}}$ if $\theta\in[0,\pi)$; in practice you can calculate this using half-angle formulas on real/imag parts:
$$ \sqrt{r+\mathbf{v}}=\sqrt{\frac{1+r}{2}}+\sqrt{\frac{1-r}{2}}\frac{\mathbf{v}}{\|\mathbf{v}\|}. $$
(I define quaternions as sums of scalars $r$ and 3D vectors $\mathbf{v}$; these are real/imaginary parts.)

Here's how to get the formula. First define the left/right multiplication maps $L_p(x)=px$ and $R_p(x)=xp$. If we can create a rotation $G(x)$ which sends $1$ to $a^{-1}b$ and is a nonrotation in the complent, then we can define $F=L_a\circ G\circ L_a^{-1}$. The effect of the $L_a^{-1}$ is to send the $ab$-plane to the $1(a^{-1}b)$-plane, and the former's complement to the latter's complement. After we rotate purely in the $1(a^{-1}b)$-plane with $G$ we go back to the $ab$-plane with $L_a$.
If $p=e^{\theta\mathbf{u}}$ then $L_p$ and $R_p$ both act as a rotation by $\theta$ in the $1\mathbf{u}$-plane and its complement, the only difference being opposite directions in the complement (this follows algebraically from the fact orthogonal vectors anticommute). Therefore if we were to compose $L_{\sqrt{p}}$ and $R_{\sqrt{p}}$ we get a rotation from $1$ to $p$ in the $1\mathbf{u}$-plane and a nonrotation in the complement.
Thus $F(x)=a\big(\sqrt{a^{-1}b}(a^{-1}x)\sqrt{a^{-1}b}\big)$. Let $q=a\sqrt{a^{-1}b}a^{-1}$. Then $q^2=ba^{-1}$ so $q=\pm\sqrt{ba^{-1}}$ (square roots of nonreal quaternions are unique up to $\pm$; you can prove this with polar form). Since $\sqrt{a^{-1}b}$ has positive real part, so does its conjugate $q$, so $q$ must be $+\sqrt{ba^{-1}}$.
Conclude $F(x)=\sqrt{ba^{-1}}x\sqrt{a^{-1}b}$.

A more symmetric version of the above is how I intuited the formula: $L_{\sqrt{ba^{-1}}}$ and $R_{\sqrt{ab^{-1}}}$ both rotate $a$ halfway to $b$, but rotate in opposite directions in the complement.
