I'll describe a computational method, assuming the rotation has been translated to be around an axis through the origin.
If you know the axis of rotation $A=(a,b,c)$, then you can find a vector orthogonal to this one. For example, if $a \ne 0$, $b \ne 0$, then $V=(b,-a,0)$ is such a vector. Or, if one component is $0$--say, for example $c=0$, then $V=(0,0,1)$ is orthogonal. And that covers all cases.
A third vector in a right-handed triple is found using the vector cross product: $W=A\times V$. Normalize $V$ and $W$ to unit vectors $\hat{V}$, $\hat{W}$.
Now, feed $\hat{V}$ into your matrix and determine the output in terms of $\hat{V}$, $\hat{W}$ using dot products
$$
R\hat{V} = \{(R\hat{V})\cdot\hat{V}\}\hat{V}+\{(R\hat{V})\cdot\hat{W}\}\hat{W}=\alpha\hat{V}+\beta\hat{W},
$$
and compare to the action of a right-handed rotation about $A$ through an angle $\theta$:
$$
\hat{V} \mapsto \cos\theta \hat{V}+\sin\theta \hat{W}.
$$
I will assume a real orthogonal matrix is involved. Further, for an orthogonal matrix to represent a "rotation" means that the determinant is 1.
1) Any (nonzero) multiple of an eigenvector is again an eigenvector, so it is not the case that eigenvectors of an orthogonal matrix must be unit vectors. However by the same token, any eigenvector can be scaled to be a vector of length one.
2) It is not generally the case that eigenvectors of a rotation matrix are mutually orthogonal, as a simple case illustrates. The trivial rotation corresponds to the identity matrix $I$, for which all nonzero vectors are eigenvectors corresponding to eigenvalue 1. Clearly we can pick two eigenvectors in this case which are not orthogonal.
More generally we can have an eigenspace for eigenvalue 1 of dimension greater than 1, and/or an eigenspace for eigenvalue -1 of even dimension greater than 0. In either of these cases we can also pick a pair of eigenvectors for the same eigenvalue that are not mutually orthogonal.
Typically a rotation matrix will have some complex eigenvalues, as for example in dimension two:
$$ \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} $$
whose characteristic polynomial is $\lambda^2 + 1$. In such cases there are no real eigenvectors corresponding to the complex eigenvalues.
A rotation preserves the length of a vector, so if real eigenvalues for a rotation matrix exist, they must have absolute value 1. Thus the only possible real eigenvalues are $\pm 1$.
Added:
The PDF linked from Comments on the Question describes $3\times 3$ orthogonal (real) matrices $A$, with special emphasis on rotations, i.e. where $\det A = 1$. Euler's Displacement Thm. is stated, to the effect that if the rotation is nonzero, there exists a unique axis of rotation, and the matrix $A$ itself can be expressed in terms of a unit vector $w$ directed along that axis and an angle of rotation $\phi$ (Rodriguez' Formula).
The immediate Question concerns conditions for which eigenvectors of $A$ will be orthogonal. Although the linked paper describes a real orthonormal basis for $\mathbb{R}^3$ in terms of unit eigenvectors, only one of the basis vectors is an eigenvector. The other two are linear combinations of complex eigenvectors.
Specifically the characteristic polynomial of $A$ will be of degree 3, and if the angle of rotation $\phi \neq 0$, then the three eigenvalues (roots of the characteristic polynomial) will be as follows: $\lambda_1 = 1, \lambda_2 = \cos \phi + i \sin \phi, \lambda_3 = \cos \phi - i \sin \phi$.
Given the unit vector $w$ as above directed along the axis of rotation, this is an eigenvector corresponding to eigenvalue $\lambda_1 = 1$:
$$ A w = w $$
Clearly the unique eigenvalue $\lambda_1 = 1$, being real, is distinct from the two complex (non-real) eigenvalues. Unless nonzero angle $\phi = \pi$ (which is the case Matt describes in one of his Comments on the Question), the eigenvalues $\lambda_2,\lambda_3$ are distinct (they are complex conjugates in any case).
Let $u,v$ be complex unit eigenvectors of $A$ corresponding to $\lambda_2,\lambda_3$ respectively. We first show that $w$ is orthogonal to both $u,v$.
Since $A^T=A^{-1}$ by definition of an orthogonal matrix, we have:
$$ w^T A = (A^T w)^T = (A^{-1} w)^T = w^T $$
That is $w^T$ is a left eigenvector corresponding to eigenvalue $\lambda_1 = 1$, just as $w$ is a right eigenvector for $\lambda_1$. This implies:
$$ \lambda_1 w^T u = w^T A u = \lambda_2 w^T u $$
Since $\lambda_2 \neq \lambda_1$, this tells us $w^T u = 0$. Note that this last is a scalar value, and regardless of the fact that $u$ is complex, it means that $w,u$ are orthogonal. Similarly $w,v$ are orthogonal.
Because $u,v$ are complex vectors, we want to clarify what it means for them to be "orthogonal".
Rather than defining the dot-product as simply $u^T v$, we need to take the complex conjugate of one vector in this product. It doesn't make a significant difference which one is conjugated; the effect of switching to conjugating the other only means the scalar dot-product is conjugated.
Thus $\overline{u}^T v = 0$ if and only if $u^T \overline{v} = 0$, so the notion of "orthgonality" is exactly the same, whichever vector in the definition gets conjugated.
Bearing in mind that conjugation does not alter a real value (matrix or vector), we can compute as follows:
$$ \overline{u}^T A = \overline{u^T A} = \overline{A^T u}^T = \overline{A^{-1} u}^T $$
Since $Au = \lambda_2 u$, we have $A^{-1} u = \lambda_2^{-1} u$. Also since $|\lambda_2|=1$, $\lambda_2^{-1} = \overline{\lambda_2}$. Putting these observations together:
$$ \overline{u}^T A = \overline{A^{-1} u}^T = \overline{\lambda_2^{-1} u^T}
= \overline{\lambda_2^{-1}} \overline{u}^T = \lambda_2 \overline{u}^T $$
We can now compute the orthgonality of $u,v$ similarly to previously for $w,u$ and $w,v$:
$$ \lambda_2 \overline{u}^T v = \overline{u}^T A v = \overline{u}^T (\lambda_3 v) = \lambda_3 \overline{u}^T v $$
That is, since $\lambda_2 \neq \lambda_3$, the above implies $\overline{u}^T v = 0$.
Therefore $\{u,v,w\}$ form an orthonormal basis that diagonalizes $A$, i.e. the unitarily similar matrix is diagonal:
$$ \begin{pmatrix} \lambda_2 & 0 & 0 \\ 0 & \lambda_3 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$
However it is more convenient to work with an orthonormal basis of real vectors (the paper outlines how to do this) $\{c_2,c_3,w\}$ with respect to which the representation is only "block" diagonal:
$$ \begin{pmatrix} \cos \phi & \sin \phi & 0 \\ -\sin \phi & \cos \phi & 0 \\ 0 & 0 & 1
\end{pmatrix} $$
Bear in mind that although this is an orthonormal basis, the first two real basis elements are not eigenvectors, but rather linear combinations of $u,v$ chosen to give us real vectors. The third basis element $w$ is the same real unit vector, in the direction of the axis of rotation, as before.
Best Answer
From $U U^T=U^T U=1$ you deduce that any eigenvalue has modulus 1. If $Ux=\lambda x$ then also $U\bar{x}=\bar{\lambda} \bar{x}$ so they come in complex pairs. In dimension 3, at least one eval has to be real and equal to +1 when looking at $SO$.
Suppose $U (x+ iy) = e^{i\theta}(x+ iy) = (\cos \theta +i \sin \theta) (x +iy) \ $ (with $\theta \neq 0$ mod $\pi$). Separating into real and imaginary parts we may write this as: $$ U [x \ \ y] = [x \ \ y] \left(\begin{matrix} \cos \theta & - \sin \theta\\ \sin \theta & \cos \theta \end{matrix} \right) $$
We claim that $x$ and $y$ are orthogonal and have the same length.
First notice that $U^T=U^{-1}$ so we also have $U^T (x+iy) = e^{-i\theta} (x+iy)$. From this we see that $$ Ux= x \cos \theta - y \sin \theta \ \ \mbox{ and } \ \ U^T x = x \cos \theta + y \sin \theta.$$ Take e.g. the scalar product with $x$ to obtain: $$ x^T U^Tx = (Ux)^T x = x^T U x \ \ \Rightarrow x^T y \sin \theta =0.$$ The assumption $\theta\neq 0$ mod $\pi$ implies that $x$ and $y$ are orthogonal. Calculating $|Ux|^2=|x|^2$ you deduce similarly that $|x|=|y|$.
This gives a very explicit geometric picture of rotation associated with the complex eigenvalues: The vectors $x$ and $y$ form an orthonormal base for a plane which is $U$ invariant and $U$ performs a rotation by angle $\theta$ in this plane.
Finally, you may notice that the discussion of complex eigenvalues applies in any dimension $\geq 2$.