Eigenvectors of a 4D rotation, and their interpretation Let us define a 4D rotation by using two unit quaternions: $$\mathring{q}_l=\frac{a+ib+jc+kd}{\left|a+ib+jc+kd\right|}$$ and $$\mathring{q}_r=\frac{e+ib+jc+kd}{\left|e+ib+jc+kd\right|}.$$ They differ only by the real term (the denominators are only for normalization). And let us consider two cases: the $a,b,c,d,e \in \mathbb{R}$ in the first case and $a,b,c,d,e \in \mathbb{C}$ in the second case - that would be biquaternions. Now let us consider a rotation of an arbitrary (non-unit) quaternion $\mathring{q}$ so that: 
$$\mathring{q}_{rotated}=\mathring{q}_l\,\mathring{q}\,\mathring{q}_r^* $$ 
(note the conjugate of $\mathring{q}_r$). I can represent this 4D rotation using a $4 \times 4$ rotation matrix, and thus I can find eigenvalues and eigenvectors of this matrix. The question is - what is the interpretation of eigenvectors of such rotation?
In 3D rotations the invariant vector (a vector that is not being rotated) is the rotation axis, and the eigenvalue has to be 1 (because it is a rotation, so no particular interpretation here). In 4D rotations I know that there  are two planes around which the rotation occurs, but are those planes the invariant of the rotation? If yes, then how can I describe this plane as an eigenvector quaternion - would such quaternion be a "normal vector" of this plane, analogous to 3D space where normal vector of a plane is just a 3D vector?
How is that different between quaternions and biquaternions?
Please note that since $e \neq a$ the normalized quaternion $\mathring{q}_r$ will be very different from $\mathring{q}_l$ - all the imaginary components of unit quaternion (those created by $b,c,d$) will be in fact different, thanks to normalization.
 A: To be honest, I have a hard time interpreting geometrically what's going on here when using quaternions or biquaternions or anything else.  All the algebra of rotations in 4d is adequately handled by a geometric algebra, with the elements of that algebra having clear geometric interpretations.  The mathematics is similar to quaternions, but differs in some conceptual ways.

Okay, Muphrid, what is geometric algebra, and how can it help us talk about rotations?
Geometric algebra is a kind of clifford algebra.  It posits a "geometric product" between vectors that is denoted by juxtaposition, so the geometric product of two vectors $a$ and $b$ is denoted $ab$.  This product has the following properties:
$$aa = |a|^2, \quad (ab)c = a(bc)$$
From these two properties, you get a wealth of useful structure to supplement traditional vector algebra.  Most relevant here are bivectors, which represent oriented planes.  Given four orthonormal basis vectors $e_1, e_2, e_3, e_4$, you get the following unit bivectors:
$$\text{Bivectors:} \, e_1 e_2, e_2 e_3, e_3 e_1, e_1 e_4, e_2 e_4, e_3 e_4$$
The geometric product of vectors also produces objects called rotors, which are analogues of quaternions in that they perform rotations.  For example, in 3d you can multiply two vectors $a$ and $b$ to get the following:
$$a b = (a^1 b^1 + a^2 b^2 + a^3 b^3) + (a^2 b^3 - a^3 b^2) e_2 e_3 + (a^3 b^1 - a^1 b^3) e_3 e_1 + (a^1 b^2 - a^2 b^1) e_1 e_2$$
This has four terms, like a quaternion.  In fact, you can make the following identifications:
$$i = -e_2 e_3, \quad j = -e_3 e_1, \quad k = - e_1 e_2$$
One advantage that geometric algebra has over quaternions is that quaternions have to perform double-duty: you use pure imaginary quats to represent vectors.  GA doesn't do this; vectors and rotors are kept clearly separated according to their geometric properties and function.  You would never confuse $e_1$--a vector--with $e_1 e_2$--a bivector.

But Muphrid, what about 4d?  Isn't that what we're interested in here?
Right, let's talk about the geometric algebra of 4d Euclidean space.  As I said, there are six unit bivectors, and you might have noticed that there were six imaginaries involved when you talk about two quaternions.  Not a coincidence.  GA lets us handle that directly, rather than hacking quaternions to make it all work.
Here's how:  there's an important concept of duality, which we represent though multiplication by the pseudoscalar, which I will call $\epsilon = e_1 e_2 e_3 e_4$.  The pseudoscalar multiplied by a bivector returns the corresponding orthogonal bivector.  Let's see how:
$$\epsilon e_1 e_2 = - e_3 e_4, \quad \epsilon e_2 e_3 = - e_1 e_4, \quad \epsilon e_3 e_1 = - e_2 e_4$$
This is the reason why you can get away with using two quaternions or biquaternions: any bivector can be expressed as a linear combination like so:
$$B = (\alpha e_1 e_2 + \beta e_2 e_3 + \gamma e_3 e_1) + \epsilon (\lambda e_1 e_2 + \mu e_2 e_3 + \nu e_3 e_1)$$
Again, $\epsilon = e_1 e_2 e_3 e_4$, and $\epsilon \epsilon = -1$, as it happens.  It acts like yet another imaginary unit, partitioning the very complicated rotors of 4d into 2 3d-like rotors.

So Muphrid, what does that tell us about eigenvectors (or eigenbivectors, or eigenrotors) of a general rotation operation in 4d?
These are most easily understood using the double-isoclinic decomposition.
Let a general rotation bivector be given as $B = U + \epsilon V$, where $U, V$ are linear combinations of $e_1 e_2, e_2 e_3, e_3 e_1$.  Let $I_\pm = (1 \pm \epsilon)/2$, then we can rewrite $B$ as
$$B = I_+ X + I_- Y, \quad X = (U+V)/2, \quad Y = (U-V)/2$$
The rotations using $I_+ X$ and $I_- Y$ are both "isoclinic", meaning that they each rotate in two orthogonal planes by the same angle.  The corresponding rotation takes the form
$$\underline R(a) = I_+ \exp(X) a \exp(-Y) + I_- \exp(Y) a \exp(-X)$$
Outside of special cases where $X$ and $Y$ are linearly dependent, I can't see any individual vectors that would, in general, be eigenvectors.
For eigenbivectors, the power series of an exponential $\exp(B)$ tells us that $B$ and $\epsilon B$ are both eigenbivectors, which is a much easier analysis.
