Is there a reference for expressing biquaternion Lorentz transformations as a matrix? Background: In 3D geometry, a suitable subset of the quaternions can be used to model rotations, namely those where $qq^* = 1$. There is a well-known direct correspondence between the the coefficients of such a quaternion $q = w + x\mathbf i + y\mathbf j + z\mathbf k$ and the coefficients of a rotation matrix $M$:
$M = \begin{bmatrix}
1 - 2(y^2 + z^2) & 2(xy - zw) & 2(xz + yw) \\
2(xy + zw) & 1 - 2(x^2 + z^2) & 2(yz - xw) \\
2(xz - yw) & 2(yz + xw) & 1 - 2(x^2 + y^2)
\end{bmatrix}$
In Minkowski space, the biquaternions can extend this to model Lorentz transformations. If we use $q^*$ for quaternion conjugation and $\overline q$ for complex conjugation, then $qg\overline q^*$ (which is equal to $qg\overline{q^*}$) will be a Lorentz transformation when $qq^* = 1$. This naturally includes the ordinary rotations SO(3) as a special case, but also Lorentz boosts etc.
The $g$ that is being transformed is a so-called minquat, a restricted subset of the biquaternions with imaginary vector components and a real scalar component. (We could multiply the coefficients by $-i$, so that the vector components are real and we have a direct extension of the vectors used with quaternion rotation, but then the scalar component is negated which is awkward.) The action of the Lorentz transformations can be shown to be closed under this space. The real scalar component corresponds to time when dealing with relativity, or simply the 4th coordinate when dealing with the Hyperboloid model of hyperbolic space, which is my interest here.
My specific question is: is there a reference for the 4x4 transformation matrix derived from the complex biquaternion coefficients, analagous to the 3x3 matrix above? I'm working it out myself by hand, but it would be nice to check my work.
 A: I've pored through https://en.wikipedia.org/wiki/History_of_Lorentz_transformations, but surprisingly this formula does not seem to be directly listed there or anywhere else obvious. But after working it out by hand, I am fairly confident that given biquaternion $q=w+x\mathbf i+y\mathbf j+z\mathbf k$, the equivalent transformation matrix will be $M=$
