Rotation interpolation

Quaternions

Let $p_1$ and $p_2$ be unit quaternions. The slerp formula for quaternions is $$p(t,p_1,p_2) = p_1^{(1-t)}p_2^t = p_1(p_1^{-1}p_2)^t = p_1\exp(t\underbrace{\log(p_1^{-1}p_2)}_{\theta v}) =p_1(\cos(t\theta) +\sin(t\theta)v).$$ Here we have used that $p_1^{-1}p_2$ is a unit quaternion and Euler's formula for unit quaternions with zero real part: $$p_1^{-1}p_2 = \exp(\theta v) = \cos(\theta) + \sin(\theta) v, \quad \|v\|=1, \quad v_w=0.$$ We can express $v$ in terms of $p_2$ from $p_2 = p_1\exp(\theta v)$ as follows: $$p_2 = p_1\exp(\theta v) = p_1(\cos(\theta)+\sin(\theta)v) \implies v = p_1^{-1}\frac{p_2-p_1\cos(\theta)}{\sin\theta}.$$ Plugging this into our expression for the slerp yields the familiar formula \begin{align}p(t,p_1,p_2) &= p_1\cos(t\theta) + \sin(t\theta)\frac{p_2-p_1\cos(\theta)}{\sin\theta} \\ &= \frac{\sin(\theta)\cos(t\theta)-\sin(t\theta)\cos(\theta)}{\sin\theta}p_1 + \frac{\sin(t\theta)}{\sin\theta)}p_2 \\ &= \frac{\sin((1-t)\theta)}{\sin\theta}p_1 + \frac{\sin(t\theta)}{\sin\theta}p_2. \end{align}

So you are indeed interpolating the angle $\theta$ but with axis $v = \frac{\log(p_1^{-1}p_2)}{\theta}$.

Matrices

If you are using matrices $R_1$ and $R_2$ to represent your rotations then the formula is almost the same $$R(t,R_1,R_2) = R_2^{t}R_1^{(1-t)}.$$ The arguments are switched because rotation with matrices is applied as $u' = Ru$ so we compose matrices on the left, while we compose quaternions on the right.

Similar to the quaternion case we can look at $$R(t,R_1,R_2) = R_2^{t}R_1^{(1-t)} = \exp(t\log(R_2R_1^{-1}))R_1 = \exp(t\theta [v]_{\times})R_1,$$ where the matrix $[v]_{\times}$ built from the unit length vector $v$ looks like this $$[v]_{\times} = \begin{bmatrix} 0 & -v_z & v_y \\ v_z & 0 & -v_x \\ -v_y & v_x & 0 \end{bmatrix}.$$
By using the series definition of the exponential and properties of the above matrix one can show that the exponential of the matrix $t\theta[v]_{\times}$ is equal to $$\exp(t\theta [v]_{\times}) = I+\sin(t\theta) [v]_{\times} + (1-\cos(t\theta)) [v]_{\times}^2.$$ Vice versa, suppose you have the matrix $P = R_2R_1^{-1}$ and you want to find the logarithm. Then use that $P-P^{\top} = 2\sin(t\theta) [v]_{\times}$ (because $[v]_{\times}$ is anti-symmetric). Then you have that $$Q = \frac{P-P^{\top}}{2} = \sin(t\theta) [v]_{\times}.$$ Then by matching the matrices and using that $\|v\|=1$ you can figure out that $\sin\theta = \pm\sqrt{q_{23}^2 + q_{13}^2 + q_{21}^2}$ and that $v_x = q_{21}/\sin\theta$, $v_y = q_{13}/\sin\theta$, $v_z=q_{23}/\sin\theta$. For more details see the article by Eberly "Interpolating Rigid Motions in 3D".

Higher Dimensions

You can find formulae for the exponential and logarithm for some higher dimensions in the paper by Kaiser "Solving the matrix exponential function for special orthogonal groups $SO(n)$ and the exceptional $G_2$", but it gets rather cumbersome. It's similarly cumbersome if you use (sums of) bivectors in higher dimensions as a generalization to the quaternion $\theta v$. In fact I think it's even more cumbersome than matrices in 5D+, we are just (un)lucky that in 3D quaternion interpolation is cheaper than matrix interpolation (however rotation matrix application is cheaper than quaternion rotation application in terms of number of ops).