So the hyperbola-preserving transformation in 2 dimensional space is given by the matrix \begin{pmatrix} \cosh(\phi) & \sinh(\phi) \\ \sinh(\phi) & \cosh(\phi) \end{pmatrix}
I'm wondering what such a matrix would be in 3 dimensional space (so that it preserves 2 dimensional hyperboloids) and 4 dimensional space (so that it preserves 3 dimensional hyperboloids). Sources or derivations would be appreciated. Thank you!
In a way, your transformation matrix is a variation of a common 2d rotation matrix
$$\begin{pmatrix}\cos\phi&-\sin\phi\\\sin\phi&\cos\phi\end{pmatrix}\;.$$
Where the above preserves the unit circle $x^2+y^2=1$, yours preserves the hyperbola $x^2-y^2=1$. The unit circle here corresponds to the unit sphere in 3d. There are many ways to describe 3d rotations, but one very common one is to describe them as a product of rotations around the coordinate axes. You can do the same for your hyperboloid as well.
For example, the one-sheeted hyperboloid $x^2+y^2-z^2=1$ has rotational symmetry around the $z$ axis. So you'd have these three “rotation” matrices:
$$ \begin{pmatrix} \cos\alpha & -\sin\alpha & 0 \\ \sin\alpha & \cos\alpha & 0 \\ 0 & 0 & 1 \end{pmatrix} \qquad \begin{pmatrix} \cosh\beta & 0 & \sinh\beta \\ 0 & 1 & 0 \\ \sinh\beta & 0 & \cosh\beta \end{pmatrix} \qquad \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cosh\gamma & \sinh\gamma \\ 0 & \sinh\gamma & \cosh\gamma \end{pmatrix} $$
Each of them preserves the hyperboloid, so a product of them will preserve it as well. The two-sheeted hyperboloid $z^2-y^2-x^2=1$ is preserved by the above matrices, too. If you want $x^2-y^2-z^2=1$ instead, you have to change coordinates, so that the rotation around $x$ becomes a regular rotation while the other two use hyperbolic functions.
I don't know whether it would make sense to translate any of the other rotation formalisms (like axis-angle or quaternions) into something preserving a hyperboloid. Probably things would become too complicated to make this useful.
In four dimensions, you can try the same approach. You have $\binom42=6$ coordinate planes, and for each of them you can give a possible rotation matrix. If the signs of the coordinates in the equation of your hyperboloid are the same, you use a regular rotation, otherwise you use hyperbolic functions. So for example if you have $x^2+y^2+z^2-w^2=1$ then any matrix modifying the $w$ coordinate would use hyperbolic functions, while those that just involve two of $x,y,z$ are regular rotations using circular trigonometric functions.