Method for finding eigenplanes of a linear transformation. Every introductory linear algebra course teaches methods for finding eigenvalues and associated eigenvectors of linear transformations $T$ acting on $\mathbb{R}^n$ and describes the geometric interpretation of such objects. My question pertains to the case of higher dimensional invariant subspaces - in particular, the case where one or more eigenvalues are complex.
For the sake of simplicity, assume $T$ is an invertible linear transformation on $\mathbb{R}^n$ with associated matrix $A_T\in Mat_{n}(\mathbb{R})$. Further, assume that the determinant $\mbox{det}(A_T)=1$. Let the linear subspace $V\subset\mathbb{R}^n$ be called an eigenplane of the transformation $T$ if $V$ is invariant under the action of $T$. That is, $T(V)=V$. Trivial eigenplanes of any transformation $T$ include the subspace $\{0\}$ and the image $\mbox{Im}\: T$ (Obviously in the case $T$ is invertible, $\mbox{Im}\: T=\mathbb{R}^n$).
Clearly, in the case of $\mbox{dim}(V)=1$, $V$ is just the subspace spanned by an eigenvector and, given a set of $k$, real-valued, eigenvectors $A=\{v_1,\ldots ,v_k\}$, the space $\mbox{span}(A)$ is a $k$-dimensional eigenplane. It follows that in the case that $T\colon\mathbb{R}^n\rightarrow\mathbb{R}^n$ has $n$ distinct, real eigenvectors, we can classify all eigenplanes of $T$ as the spans of all $2^n$ possible combinations of eigenvectors.
Now, I wonder if a similar classification of all eigenplanes of a transformation $T$ can be done when $T$ doesn't admit such a nice set of eigenvectors. In particular, how would an eigenvector of $T$ with complex entries be interpretted in such a setting? It's possible that I'm missing some very obvious linear algebra which makes finding eigenplanes simple, but in any case I would appreciate any input. I would also like to stress the generality of the dimension $n$. The question seems rather easy in low dimensions ($n\leq 3$).
 A: You're not missing anything; this is something that traditional linear algebra frankly fails at describing in an integrated, sensible way.