\begin{bmatrix}
|w|^2+|x|^2+|y|^2+|z|^2 &&
2\operatorname{Im}(x\overline w+y\overline z) &&
2\operatorname{Im}(y\overline w+z\overline x) &&
2\operatorname{Im}(z\overline w+x\overline y) \\
2\operatorname{Im}(x\overline w-y\overline z) &&
|w|^2+|x|^2-|y|^2-|z|^2 &&
2\operatorname{Re}(x\overline y-z\overline w) &&
2\operatorname{Re}(z\overline x+y\overline w) \\
2\operatorname{Im}(y\overline w-z\overline x) &&
2\operatorname{Re}(x\overline y+z\overline w) &&
|w|^2-|x|^2+|y|^2-|z|^2 &&
2\operatorname{Re}(y\overline z-x\overline w) \\
2\operatorname{Im}(z\overline w-x\overline y) &&
2\operatorname{Re}(z\overline x-y\overline w) &&
2\operatorname{Re}(y\overline z+x\overline w) &&
|w|^2-|x|^2-|y|^2+|z|^2 \\
\end{bmatrix}
This reduces to $\begin{bmatrix}1 && 0 \\ 0 && R \end{bmatrix}$, where $R$ is a rotation matrix equivalent to the form given in the question, when the coefficients are real - so that's a promising sign. I've also verified that the first two columns are orthogonal, as well as columns 2+3. Given the nature of the matrix, I'm confident this orthogonality extends to the rest of the pairs. The verification/calculation of orthogonality is rather long and tedious, so I'm omitting it here.
It's also easy to verify that Lorentz boosts work properly. I'm adopting the biquaternion convention of $\sqrt{-1} = \mathbf h$, to avoid confusion with the quaternion unit $\mathbf i$. Given $q = \cosh(\alpha/2)+\mathbf{h}\sinh(\alpha/2)\mathbf i$, the resulting transformation will be
\begin{bmatrix}
\cosh(\alpha/2)^2+\sinh(\alpha/2)^2 && 2\sinh(\alpha/2)\cosh(\alpha/2) && 0 && 0 \\
2\sinh(\alpha/2)\cosh(\alpha/2) && \cosh(\alpha/2)^2+\sinh(\alpha/2)^2 && 0 && 0 \\
0 && 0 && 1 && 0 \\
0 && 0 && 0 && 1 \\
\end{bmatrix}
Using the double-angle identities:
\begin{bmatrix}
\cosh(\alpha) && \sinh(\alpha) && 0 && 0 \\
\sinh(\alpha) && \cosh(\alpha) && 0 && 0 \\
0 && 0 && 1 && 0 \\
0 && 0 && 0 && 1 \\
\end{bmatrix}
Exactly as expected for an X-axis boost. (And exactly analogously to the way rotations work, except with hyperbolic functions.)
Finally, here's the derivation of the formula. I'm following the strategy and notation used at https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion-derived_rotation_matrix, with the additional wrinkle that "$\cdot$" means the "ordinary" (symmetric) dot-product, not the anti-symmetric inner-product that is often used with complex vector algebra. This makes the calculations simpler and maintains the correspondence to the quaternion formulas.
Let $q = (w, \mathbf v) = w + x\mathbf i + y\mathbf j + z \mathbf k$. Then for arbitrary minquats $g = (p_w, \mathbf h \mathbf p)$ where $p_w$ and $\mathbf p$ are real, we want to  find the result of $qg\overline q^*$. I.e.
\begin{align}({p_w}', \mathbf h \mathbf p')
=\ (&(w,\mathbf v)(p_w, \mathbf h \mathbf p))(\overline w, -\mathbf{\overline v})\\
=\ (&w p_w - \mathbf h \mathbf v\cdot\mathbf p, \mathbf h w\mathbf p + \mathbf v p_w + \mathbf h\mathbf v \times \mathbf p)(\overline w, -\mathbf{\overline v})\\
=\ (&w \overline w p_w - \mathbf h \overline w \mathbf v\cdot\mathbf p + \mathbf h w \mathbf p\cdot\mathbf{\overline v} + (\mathbf v \cdot \mathbf{\overline v}) p_w + \mathbf h\mathbf{\overline v} \cdot(\mathbf v \times \mathbf p),\\
& \mathbf h w \overline w \mathbf p + \overline w \mathbf v p_w + \mathbf h \overline w \mathbf v \times \mathbf p - w \mathbf{\overline v} p_w + \mathbf h \mathbf{\overline v}(\mathbf v\cdot\mathbf p) - (\mathbf h w\mathbf p + \mathbf v p_w + \mathbf h\mathbf v \times \mathbf p) \times \mathbf{\overline v})\\
=\ (&|w|^2 p_w + \mathbf h(w \mathbf{\overline v}-\overline w \mathbf v)\cdot\mathbf p + |\mathbf v|^2 p_w + \mathbf h \mathbf p \cdot (\mathbf{\overline v} \times \mathbf v),\\
& \mathbf h |w|^2 \mathbf p + (\overline w \mathbf v - w \mathbf{\overline v}) p_w + \mathbf h \overline w \mathbf v \times \mathbf p + \mathbf h w \mathbf{\overline v} \times \mathbf p + \mathbf h \mathbf{\overline v}(\mathbf v\cdot\mathbf p) + \mathbf{\overline v} \times \mathbf v p_w + \mathbf h\mathbf{\overline v} \times (\mathbf v \times \mathbf p))\\
=\ (&(|w|^2 + |\mathbf v|^2)p_w + \mathbf h \cdot 2 \mathbf h \operatorname{Im}(w \mathbf{\overline v})\cdot\mathbf p + \mathbf h \mathbf p \cdot (\mathbf{\overline v} \times \mathbf v),\\
& \mathbf h |w|^2 \mathbf p + \mathbf h 2\operatorname{Im}(\overline w \mathbf v) p_w + \mathbf h 2 \operatorname{Re}(\overline w \mathbf v) \times \mathbf p + \mathbf h \mathbf{\overline v}(\mathbf v\cdot\mathbf p) - (\mathbf h \cdot \mathbf h) \mathbf{\overline v} \times \mathbf v p_w + \mathbf h\mathbf{\overline v} \times (\mathbf v \times \mathbf p))\\
=\ (&(|w|^2 + |\mathbf v|^2)p_w + 2\operatorname{Im}(\overline w \mathbf v)\cdot\mathbf p + \mathbf h (\mathbf{\overline v} \times \mathbf v) \cdot \mathbf p,\\
& \mathbf h \left(|w|^2 \mathbf p + 2\operatorname{Im}(\overline w \mathbf v) p_w - \mathbf h \mathbf{\overline v} \times \mathbf v p_w + 2 \operatorname{Re}(\overline w \mathbf v) \times \mathbf p + \mathbf{\overline v}(\mathbf v\cdot\mathbf p) + \mathbf{\overline v} \times (\mathbf v \times \mathbf p)\right))\\
=\ (&(|w|^2 + |\mathbf v|^2)p_w + (2\operatorname{Im}(\overline w \mathbf v) + \mathbf h (\mathbf{\overline v} \times \mathbf v)) \mathbf p,\\
& \mathbf h \left( \left(2\operatorname{Im}(\overline w \mathbf v) - \mathbf h (\mathbf{\overline v} \times \mathbf v) \right) p_w + \left(|w|^2 \mathbf I + 2 [\operatorname{Re}(\overline w \mathbf v)]_\times + \mathbf{\overline v} \otimes \mathbf v + [\mathbf{\overline v}]_\times [\mathbf v]_\times \right) \mathbf p \right))\\
\end{align}
Then we need to compute the (sub)matrices that show up in that result.
$[\mathbf v]_\times = \begin{bmatrix}
0 && -z && y \\
z && 0 && -x \\
-y && x && 0 \\
\end{bmatrix}, [\mathbf{\overline v}]_\times = \begin{bmatrix}
0 && -\overline z && \overline y \\
\overline z && 0 && -\overline x \\
-\overline y && \overline x && 0 \\
\end{bmatrix},\\
[\mathbf{\overline v}]_\times [\mathbf v]_\times = \begin{bmatrix}
-|z|^2-|y|^2 && x \overline y && \overline z x\\
\overline x y && -|x|^2-|z|^2 && y \overline z \\
z \overline x && \overline y z && -|y|^2-|x|^2 \\
\end{bmatrix},\\
\mathbf{\overline v} \times \mathbf v = [\mathbf{\overline v}]_\times \mathbf v = \begin{bmatrix}
\overline y z - y \overline z \\
\overline z x - z \overline x \\
\overline x y - x \overline y \\
\end{bmatrix} = -2\mathbf h\begin{bmatrix}
\operatorname{Im}(y \overline z) \\
\operatorname{Im}(z \overline x) \\
\operatorname{Im}(x \overline y) \\
\end{bmatrix},\\
2 [\operatorname{Re}(\overline w \mathbf v)]_\times =
\begin{bmatrix}
0 && -2 \operatorname{Re}(z \overline w) && 2 \operatorname{Re}(y \overline w) \\
2 \operatorname{Re}(z \overline w) && 0 && -2 \operatorname{Re}(x \overline w) \\
-2 \operatorname{Re}(y \overline w) && 2 \operatorname{Re}(x \overline w) && 0 \\
\end{bmatrix},\\
\mathbf{\overline v} \otimes \mathbf v = \begin{bmatrix}
|x|^2 && \overline x y && z \overline x \\
x \overline y && |y|^2 && \overline y z \\
\overline z x && y \overline z && |z|^2 \\
\end{bmatrix}$
Shoving everything back in to the last equation gives us
$({p_w}', \mathbf h \mathbf p') = \left( (|w|^2 + |\mathbf v|^2)p_w + \left(2\begin{bmatrix}\operatorname{Im}(x \overline w) \\
\operatorname{Im}(y \overline w) \\
\operatorname{Im}(z \overline w)\end{bmatrix}^T +
\mathbf h \cdot -2 \mathbf h \begin{bmatrix}\operatorname{Im}(y \overline z) \\
\operatorname{Im}(z \overline x) \\
\operatorname{Im}(x \overline y)\end{bmatrix}^T
\right) \mathbf p,\\
 \mathbf h \left( \left(2\begin{bmatrix}\operatorname{Im}(x \overline w) \\
\operatorname{Im}(y \overline w) \\
\operatorname{Im}(z \overline w)\end{bmatrix} - \mathbf h \cdot -2 \mathbf h \begin{bmatrix}\operatorname{Im}(y \overline z) \\
\operatorname{Im}(z \overline x) \\
\operatorname{Im}(x \overline y)\end{bmatrix} \right) p_w + \left(|w|^2 \mathbf I + \begin{bmatrix}
0 && -2 \operatorname{Re}(z \overline w) && 2 \operatorname{Re}(y \overline w) \\
2 \operatorname{Re}(z \overline w) && 0 && -2 \operatorname{Re}(x \overline w) \\
-2 \operatorname{Re}(y \overline w) && 2 \operatorname{Re}(x \overline w) && 0 \\
\end{bmatrix} +\\
\begin{bmatrix}
|x|^2 && \overline x y && z \overline x \\
x \overline y && |y|^2 && \overline y z \\
\overline z x && y \overline z && |z|^2 \\
\end{bmatrix} + \begin{bmatrix}
-|z|^2-|y|^2 && x \overline y && \overline z x\\
\overline x y && -|x|^2-|z|^2 && y \overline z \\
z \overline x && \overline y z && -|y|^2-|x|^2 \\
\end{bmatrix} \right) \mathbf p \right) \right)$
which simplifies to
$\left( (|w|^2 + |\mathbf v|^2)p_w + 2\begin{bmatrix}
\operatorname{Im}(x \overline w + y \overline z)\\
\operatorname{Im}(y \overline w + z \overline x)\\
\operatorname{Im}(z \overline w + x \overline y)\end{bmatrix}^T \mathbf p,\\
 \mathbf h \left( 2\begin{bmatrix}
\operatorname{Im}(x \overline w - y \overline z)\\
\operatorname{Im}(y \overline w - z \overline x)\\
\operatorname{Im}(z \overline w - x \overline y)\end{bmatrix} p_w +
 \left(\begin{bmatrix}
|w|^2 && -2 \operatorname{Re}(z \overline w) && 2 \operatorname{Re}(y \overline w) \\
2 \operatorname{Re}(z \overline w) && |w|^2 && -2 \operatorname{Re}(x \overline w) \\
-2 \operatorname{Re}(y \overline w) && 2 \operatorname{Re}(x \overline w) && |w|^2 \\
\end{bmatrix} +\\
\begin{bmatrix}
|x|^2 -|y|^2-|z|^2 && x \overline y + \overline x y && z \overline x + \overline z x \\
x \overline y + \overline x y && -|x|^2+|y|^2-|z|^2 &&  y \overline z + \overline y z \\
z \overline x + \overline z x && y \overline z + \overline y z&& -|x|^2-|y|^2 +|z|^2 \\
\end{bmatrix} \right) \mathbf p \right) \right)$
From there the final simplifications and gathering it into the 4x4 matrix form is straightforward.