However, it can be described quite well using geometric algebra, which builds upon exterior algebra and gives great insights into the theory of linear operators.
Key to geometric algebra with linear operators is the wedge product of vectors.  It has the following properties:
$$a \wedge b= -b \wedge a \implies a \wedge a =0$$
This is similar to the cross product, but the result is interpreted as a bivector, an oriented plane.  Just as a vector has a linear (1d) subspace that it is associated with, a bivector has a 2d subspace that it is associated with.  This makes bivectors the natural objects for dealing with planes in any space (more general than 3d with normal vectors).
We define the action of a linear operator $\underline T$ on wedge products as follows:
$$\underline T(a \wedge b) \equiv \underline T(a) \wedge \underline T(b)$$
The wedge product is associative, so given a set of orthonormal basis vectors $\{e_1, e_2, \ldots, e_n\} \in \mathbb R^n$, the wedge product of these vectors determines a pseudoscalar:
$$i \equiv e_1 \wedge e_2 \wedge \ldots \wedge e_n$$
The use of $i$ to denote this object is meant to be suggestive.  In 2d, for instance, it has many of the mathematical properties of the complex imaginary, but without an actual complexification of the space.
The unit pseudoscalar is the unit object associated with the subspace of hypervolumes, a 1d vector space.  All other pseudoscalar quantities are scalar multiples of the pseudoscalar.  In particular, this means that
$$\underline T(i) = \alpha i$$
for some scalar $\alpha$.  This is a definition of the determinant.  $\alpha = \det \underline T$.
Now then, let's consider explicitly a rotation operator in 3d.  This operator will be called $\underline R$:
$$\begin{align*} \underline R(e_1) &= e_1 \cos \theta - e_2 \sin \theta \\ \underline R(e_2) &= e_2 \cos \theta + e_1 \sin \theta \\ \underline R(e_3) &= e_3 \end{align*}$$
As is usual, we can find $\underline R- \lambda\underline I$:
$$\begin{align*} (\underline R - \lambda \underline I)(e_1) &= e_1 (\cos \theta-\lambda) - e_2 \sin \theta \\ (\underline R - \lambda \underline I)(e_2) &= e_1 \sin \theta + e_2 (\cos \theta - \lambda) \\ (\underline R - \lambda \underline I)(e_3) &= (1- \lambda)e_3 \end{align*}$$
And we can read off the determinant by wedging these vectors.
$$(\underline R - \lambda \underline I)(i) = [(\cos \theta - \lambda)^2 (1-\lambda) + \sin^2 \theta (1-\lambda)] i$$
And we would conclude that $1 - \lambda = 0$ or $1- 2 \lambda \cos \theta + \lambda^2 = 0$.  The second equation usually has no real solutions, and traditionally, we would say that there are complex eigenvalues, and so on and so forth.
Now, actually backing it up to say, okay, what planes are actually eigenplanes?  That's a bit more difficult.  The characteristic polynomial has long been the main weapon in finding eigenvectors, but when it comes to eigenmultivectors, it's a much less useful tool because, even if you know that $\lambda$ cannot be a real scalar, it's not clear what it actually must be (just saying it's a complex number doesn't help, as there are lots of objects with the same algebraic properties in geometric algebra).
An approach that may work is to replace $\lambda$ by a generalized object $\lambda + \omega$, where $\omega$ is a bivector, not necessarily commuting with vectors.  Then, $\omega$ would clearly be $e_1 \wedge e_2 \sin \theta$, as this is what corresponds to the imaginary part of one of the eigenvalues.  But how one would solve this in a more general setting is unclear to me.  I have yet to see a reference on the subject, at least, and it's a problem of great importance in, for example, moment of inertia problems.
A: Since you give as an example eigenspaces, I assume you are talking about theoretical rather than computational eigenplanes. So assume you have factored the characteristic polynomial into a product of linear and irreducible quadratic factors. (This will also work for any field as you can factor the characteristic polynomial into irreducible factors but they may have degree > 2.) For an invariant $k$-dimensional subspace $V$, the matrix for $T$ is $$\left[\begin{array}{} A & B \\ 0 & D \end{array} \right]$$ where the first $k$ columns form a basis for $V$. The block diagonal form shows that the characteristic polynomial for $A$ is a factor of the characteristic polynomial for $T$. Now collect the invariant subspaces for all powers of each irreducible  factor and put the matrix of $T$ in Jordan canonical form over ${\mathbb C}$ using $2\times 2$ companion matrix blocks rather than constants on the diagonal. That this will give all invariant subspaces computationally follows from a proof of Basis Theorem for Finitely Generated Modules over a P.I.D. and the Cayley Hamilton Theorem that a matrix (or linear transformation) satisfies its characteristic polynomial provided you can factor the characteristic polynomial.
A: Let $V$ be an $n$-dimensional vector space over a field $F$, and let $T:V\to V$ be a linear map. There is actually a very nice way of finding invariant $k$-dimensional subspaces of $V$. It is based on the structure theorem for finitely generated modules over a PID. 
The vector space can be thought of as a $F[x]$-module with scalars acting on $V$ by scalar mutliplication, and $x$ acting on $V$ as $x\cdot v = T(v)$. Since $V$ is finite dimensional, this module is finitely generated. It is also a torison module; for example, the Cayley-Hamilton theorem says that the characteristic polynomial of $T$ annihilates $V$. 
By the structure theorem, we have find invariant factors for $V$ as a $F[x]$-module: there are unique monic polynomials $a_1(x),\ldots, a_t(x)$ such that $a_i|a_{i+1}$, and 
$$V \simeq F[x]/(a_1(x))\oplus\cdots\oplus F[x]/(a_t(x)).$$ 
Note that each component of this sum is fixed by $x$, i.e. it is invariant under $T$, and the dimension of the $i$th omponent is the degree of the polynomial $a_i$. Then the $k$-dimensional subspaces are those sums 
$$F[x]/(a_{i_1})\oplus \cdots \oplus F[x]/(a_{i_s})$$
where $\deg(a_{i_1})+\cdots+\deg(a_{i_s}) = k.$
EDIT:
Also, you can explicitly find the invariant factors by computing the Smtih normal form of $xI - A$, which is an $n\times n$ matrix with coefficients in $F[x]$, where $A$ is the matrix of the linear transformation $T$.
A: Let $T$ be a linear transformation acting on $\mathbb{R}^n$ with complex eigenvalues. Take two of the complex eigenvectors, $v_1$ and $v_2$ such that some linear combination of those, $\alpha v_1 + \beta v_2$, is a real vector. (The easiest vectors to pick are complex conjugates of eachother.) Then $T$ applied to this real vector will produce another real vector, however since $v_1$ and $v_2$ are eigenvectors, the vector $T(\alpha v_1 + \beta v_2)$ is still a linear combination of $v_1$ and $v_2$. Therefore, $\alpha v_1 + \beta v_2$ and $T(\alpha v_1 + \beta v_2)$ both lie within a two-dimensional subspace of $\mathbb{R}^n$ that is generated by $v_1$ and $v_2$, and is therefore invariant under $T$. If $T(\alpha v_1 + \beta v_2)$ is linearly independent of $\alpha v_1 + \beta v_2$, then those two vectors also generate the subspace. In this way, you can find two-dimensional real invariant subspaces of $\mathbb{R}^n$ for transformations with complex eigenvalues. You can generalize this to higher dimensional subspaces similarly to the real eigenvalue case. 
